Actual source code: ex6f.F
slepc-3.6.1 2015-09-03
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
4: !
5: ! This file is part of SLEPc.
6: !
7: ! SLEPc is free software: you can redistribute it and/or modify it under the
8: ! terms of version 3 of the GNU Lesser General Public License as published by
9: ! the Free Software Foundation.
10: !
11: ! SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
12: ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
13: ! FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
14: ! more details.
15: !
16: ! You should have received a copy of the GNU Lesser General Public License
17: ! along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: !
20: ! Program usage: mpirun -np n ex6f [-help] [-m <m>] [all SLEPc options]
21: !
22: ! Description: This example solves the eigensystem arising in the Ising
23: ! model for ferromagnetic materials. The file mvmisg.f must be linked
24: ! together. Information about the model can be found at the following
25: ! site http://math.nist.gov/MatrixMarket/data/NEP
26: !
27: ! The command line options are:
28: ! -m <m>, where <m> is the number of 2x2 blocks, i.e. matrix size N=2*m
29: !
30: ! ----------------------------------------------------------------------
31: !
32: program main
33: implicit none
35: #include <petsc/finclude/petscsys.h>
36: #include <petsc/finclude/petscvec.h>
37: #include <petsc/finclude/petscmat.h>
38: #include <petsc/finclude/petscviewer.h>
39: #include <slepc/finclude/slepcsys.h>
40: #include <slepc/finclude/slepceps.h>
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Declarations
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: !
46: ! Variables:
47: ! A operator matrix
48: ! eps eigenproblem solver context
50: Mat A
51: EPS eps
52: EPSType tname
53: PetscReal tol
54: PetscInt N, m
55: PetscInt nev, maxit, its
56: PetscMPIInt sz, rank
57: PetscErrorCode ierr
58: PetscBool flg, terse
60: ! This is the routine to use for matrix-free approach
61: !
62: external MatIsing_Mult
64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
65: ! Beginning of program
66: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
69: #if defined(PETSC_USE_COMPLEX)
70: write(*,*) 'This example requires real numbers.'
71: goto 999
72: #endif
73: call MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr)
74: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
75: if (sz .ne. 1) then
76: if (rank .eq. 0) then
77: write(*,*) 'This is a uniprocessor example only!'
78: endif
79: SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
80: endif
81: m = 30
82: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
83: N = 2*m
85: if (rank .eq. 0) then
86: write(*,*)
87: write(*,'(A,I6,A)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
88: write(*,*)
89: endif
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: ! Register the matrix-vector subroutine for the operator that defines
93: ! the eigensystem, Ax=kx
94: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: call MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_OBJECT, &
97: & A,ierr)
98: call MatShellSetOperation(A,MATOP_MULT,MatIsing_Mult,ierr)
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Create the eigensolver and display info
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: ! ** Create eigensolver context
105: call EPSCreate(PETSC_COMM_WORLD,eps,ierr)
107: ! ** Set operators. In this case, it is a standard eigenvalue problem
108: call EPSSetOperators(eps,A,PETSC_NULL_OBJECT,ierr)
109: call EPSSetProblemType(eps,EPS_NHEP,ierr)
111: ! ** Set solver parameters at runtime
112: call EPSSetFromOptions(eps,ierr)
114: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: ! Solve the eigensystem
116: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: call EPSSolve(eps,ierr)
119: call EPSGetIterationNumber(eps,its,ierr)
120: if (rank .eq. 0) then
121: write(*,'(A,I4)') ' Number of iterations of the method: ', its
122: endif
124: ! ** Optional: Get some information from the solver and display it
125: call EPSGetType(eps,tname,ierr)
126: if (rank .eq. 0) then
127: write(*,'(A,A)') ' Solution method: ', tname
128: endif
129: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER, &
130: & PETSC_NULL_INTEGER,ierr)
131: if (rank .eq. 0) then
132: write(*,'(A,I2)') ' Number of requested eigenvalues:', nev
133: endif
134: call EPSGetTolerances(eps,tol,maxit,ierr)
135: if (rank .eq. 0) then
136: write(*,'(A,1PE10.4,A,I6)') ' Stopping condition: tol=', tol, &
137: & ', maxit=', maxit
138: endif
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: ! Display solution and clean up
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: ! ** show detailed info unless -terse option is given by user
145: call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-terse',terse,ierr)
146: if (terse) then
147: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_OBJECT,ierr)
148: else
149: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, &
150: & PETSC_VIEWER_ASCII_INFO_DETAIL,ierr)
151: call EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr)
152: call EPSErrorView(eps,EPS_ERROR_RELATIVE, &
153: & PETSC_VIEWER_STDOUT_WORLD,ierr)
154: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr)
155: endif
156: call EPSDestroy(eps,ierr)
157: call MatDestroy(A,ierr)
159: #if defined(PETSC_USE_COMPLEX)
160: 999 continue
161: #endif
162: call SlepcFinalize(ierr)
163: end
165: ! -------------------------------------------------------------------
166: !
167: ! MatIsing_Mult - user provided matrix-vector multiply
168: !
169: ! Input Parameters:
170: ! A - matrix
171: ! x - input vector
172: !
173: ! Output Parameter:
174: ! y - output vector
175: !
176: subroutine MatIsing_Mult(A,x,y,ierr)
177: implicit none
179: #include <petsc/finclude/petscsys.h>
180: #include <petsc/finclude/petscvec.h>
181: #include <petsc/finclude/petscmat.h>
183: Mat A
184: Vec x,y
185: PetscInt trans,one,N
186: PetscScalar x_array(1),y_array(1)
187: PetscOffset i_x,i_y
188: PetscErrorCode ierr
190: ! The actual routine for the matrix-vector product
191: external mvmisg
193: call MatGetSize(A,N,PETSC_NULL_INTEGER,ierr)
194: call VecGetArrayRead(x,x_array,i_x,ierr)
195: call VecGetArray(y,y_array,i_y,ierr)
197: trans = 0
198: one = 1
199: call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)
201: call VecRestoreArrayRead(x,x_array,i_x,ierr)
202: call VecRestoreArray(y,y_array,i_y,ierr)
204: return
205: end
207: ! -------------------------------------------------------------------
208: ! The actual routine for the matrix-vector product
209: ! See http://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
211: SUBROUTINE MVMISG( TRANS, N, M, X, LDX, Y, LDY )
212: ! ..
213: ! .. Scalar Arguments ..
214: PetscInt LDY, LDX, M, N, TRANS
215: ! ..
216: ! .. Array Arguments ..
217: PetscScalar Y( LDY, * ), X( LDX, * )
218: ! ..
219: !
220: ! Purpose
221: ! =======
222: !
223: ! Compute
224: !
225: ! Y(:,1:M) = op(A)*X(:,1:M)
226: !
227: ! where op(A) is A or A' (the transpose of A). The A is the Ising
228: ! matrix.
229: !
230: ! Arguments
231: ! =========
232: !
233: ! TRANS (input) INTEGER
234: ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
235: ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
236: !
237: ! N (input) INTEGER
238: ! The order of the matrix A. N has to be an even number.
239: !
240: ! M (input) INTEGER
241: ! The number of columns of X to multiply.
242: !
243: ! X (input) DOUBLE PRECISION array, dimension ( LDX, M )
244: ! X contains the matrix (vectors) X.
245: !
246: ! LDX (input) INTEGER
247: ! The leading dimension of array X, LDX >= max( 1, N )
248: !
249: ! Y (output) DOUBLE PRECISION array, dimension (LDX, M )
250: ! contains the product of the matrix op(A) with X.
251: !
252: ! LDY (input) INTEGER
253: ! The leading dimension of array Y, LDY >= max( 1, N )
254: !
255: ! ===================================================================
256: !
257: !
258: #include <petsc/finclude/petscsys.h>
260: ! .. Local Variables ..
261: PetscInt I, K
262: PetscReal ALPHA, BETA
263: PetscReal COSA, COSB, SINA
264: PetscReal SINB, TEMP, TEMP1
265: !
266: ! .. Intrinsic functions ..
267: INTRINSIC COS, SIN
268: !
269: ALPHA = PETSC_PI/4
270: BETA = PETSC_PI/4
271: COSA = COS( ALPHA )
272: SINA = SIN( ALPHA )
273: COSB = COS( BETA )
274: SINB = SIN( BETA )
275: !
276: IF ( TRANS.EQ.0 ) THEN
277: !
278: ! Compute Y(:,1:M) = A*X(:,1:M)
280: DO 30 K = 1, M
281: !
282: Y( 1, K ) = COSB*X( 1, K ) - SINB*X( N, K )
283: DO 10 I = 2, N-1, 2
284: Y( I, K ) = COSB*X( I, K ) + SINB*X( I+1, K )
285: Y( I+1, K ) = -SINB*X( I, K ) + COSB*X( I+1, K )
286: 10 CONTINUE
287: Y( N, K ) = SINB*X( 1, K ) + COSB*X( N, K )
288: !
289: DO 20 I = 1, N, 2
290: TEMP = COSA*Y( I, K ) + SINA*Y( I+1, K )
291: Y( I+1, K ) = -SINA*Y( I, K ) + COSA*Y( I+1, K )
292: Y( I, K ) = TEMP
293: 20 CONTINUE
294: !
295: 30 CONTINUE
296: !
297: ELSE IF ( TRANS.EQ.1 ) THEN
298: !
299: ! Compute Y(:1:M) = A'*X(:,1:M)
300: !
301: DO 60 K = 1, M
302: !
303: DO 40 I = 1, N, 2
304: Y( I, K ) = COSA*X( I, K ) - SINA*X( I+1, K )
305: Y( I+1, K ) = SINA*X( I, K ) + COSA*X( I+1, K )
306: 40 CONTINUE
307: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
308: DO 50 I = 2, N-1, 2
309: TEMP1 = COSB*Y( I, K ) - SINB*Y( I+1, K )
310: Y( I+1, K ) = SINB*Y( I, K ) + COSB*Y( I+1, K )
311: Y( I, K ) = TEMP1
312: 50 CONTINUE
313: Y( N, K ) = -SINB*Y( 1, K ) + COSB*Y( N, K )
314: Y( 1, K ) = TEMP
315: !
316: 60 CONTINUE
317: !
318: END IF
319: !
320: RETURN
321: !
322: ! END OF MVMISG
323: END