Actual source code: ex6f90.F90
slepc-3.20.0 2023-09-29
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: return
170: end
172: ! -------------------------------------------------------------------
173: ! The actual routine for the matrix-vector product
174: ! See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
176: SUBROUTINE MVMISG(TRANS, N, M, X, LDX, Y, LDY)
177: #include <petsc/finclude/petscsys.h>
178: use petscsys
179: ! ..
180: ! .. Scalar Arguments ..
181: PetscInt LDY, LDX, M, N, TRANS
182: ! ..
183: ! .. Array Arguments ..
184: PetscScalar Y(LDY, *), X(LDX, *)
185: ! ..
186: !
187: ! Purpose
188: ! =======
189: !
190: ! Compute
191: !
192: ! Y(:,1:M) = op(A)*X(:,1:M)
193: !
194: ! where op(A) is A or A' (the transpose of A). The A is the Ising
195: ! matrix.
196: !
197: ! Arguments
198: ! =========
199: !
200: ! TRANS (input) INTEGER
201: ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
202: ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
203: !
204: ! N (input) INTEGER
205: ! The order of the matrix A. N has to be an even number.
206: !
207: ! M (input) INTEGER
208: ! The number of columns of X to multiply.
209: !
210: ! X (input) DOUBLE PRECISION array, dimension (LDX, M)
211: ! X contains the matrix (vectors) X.
212: !
213: ! LDX (input) INTEGER
214: ! The leading dimension of array X, LDX >= max(1, N)
215: !
216: ! Y (output) DOUBLE PRECISION array, dimension (LDX, M)
217: ! contains the product of the matrix op(A) with X.
218: !
219: ! LDY (input) INTEGER
220: ! The leading dimension of array Y, LDY >= max(1, N)
221: !
222: ! ===================================================================
223: !
224: !
226: ! .. Local Variables ..
227: PetscInt I, K
228: PetscReal ALPHA, BETA
229: PetscReal COSA, COSB, SINA, SINB
230: PetscScalar TEMP, TEMP1
231: !
232: ! .. Intrinsic functions ..
233: INTRINSIC COS, SIN
234: !
235: ALPHA = PETSC_PI/4
236: BETA = PETSC_PI/4
237: COSA = COS(ALPHA)
238: SINA = SIN(ALPHA)
239: COSB = COS(BETA)
240: SINB = SIN(BETA)
241: !
242: IF (TRANS.EQ.0) THEN
243: !
244: ! Compute Y(:,1:M) = A*X(:,1:M)
246: DO 30 K = 1, M
247: !
248: Y(1, K) = COSB*X(1, K) - SINB*X(N, K)
249: DO 10 I = 2, N-1, 2
250: Y(I, K) = COSB*X(I, K) + SINB*X(I+1, K)
251: Y(I+1, K) = -SINB*X(I, K) + COSB*X(I+1, K)
252: 10 CONTINUE
253: Y(N, K) = SINB*X(1, K) + COSB*X(N, K)
254: !
255: DO 20 I = 1, N, 2
256: TEMP = COSA*Y(I, K) + SINA*Y(I+1, K)
257: Y(I+1, K) = -SINA*Y(I, K) + COSA*Y(I+1, K)
258: Y(I, K) = TEMP
259: 20 CONTINUE
260: !
261: 30 CONTINUE
262: !
263: ELSE IF (TRANS.EQ.1) THEN
264: !
265: ! Compute Y(:1:M) = A'*X(:,1:M)
266: !
267: DO 60 K = 1, M
268: !
269: DO 40 I = 1, N, 2
270: Y(I, K) = COSA*X(I, K) - SINA*X(I+1, K)
271: Y(I+1, K) = SINA*X(I, K) + COSA*X(I+1, K)
272: 40 CONTINUE
273: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
274: DO 50 I = 2, N-1, 2
275: TEMP1 = COSB*Y(I, K) - SINB*Y(I+1, K)
276: Y(I+1, K) = SINB*Y(I, K) + COSB*Y(I+1, K)
277: Y(I, K) = TEMP1
278: 50 CONTINUE
279: Y(N, K) = -SINB*Y(1, K) + COSB*Y(N, K)
280: Y(1, K) = TEMP
281: !
282: 60 CONTINUE
283: !
284: END IF
285: !
286: RETURN
287: !
288: ! END OF MVMISG
289: END
291: !/*TEST
292: !
293: ! test:
294: ! suffix: 1
295: ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
296: ! requires: !complex
297: !
298: !TEST*/