Actual source code: ex6f90.F90
slepc-3.11.2 2019-07-30
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2019, 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: This example solves the eigensystem arising in the Ising
13: ! model for ferromagnetic materials. The file mvmisg.f must be linked
14: ! together. Information about the model can be found at the following
15: ! site http://math.nist.gov/MatrixMarket/data/NEP
16: !
17: ! The command line options are:
18: ! -m <m>, where <m> is the number of 2x2 blocks, i.e. matrix size N=2*m
19: !
20: ! ----------------------------------------------------------------------
21: !
22: program main
23: #include <slepc/finclude/slepceps.h>
24: use slepceps
25: implicit none
27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
28: ! Declarations
29: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: !
31: ! Variables:
32: ! A operator matrix
33: ! eps eigenproblem solver context
35: Mat A
36: EPS eps
37: EPSType tname
38: PetscReal tol
39: PetscInt N, m
40: PetscInt nev, maxit, its
41: PetscMPIInt sz, rank
42: PetscErrorCode ierr
43: PetscBool flg, terse
45: ! This is the routine to use for matrix-free approach
46: !
47: external MatIsing_Mult
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Beginning of program
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
54: if (ierr .ne. 0) then
55: print*,'SlepcInitialize failed'
56: stop
57: endif
58: #if defined(PETSC_USE_COMPLEX)
59: call PetscError(PETSC_COMM_SELF,1,0,'This example requires real numbers')
60: call MPIU_Abort(PETSC_COMM_SELF,ierr)
61: #endif
62: call MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr);CHKERRA(ierr)
63: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
64: if (sz .ne. 1) then; SETERRA(PETSC_COMM_SELF,1,'This is a uniprocessor example only!'); endif
65: m = 30
66: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr);CHKERRA(ierr)
67: N = 2*m
69: if (rank .eq. 0) then
70: write(*,'(/A,I6,A/)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
71: endif
73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: ! Register the matrix-vector subroutine for the operator that defines
75: ! the eigensystem, Ax=kx
76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: call MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_INTEGER,A,ierr);CHKERRA(ierr)
79: call MatShellSetOperation(A,MATOP_MULT,MatIsing_Mult,ierr);CHKERRA(ierr)
81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: ! Create the eigensolver and display info
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: ! ** Create eigensolver context
86: call EPSCreate(PETSC_COMM_WORLD,eps,ierr);CHKERRA(ierr)
88: ! ** Set operators. In this case, it is a standard eigenvalue problem
89: call EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr);CHKERRA(ierr)
90: call EPSSetProblemType(eps,EPS_NHEP,ierr);CHKERRA(ierr)
92: ! ** Set solver parameters at runtime
93: call EPSSetFromOptions(eps,ierr);CHKERRA(ierr)
95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: ! Solve the eigensystem
97: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: call EPSSolve(eps,ierr);CHKERRA(ierr)
100: call EPSGetIterationNumber(eps,its,ierr);CHKERRA(ierr)
101: if (rank .eq. 0) then
102: write(*,'(A,I4)') ' Number of iterations of the method: ',its
103: endif
105: ! ** Optional: Get some information from the solver and display it
106: call EPSGetType(eps,tname,ierr);CHKERRA(ierr)
107: if (rank .eq. 0) then
108: write(*,'(A,A)') ' Solution method: ', tname
109: endif
110: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
111: if (rank .eq. 0) then
112: write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
113: endif
114: call EPSGetTolerances(eps,tol,maxit,ierr);CHKERRA(ierr)
115: if (rank .eq. 0) then
116: write(*,'(A,1PE11.4,A,I6)') ' Stopping condition: tol=',tol,', maxit=', maxit
117: endif
119: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: ! Display solution and clean up
121: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: ! ** show detailed info unless -terse option is given by user
124: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
125: if (terse) then
126: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
127: else
128: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
129: call EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
130: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
131: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
132: endif
133: call EPSDestroy(eps,ierr);CHKERRA(ierr)
134: call MatDestroy(A,ierr);CHKERRA(ierr)
135: call SlepcFinalize(ierr)
136: end
138: ! -------------------------------------------------------------------
139: !
140: ! MatIsing_Mult - user provided matrix-vector multiply
141: !
142: ! Input Parameters:
143: ! A - matrix
144: ! x - input vector
145: !
146: ! Output Parameter:
147: ! y - output vector
148: !
149: subroutine MatIsing_Mult(A,x,y,ierr)
150: #include <petsc/finclude/petscmat.h>
151: use petscmat
152: implicit none
154: Mat A
155: Vec x,y
156: PetscInt trans,one,N
157: PetscScalar x_array(1),y_array(1)
158: PetscOffset i_x,i_y
159: PetscErrorCode ierr
161: ! The actual routine for the matrix-vector product
162: external mvmisg
164: call MatGetSize(A,N,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
165: call VecGetArrayRead(x,x_array,i_x,ierr);CHKERRQ(ierr)
166: call VecGetArray(y,y_array,i_y,ierr);CHKERRQ(ierr)
168: trans = 0
169: one = 1
170: call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)
172: call VecRestoreArrayRead(x,x_array,i_x,ierr);CHKERRQ(ierr)
173: call VecRestoreArray(y,y_array,i_y,ierr);CHKERRQ(ierr)
175: return
176: end
178: ! -------------------------------------------------------------------
179: ! The actual routine for the matrix-vector product
180: ! See http://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
182: SUBROUTINE MVMISG( TRANS, N, M, X, LDX, Y, LDY )
183: #include <petsc/finclude/petscsys.h>
184: use petscsys
185: ! ..
186: ! .. Scalar Arguments ..
187: PetscInt LDY, LDX, M, N, TRANS
188: ! ..
189: ! .. Array Arguments ..
190: PetscScalar Y( LDY, * ), X( LDX, * )
191: ! ..
192: !
193: ! Purpose
194: ! =======
195: !
196: ! Compute
197: !
198: ! Y(:,1:M) = op(A)*X(:,1:M)
199: !
200: ! where op(A) is A or A' (the transpose of A). The A is the Ising
201: ! matrix.
202: !
203: ! Arguments
204: ! =========
205: !
206: ! TRANS (input) INTEGER
207: ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
208: ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
209: !
210: ! N (input) INTEGER
211: ! The order of the matrix A. N has to be an even number.
212: !
213: ! M (input) INTEGER
214: ! The number of columns of X to multiply.
215: !
216: ! X (input) DOUBLE PRECISION array, dimension ( LDX, M )
217: ! X contains the matrix (vectors) X.
218: !
219: ! LDX (input) INTEGER
220: ! The leading dimension of array X, LDX >= max( 1, N )
221: !
222: ! Y (output) DOUBLE PRECISION array, dimension (LDX, M )
223: ! contains the product of the matrix op(A) with X.
224: !
225: ! LDY (input) INTEGER
226: ! The leading dimension of array Y, LDY >= max( 1, N )
227: !
228: ! ===================================================================
229: !
230: !
232: ! .. Local Variables ..
233: PetscInt I, K
234: PetscReal ALPHA, BETA
235: PetscReal COSA, COSB, SINA, SINB
236: PetscScalar TEMP, TEMP1
237: !
238: ! .. Intrinsic functions ..
239: INTRINSIC COS, SIN
240: !
241: ALPHA = PETSC_PI/4
242: BETA = PETSC_PI/4
243: COSA = COS( ALPHA )
244: SINA = SIN( ALPHA )
245: COSB = COS( BETA )
246: SINB = SIN( BETA )
247: !
248: IF ( TRANS.EQ.0 ) THEN
249: !
250: ! Compute Y(:,1:M) = A*X(:,1:M)
252: DO 30 K = 1, M
253: !
254: Y( 1, K ) = COSB*X( 1, K ) - SINB*X( N, K )
255: DO 10 I = 2, N-1, 2
256: Y( I, K ) = COSB*X( I, K ) + SINB*X( I+1, K )
257: Y( I+1, K ) = -SINB*X( I, K ) + COSB*X( I+1, K )
258: 10 CONTINUE
259: Y( N, K ) = SINB*X( 1, K ) + COSB*X( N, K )
260: !
261: DO 20 I = 1, N, 2
262: TEMP = COSA*Y( I, K ) + SINA*Y( I+1, K )
263: Y( I+1, K ) = -SINA*Y( I, K ) + COSA*Y( I+1, K )
264: Y( I, K ) = TEMP
265: 20 CONTINUE
266: !
267: 30 CONTINUE
268: !
269: ELSE IF ( TRANS.EQ.1 ) THEN
270: !
271: ! Compute Y(:1:M) = A'*X(:,1:M)
272: !
273: DO 60 K = 1, M
274: !
275: DO 40 I = 1, N, 2
276: Y( I, K ) = COSA*X( I, K ) - SINA*X( I+1, K )
277: Y( I+1, K ) = SINA*X( I, K ) + COSA*X( I+1, K )
278: 40 CONTINUE
279: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
280: DO 50 I = 2, N-1, 2
281: TEMP1 = COSB*Y( I, K ) - SINB*Y( I+1, K )
282: Y( I+1, K ) = SINB*Y( I, K ) + COSB*Y( I+1, K )
283: Y( I, K ) = TEMP1
284: 50 CONTINUE
285: Y( N, K ) = -SINB*Y( 1, K ) + COSB*Y( N, K )
286: Y( 1, K ) = TEMP
287: !
288: 60 CONTINUE
289: !
290: END IF
291: !
292: RETURN
293: !
294: ! END OF MVMISG
295: END
297: !/*TEST
298: !
299: ! test:
300: ! suffix: 1
301: ! args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
302: ! requires: !complex
303: !
304: !TEST*/