Actual source code: ex6f.F

slepc-3.6.3 2016-03-29
Report Typos and Errors
  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