Actual source code: ex6f.F90
slepc-3.23.3 2025-09-08
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5: !
6: ! This file is part of SLEPc.
7: ! SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: !
10: ! Program usage: mpiexec -n <np> ./ex6f [-help] [-m <m>] [all SLEPc options]
11: !
12: ! Description: Eigensystem from the Ising model for ferromagnetic materials.
13: ! Information about the model can be found at the following
14: ! site https://math.nist.gov/MatrixMarket/data/NEP
15: !
16: ! The command line options are:
17: ! -m <m>, where <m> is the number of 2x2 blocks, i.e. matrix size N=2*m
18: !
19: ! ----------------------------------------------------------------------
20: !
21: program main
22: #include <slepc/finclude/slepceps.h>
23: use slepceps
24: implicit none
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: ! Declarations
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: !
30: ! Variables:
31: ! A operator matrix
32: ! eps eigenproblem solver context
34: Mat A
35: EPS eps
36: EPSType tname
37: PetscReal tol
38: PetscInt N, m
39: PetscInt nev, maxit, its
40: PetscMPIInt sz, rank
41: PetscErrorCode ierr
42: PetscBool flg, terse
44: ! This is the routine to use for matrix-free approach
45: !
46: external MatMult_Ising, MatMultTranspose_Ising
48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: ! Beginning of program
50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
53: #if defined(PETSC_USE_COMPLEX)
54: SETERRA(PETSC_COMM_SELF,PETSC_ERR_SUP,'This example requires real numbers')
55: #endif
56: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr))
57: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
58: if (sz .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only!'); endif
59: m = 30
60: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
61: N = 2*m
63: if (rank .eq. 0) then
64: write(*,'(/A,I6,A/)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
65: endif
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: ! Register the matrix-vector subroutine for the operator that defines
69: ! the eigensystem, Ax=kx
70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: PetscCallA(MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_INTEGER,A,ierr))
73: PetscCallA(MatShellSetOperation(A,MATOP_MULT,MatMult_Ising,ierr))
74: PetscCallA(MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,MatMultTranspose_Ising,ierr))
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: ! Create the eigensolver and display info
78: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: ! ** Create eigensolver context
81: PetscCallA(EPSCreate(PETSC_COMM_WORLD,eps,ierr))
83: ! ** Set operators. In this case, it is a standard eigenvalue problem
84: PetscCallA(EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr))
85: PetscCallA(EPSSetProblemType(eps,EPS_NHEP,ierr))
87: ! ** Set solver parameters at runtime
88: PetscCallA(EPSSetFromOptions(eps,ierr))
90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Solve the eigensystem
92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: PetscCallA(EPSSolve(eps,ierr))
95: PetscCallA(EPSGetIterationNumber(eps,its,ierr))
96: if (rank .eq. 0) then
97: write(*,'(A,I4)') ' Number of iterations of the method: ',its
98: endif
100: ! ** Optional: Get some information from the solver and display it
101: PetscCallA(EPSGetType(eps,tname,ierr))
102: if (rank .eq. 0) then
103: write(*,'(A,A)') ' Solution method: ', tname
104: endif
105: PetscCallA(EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
106: if (rank .eq. 0) then
107: write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
108: endif
109: PetscCallA(EPSGetTolerances(eps,tol,maxit,ierr))
110: if (rank .eq. 0) then
111: write(*,'(A,1PE11.4,A,I6)') ' Stopping condition: tol=',tol,', maxit=', maxit
112: endif
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: ! Display solution and clean up
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: ! ** show detailed info unless -terse option is given by user
119: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
120: if (terse) then
121: PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
122: else
123: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
124: PetscCallA(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr))
125: PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
126: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
127: endif
128: PetscCallA(EPSDestroy(eps,ierr))
129: PetscCallA(MatDestroy(A,ierr))
130: PetscCallA(SlepcFinalize(ierr))
131: end
133: ! -------------------------------------------------------------------
134: !
135: ! MatMult_Ising - user provided matrix-vector multiply
136: !
137: ! Input Parameters:
138: ! A - matrix
139: ! x - input vector
140: !
141: ! Output Parameter:
142: ! y - output vector
143: !
144: subroutine MatMult_Ising(A,x,y,ierr)
145: #include <petsc/finclude/petscmat.h>
146: use petscmat
147: implicit none
149: Mat A
150: Vec x,y
151: PetscInt trans,one,N
152: PetscScalar, pointer :: xx(:),yy(:)
153: PetscErrorCode ierr
155: ! The actual routine for the matrix-vector product
156: external mvmisg
158: PetscCall(MatGetSize(A,N,PETSC_NULL_INTEGER,ierr))
159: PetscCall(VecGetArrayRead(x,xx,ierr))
160: PetscCall(VecGetArray(y,yy,ierr))
162: trans = 0
163: one = 1
164: call mvmisg(trans,N,one,xx,N,yy,N)
166: PetscCall(VecRestoreArrayRead(x,xx,ierr))
167: PetscCall(VecRestoreArray(y,yy,ierr))
169: end
171: ! -------------------------------------------------------------------
172: !
173: ! MatMultTranspose_Ising - user provided transpose matrix-vector multiply
174: !
175: ! Input Parameters:
176: ! A - matrix
177: ! x - input vector
178: !
179: ! Output Parameter:
180: ! y - output vector
181: !
182: subroutine MatMultTranspose_Ising(A,x,y,ierr)
183: #include <petsc/finclude/petscmat.h>
184: use petscmat
185: implicit none
187: Mat A
188: Vec x,y
189: PetscInt trans,one,N
190: PetscScalar, pointer :: xx(:),yy(:)
191: PetscErrorCode ierr
193: ! The actual routine for the transpose matrix-vector product
194: external mvmisg
196: PetscCall(MatGetSize(A,N,PETSC_NULL_INTEGER,ierr))
197: PetscCall(VecGetArrayRead(x,xx,ierr))
198: PetscCall(VecGetArray(y,yy,ierr))
200: trans = 1
201: one = 1
202: call mvmisg(trans,N,one,xx,N,yy,N)
204: PetscCall(VecRestoreArrayRead(x,xx,ierr))
205: PetscCall(VecRestoreArray(y,yy,ierr))
207: end
209: ! -------------------------------------------------------------------
210: ! The actual routine for the matrix-vector product
211: ! See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
213: SUBROUTINE MVMISG(TRANS, N, M, X, LDX, Y, LDY)
214: #include <petsc/finclude/petscsys.h>
215: use petscsys
216: ! ..
217: ! .. Scalar Arguments ..
218: PetscInt LDY, LDX, M, N, TRANS
219: ! ..
220: ! .. Array Arguments ..
221: PetscScalar Y(LDY, *), X(LDX, *)
222: ! ..
223: !
224: ! Purpose
225: ! =======
226: !
227: ! Compute
228: !
229: ! Y(:,1:M) = op(A)*X(:,1:M)
230: !
231: ! where op(A) is A or A' (the transpose of A). The A is the Ising
232: ! matrix.
233: !
234: ! Arguments
235: ! =========
236: !
237: ! TRANS (input) INTEGER
238: ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
239: ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
240: !
241: ! N (input) INTEGER
242: ! The order of the matrix A. N has to be an even number.
243: !
244: ! M (input) INTEGER
245: ! The number of columns of X to multiply.
246: !
247: ! X (input) DOUBLE PRECISION array, dimension (LDX, M)
248: ! X contains the matrix (vectors) X.
249: !
250: ! LDX (input) INTEGER
251: ! The leading dimension of array X, LDX >= max(1, N)
252: !
253: ! Y (output) DOUBLE PRECISION array, dimension (LDX, M)
254: ! contains the product of the matrix op(A) with X.
255: !
256: ! LDY (input) INTEGER
257: ! The leading dimension of array Y, LDY >= max(1, N)
258: !
259: ! ===================================================================
260: !
261: !
263: ! .. Local Variables ..
264: PetscInt I, K
265: PetscReal ALPHA, BETA
266: PetscReal COSA, COSB, SINA, SINB
267: PetscScalar TEMP, TEMP1
268: !
269: ! .. Intrinsic functions ..
270: INTRINSIC COS, SIN
271: !
272: ALPHA = PETSC_PI/4
273: BETA = PETSC_PI/4
274: COSA = COS(ALPHA)
275: SINA = SIN(ALPHA)
276: COSB = COS(BETA)
277: SINB = SIN(BETA)
278: !
279: IF (TRANS.EQ.0) THEN
280: !
281: ! Compute Y(:,1:M) = A*X(:,1:M)
283: DO 30 K = 1, M
284: !
285: Y(1, K) = COSB*X(1, K) - SINB*X(N, K)
286: DO 10 I = 2, N-1, 2
287: Y(I, K) = COSB*X(I, K) + SINB*X(I+1, K)
288: Y(I+1, K) = -SINB*X(I, K) + COSB*X(I+1, K)
289: 10 CONTINUE
290: Y(N, K) = SINB*X(1, K) + COSB*X(N, K)
291: !
292: DO 20 I = 1, N, 2
293: TEMP = COSA*Y(I, K) + SINA*Y(I+1, K)
294: Y(I+1, K) = -SINA*Y(I, K) + COSA*Y(I+1, K)
295: Y(I, K) = TEMP
296: 20 CONTINUE
297: !
298: 30 CONTINUE
299: !
300: ELSE IF (TRANS.EQ.1) THEN
301: !
302: ! Compute Y(:1:M) = A'*X(:,1:M)
303: !
304: DO 60 K = 1, M
305: !
306: DO 40 I = 1, N, 2
307: Y(I, K) = COSA*X(I, K) - SINA*X(I+1, K)
308: Y(I+1, K) = SINA*X(I, K) + COSA*X(I+1, K)
309: 40 CONTINUE
310: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
311: DO 50 I = 2, N-1, 2
312: TEMP1 = COSB*Y(I, K) - SINB*Y(I+1, K)
313: Y(I+1, K) = SINB*Y(I, K) + COSB*Y(I+1, K)
314: Y(I, K) = TEMP1
315: 50 CONTINUE
316: Y(N, K) = -SINB*Y(1, K) + COSB*Y(N, K)
317: Y(1, K) = TEMP
318: !
319: 60 CONTINUE
320: !
321: END IF
322: !
323: RETURN
324: !
325: ! END OF MVMISG
326: END
328: !/*TEST
329: !
330: ! testset:
331: ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
332: ! requires: !complex
333: ! output_file: output/ex6f_1.out
334: ! filter: grep -v iterations | sed -e 's/-0.00000/0.00000/g'
335: ! test:
336: ! suffix: 1
337: ! test:
338: ! suffix: 1_ts
339: ! args: -eps_two_sided
340: ! requires: !single
341: !
342: !TEST*/