Actual source code: ex6f.F90
slepc-main 2024-11-09
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 MatIsing_Mult
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,MatIsing_Mult,ierr))
75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: ! Create the eigensolver and display info
77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: ! ** Create eigensolver context
80: PetscCallA(EPSCreate(PETSC_COMM_WORLD,eps,ierr))
82: ! ** Set operators. In this case, it is a standard eigenvalue problem
83: PetscCallA(EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr))
84: PetscCallA(EPSSetProblemType(eps,EPS_NHEP,ierr))
86: ! ** Set solver parameters at runtime
87: PetscCallA(EPSSetFromOptions(eps,ierr))
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: ! Solve the eigensystem
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: PetscCallA(EPSSolve(eps,ierr))
94: PetscCallA(EPSGetIterationNumber(eps,its,ierr))
95: if (rank .eq. 0) then
96: write(*,'(A,I4)') ' Number of iterations of the method: ',its
97: endif
99: ! ** Optional: Get some information from the solver and display it
100: PetscCallA(EPSGetType(eps,tname,ierr))
101: if (rank .eq. 0) then
102: write(*,'(A,A)') ' Solution method: ', tname
103: endif
104: PetscCallA(EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
105: if (rank .eq. 0) then
106: write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
107: endif
108: PetscCallA(EPSGetTolerances(eps,tol,maxit,ierr))
109: if (rank .eq. 0) then
110: write(*,'(A,1PE11.4,A,I6)') ' Stopping condition: tol=',tol,', maxit=', maxit
111: endif
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: ! Display solution and clean up
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: ! ** show detailed info unless -terse option is given by user
118: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
119: if (terse) then
120: PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
121: else
122: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
123: PetscCallA(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr))
124: PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
125: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
126: endif
127: PetscCallA(EPSDestroy(eps,ierr))
128: PetscCallA(MatDestroy(A,ierr))
129: PetscCallA(SlepcFinalize(ierr))
130: end
132: ! -------------------------------------------------------------------
133: !
134: ! MatIsing_Mult - user provided matrix-vector multiply
135: !
136: ! Input Parameters:
137: ! A - matrix
138: ! x - input vector
139: !
140: ! Output Parameter:
141: ! y - output vector
142: !
143: subroutine MatIsing_Mult(A,x,y,ierr)
144: #include <petsc/finclude/petscmat.h>
145: use petscmat
146: implicit none
148: Mat A
149: Vec x,y
150: PetscInt trans,one,N
151: PetscScalar x_array(1),y_array(1)
152: PetscOffset i_x,i_y
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,x_array,i_x,ierr))
160: PetscCall(VecGetArray(y,y_array,i_y,ierr))
162: trans = 0
163: one = 1
164: call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)
166: PetscCall(VecRestoreArrayRead(x,x_array,i_x,ierr))
167: PetscCall(VecRestoreArray(y,y_array,i_y,ierr))
169: end
171: ! -------------------------------------------------------------------
172: ! The actual routine for the matrix-vector product
173: ! See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
175: SUBROUTINE MVMISG(TRANS, N, M, X, LDX, Y, LDY)
176: #include <petsc/finclude/petscsys.h>
177: use petscsys
178: ! ..
179: ! .. Scalar Arguments ..
180: PetscInt LDY, LDX, M, N, TRANS
181: ! ..
182: ! .. Array Arguments ..
183: PetscScalar Y(LDY, *), X(LDX, *)
184: ! ..
185: !
186: ! Purpose
187: ! =======
188: !
189: ! Compute
190: !
191: ! Y(:,1:M) = op(A)*X(:,1:M)
192: !
193: ! where op(A) is A or A' (the transpose of A). The A is the Ising
194: ! matrix.
195: !
196: ! Arguments
197: ! =========
198: !
199: ! TRANS (input) INTEGER
200: ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
201: ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
202: !
203: ! N (input) INTEGER
204: ! The order of the matrix A. N has to be an even number.
205: !
206: ! M (input) INTEGER
207: ! The number of columns of X to multiply.
208: !
209: ! X (input) DOUBLE PRECISION array, dimension (LDX, M)
210: ! X contains the matrix (vectors) X.
211: !
212: ! LDX (input) INTEGER
213: ! The leading dimension of array X, LDX >= max(1, N)
214: !
215: ! Y (output) DOUBLE PRECISION array, dimension (LDX, M)
216: ! contains the product of the matrix op(A) with X.
217: !
218: ! LDY (input) INTEGER
219: ! The leading dimension of array Y, LDY >= max(1, N)
220: !
221: ! ===================================================================
222: !
223: !
225: ! .. Local Variables ..
226: PetscInt I, K
227: PetscReal ALPHA, BETA
228: PetscReal COSA, COSB, SINA, SINB
229: PetscScalar TEMP, TEMP1
230: !
231: ! .. Intrinsic functions ..
232: INTRINSIC COS, SIN
233: !
234: ALPHA = PETSC_PI/4
235: BETA = PETSC_PI/4
236: COSA = COS(ALPHA)
237: SINA = SIN(ALPHA)
238: COSB = COS(BETA)
239: SINB = SIN(BETA)
240: !
241: IF (TRANS.EQ.0) THEN
242: !
243: ! Compute Y(:,1:M) = A*X(:,1:M)
245: DO 30 K = 1, M
246: !
247: Y(1, K) = COSB*X(1, K) - SINB*X(N, K)
248: DO 10 I = 2, N-1, 2
249: Y(I, K) = COSB*X(I, K) + SINB*X(I+1, K)
250: Y(I+1, K) = -SINB*X(I, K) + COSB*X(I+1, K)
251: 10 CONTINUE
252: Y(N, K) = SINB*X(1, K) + COSB*X(N, K)
253: !
254: DO 20 I = 1, N, 2
255: TEMP = COSA*Y(I, K) + SINA*Y(I+1, K)
256: Y(I+1, K) = -SINA*Y(I, K) + COSA*Y(I+1, K)
257: Y(I, K) = TEMP
258: 20 CONTINUE
259: !
260: 30 CONTINUE
261: !
262: ELSE IF (TRANS.EQ.1) THEN
263: !
264: ! Compute Y(:1:M) = A'*X(:,1:M)
265: !
266: DO 60 K = 1, M
267: !
268: DO 40 I = 1, N, 2
269: Y(I, K) = COSA*X(I, K) - SINA*X(I+1, K)
270: Y(I+1, K) = SINA*X(I, K) + COSA*X(I+1, K)
271: 40 CONTINUE
272: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
273: DO 50 I = 2, N-1, 2
274: TEMP1 = COSB*Y(I, K) - SINB*Y(I+1, K)
275: Y(I+1, K) = SINB*Y(I, K) + COSB*Y(I+1, K)
276: Y(I, K) = TEMP1
277: 50 CONTINUE
278: Y(N, K) = -SINB*Y(1, K) + COSB*Y(N, K)
279: Y(1, K) = TEMP
280: !
281: 60 CONTINUE
282: !
283: END IF
284: !
285: RETURN
286: !
287: ! END OF MVMISG
288: END
290: !/*TEST
291: !
292: ! test:
293: ! suffix: 1
294: ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
295: ! requires: !complex
296: !
297: !TEST*/