Actual source code: ex6f.F
1: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
3: ! Copyright (c) 2002-2011, 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 <finclude/petscsys.h>
36: #include <finclude/petscvec.h>
37: #include <finclude/petscmat.h>
38: #include <finclude/slepcsys.h>
39: #include <finclude/slepceps.h>
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: ! Declarations
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: !
45: ! Variables:
46: ! A operator matrix
47: ! eps eigenproblem solver context
49: Mat A
50: EPS eps
51: EPSType tname
52: PetscReal tol
53: PetscInt N, m
54: PetscInt nev, maxit, its
55: PetscMPIInt sz, rank
56: PetscErrorCode ierr
57: PetscBool flg
59: ! This is the routine to use for matrix-free approach
60: !
61: external MatIsing_Mult
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: ! Beginning of program
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
68: #if defined(PETSC_USE_COMPLEX)
69: write(*,*) 'This example requires real numbers.'
70: goto 999
71: #endif
72: call MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr)
73: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
74: if (sz .ne. 1) then
75: if (rank .eq. 0) then
76: write(*,*) 'This is a uniprocessor example only!'
77: endif
78: SETERRQ(PETSC_COMM_WORLD,1,' ',ierr)
79: endif
80: m = 30
81: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
82: N = 2*m
84: if (rank .eq. 0) then
85: write(*,*)
86: write(*,'(A,I6,A)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
87: write(*,*)
88: endif
90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: ! Register the matrix-vector subroutine for the operator that defines
92: ! the eigensystem, Ax=kx
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: call MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_OBJECT, &
96: & A,ierr)
97: call MatShellSetOperation(A,MATOP_MULT,MatIsing_Mult,ierr)
99: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: ! Create the eigensolver and display info
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: ! ** Create eigensolver context
104: call EPSCreate(PETSC_COMM_WORLD,eps,ierr)
106: ! ** Set operators. In this case, it is a standard eigenvalue problem
107: call EPSSetOperators(eps,A,PETSC_NULL_OBJECT,ierr)
108: call EPSSetProblemType(eps,EPS_NHEP,ierr)
110: ! ** Set solver parameters at runtime
111: call EPSSetFromOptions(eps,ierr)
113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: ! Solve the eigensystem
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: call EPSSolve(eps,ierr)
118: call EPSGetIterationNumber(eps,its,ierr)
119: if (rank .eq. 0) then
120: write(*,'(A,I4)') ' Number of iterations of the method: ', its
121: endif
123: ! ** Optional: Get some information from the solver and display it
124: call EPSGetType(eps,tname,ierr)
125: if (rank .eq. 0) then
126: write(*,'(A,A)') ' Solution method: ', tname
127: endif
128: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER, &
129: & PETSC_NULL_INTEGER,ierr)
130: if (rank .eq. 0) then
131: write(*,'(A,I2)') ' Number of requested eigenvalues:', nev
132: endif
133: call EPSGetTolerances(eps,tol,maxit,ierr)
134: if (rank .eq. 0) then
135: write(*,'(A,1PE10.4,A,I6)') ' Stopping condition: tol=', tol, &
136: & ', maxit=', maxit
137: endif
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: ! Display solution and clean up
141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: call EPSPrintSolution(eps,PETSC_NULL_OBJECT,ierr)
144: call EPSDestroy(eps,ierr)
145: call MatDestroy(A,ierr)
147: #if defined(PETSC_USE_COMPLEX)
148: 999 continue
149: #endif
150: call SlepcFinalize(ierr)
151: end
153: ! -------------------------------------------------------------------
154: !
155: ! MatIsing_Mult - user provided matrix-vector multiply
156: !
157: ! Input Parameters:
158: ! A - matrix
159: ! x - input vector
160: !
161: ! Output Parameter:
162: ! y - output vector
163: !
164: subroutine MatIsing_Mult(A,x,y,ierr)
165: implicit none
167: #include <finclude/petscsys.h>
168: #include <finclude/petscvec.h>
169: #include <finclude/petscmat.h>
171: Mat A
172: Vec x,y
173: PetscInt trans,one,N
174: PetscScalar x_array(1),y_array(1)
175: PetscOffset i_x,i_y
176: PetscErrorCode ierr
178: ! The actual routine for the matrix-vector product
179: external mvmisg
181: call MatGetSize(A,N,PETSC_NULL_INTEGER,ierr)
182: call VecGetArray(x,x_array,i_x,ierr)
183: call VecGetArray(y,y_array,i_y,ierr)
185: trans = 0
186: one = 1
187: call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)
189: call VecRestoreArray(x,x_array,i_x,ierr)
190: call VecRestoreArray(y,y_array,i_y,ierr)
192: return
193: end
195: ! -------------------------------------------------------------------
196: ! The actual routine for the matrix-vector product
197: ! See http://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html
199: SUBROUTINE MVMISG( TRANS, N, M, X, LDX, Y, LDY )
200: ! ..
201: ! .. Scalar Arguments ..
202: PetscInt LDY, LDX, M, N, TRANS
203: ! ..
204: ! .. Array Arguments ..
205: PetscScalar Y( LDY, * ), X( LDX, * )
206: ! ..
207: !
208: ! Purpose
209: ! =======
210: !
211: ! Compute
212: !
213: ! Y(:,1:M) = op(A)*X(:,1:M)
214: !
215: ! where op(A) is A or A' (the transpose of A). The A is the Ising
216: ! matrix.
217: !
218: ! Arguments
219: ! =========
220: !
221: ! TRANS (input) INTEGER
222: ! If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
223: ! If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
224: !
225: ! N (input) INTEGER
226: ! The order of the matrix A. N has to be an even number.
227: !
228: ! M (input) INTEGER
229: ! The number of columns of X to multiply.
230: !
231: ! X (input) DOUBLE PRECISION array, dimension ( LDX, M )
232: ! X contains the matrix (vectors) X.
233: !
234: ! LDX (input) INTEGER
235: ! The leading dimension of array X, LDX >= max( 1, N )
236: !
237: ! Y (output) DOUBLE PRECISION array, dimension (LDX, M )
238: ! contains the product of the matrix op(A) with X.
239: !
240: ! LDY (input) INTEGER
241: ! The leading dimension of array Y, LDY >= max( 1, N )
242: !
243: ! ===================================================================
244: !
245: !
246: ! .. PARAMETERS ..
247: PetscReal PI
248: PARAMETER ( PI = 3.141592653589793D+00 )
249: PetscReal ALPHA, BETA
250: PARAMETER ( ALPHA = PI/4, BETA = PI/4 )
251: !
252: ! .. Local Variables ..
253: PetscInt I, K
254: PetscReal COSA, COSB, SINA
255: PetscReal SINB, TEMP, TEMP1
256: !
257: ! .. Intrinsic functions ..
258: INTRINSIC COS, SIN
259: !
260: COSA = COS( ALPHA )
261: SINA = SIN( ALPHA )
262: COSB = COS( BETA )
263: SINB = SIN( BETA )
264: !
265: IF ( TRANS.EQ.0 ) THEN
266: !
267: ! Compute Y(:,1:M) = A*X(:,1:M)
269: DO 30 K = 1, M
270: !
271: Y( 1, K ) = COSB*X( 1, K ) - SINB*X( N, K )
272: DO 10 I = 2, N-1, 2
273: Y( I, K ) = COSB*X( I, K ) + SINB*X( I+1, K )
274: Y( I+1, K ) = -SINB*X( I, K ) + COSB*X( I+1, K )
275: 10 CONTINUE
276: Y( N, K ) = SINB*X( 1, K ) + COSB*X( N, K )
277: !
278: DO 20 I = 1, N, 2
279: TEMP = COSA*Y( I, K ) + SINA*Y( I+1, K )
280: Y( I+1, K ) = -SINA*Y( I, K ) + COSA*Y( I+1, K )
281: Y( I, K ) = TEMP
282: 20 CONTINUE
283: !
284: 30 CONTINUE
285: !
286: ELSE IF ( TRANS.EQ.1 ) THEN
287: !
288: ! Compute Y(:1:M) = A'*X(:,1:M)
289: !
290: DO 60 K = 1, M
291: !
292: DO 40 I = 1, N, 2
293: Y( I, K ) = COSA*X( I, K ) - SINA*X( I+1, K )
294: Y( I+1, K ) = SINA*X( I, K ) + COSA*X( I+1, K )
295: 40 CONTINUE
296: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
297: DO 50 I = 2, N-1, 2
298: TEMP1 = COSB*Y( I, K ) - SINB*Y( I+1, K )
299: Y( I+1, K ) = SINB*Y( I, K ) + COSB*Y( I+1, K )
300: Y( I, K ) = TEMP1
301: 50 CONTINUE
302: Y( N, K ) = -SINB*Y( 1, K ) + COSB*Y( N, K )
303: Y( 1, K ) = TEMP
304: !
305: 60 CONTINUE
306: !
307: END IF
308: !
309: RETURN
310: !
311: ! END OF MVMISG
312: END