Actual source code: ex6f.F

slepc-3.7.2 2016-07-19
Report Typos and Errors
  1: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  2: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  3: !  Copyright (c) 2002-2016, 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: mpiexec -n <np> ./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_OBJECT,PETSC_NULL_CHARACTER,   &
 83:      &                        '-m',m,flg,ierr)
 84:       N = 2*m

 86:       if (rank .eq. 0) then
 87:         write(*,*)
 88:         write(*,'(A,I6,A)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
 89:         write(*,*)
 90:       endif

 92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 93: !     Register the matrix-vector subroutine for the operator that defines
 94: !     the eigensystem, Ax=kx
 95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

 97:       call MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_OBJECT,   &
 98:      &                    A,ierr)
 99:       call MatShellSetOperation(A,MATOP_MULT,MatIsing_Mult,ierr)

101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
102: !     Create the eigensolver and display info
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

105: !     ** Create eigensolver context
106:       call EPSCreate(PETSC_COMM_WORLD,eps,ierr)

108: !     ** Set operators. In this case, it is a standard eigenvalue problem
109:       call EPSSetOperators(eps,A,PETSC_NULL_OBJECT,ierr)
110:       call EPSSetProblemType(eps,EPS_NHEP,ierr)

112: !     ** Set solver parameters at runtime
113:       call EPSSetFromOptions(eps,ierr)

115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
116: !     Solve the eigensystem
117: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

119:       call EPSSolve(eps,ierr) 
120:       call EPSGetIterationNumber(eps,its,ierr)
121:       if (rank .eq. 0) then
122:         write(*,'(A,I4)') ' Number of iterations of the method: ', its
123:       endif

125: !     ** Optional: Get some information from the solver and display it
126:       call EPSGetType(eps,tname,ierr)
127:       if (rank .eq. 0) then
128:         write(*,'(A,A)') ' Solution method: ', tname
129:       endif
130:       call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,                 &
131:      &                      PETSC_NULL_INTEGER,ierr)
132:       if (rank .eq. 0) then
133:         write(*,'(A,I2)') ' Number of requested eigenvalues:', nev
134:       endif
135:       call EPSGetTolerances(eps,tol,maxit,ierr)
136:       if (rank .eq. 0) then
137:         write(*,'(A,1PE10.4,A,I6)') ' Stopping condition: tol=', tol,   &
138:      &                              ', maxit=', maxit
139:       endif

141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
142: !     Display solution and clean up
143: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

145: !     ** show detailed info unless -terse option is given by user
146:       call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER,  &
147:      &                        '-terse',terse,ierr)
148:       if (terse) then
149:         call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_OBJECT,ierr)
150:       else
151:         call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,           &
152:      &                   PETSC_VIEWER_ASCII_INFO_DETAIL,ierr)
153:         call EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr)
154:         call EPSErrorView(eps,EPS_ERROR_RELATIVE,                       &
155:      &                   PETSC_VIEWER_STDOUT_WORLD,ierr)
156:         call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr)
157:       endif
158:       call EPSDestroy(eps,ierr)
159:       call MatDestroy(A,ierr)

161: #if defined(PETSC_USE_COMPLEX)
162:  999  continue
163: #endif
164:       call SlepcFinalize(ierr)
165:       end

167: ! ------------------------------------------------------------------- 
168: ! 
169: !   MatIsing_Mult - user provided matrix-vector multiply 
170: !
171: !   Input Parameters:
172: !   A - matrix
173: !   x - input vector
174: !
175: !   Output Parameter:
176: !   y - output vector
177: ! 
178:       subroutine MatIsing_Mult(A,x,y,ierr)
179:       implicit none

181: #include <petsc/finclude/petscsys.h>
182: #include <petsc/finclude/petscvec.h>
183: #include <petsc/finclude/petscmat.h>

185:       Mat            A
186:       Vec            x,y
187:       PetscInt       trans,one,N
188:       PetscScalar    x_array(1),y_array(1)
189:       PetscOffset    i_x,i_y
190:       PetscErrorCode ierr

192: !     The actual routine for the matrix-vector product
193:       external mvmisg

195:       call MatGetSize(A,N,PETSC_NULL_INTEGER,ierr)
196:       call VecGetArrayRead(x,x_array,i_x,ierr)
197:       call VecGetArray(y,y_array,i_y,ierr)

199:       trans = 0
200:       one = 1
201:       call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)

203:       call VecRestoreArrayRead(x,x_array,i_x,ierr)
204:       call VecRestoreArray(y,y_array,i_y,ierr)

206:       return
207:       end

209: ! ------------------------------------------------------------------- 
210: !     The actual routine for the matrix-vector product
211: !     See http://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html

213:       SUBROUTINE MVMISG( TRANS, N, M, X, LDX, Y, LDY )
214: !     ..
215: !     .. Scalar Arguments ..
216:       PetscInt     LDY, LDX, M, N, TRANS
217: !     ..
218: !     .. Array Arguments ..
219:       PetscScalar  Y( LDY, * ), X( LDX, * )
220: !     ..
221: !
222: !  Purpose
223: !  =======
224: !
225: !  Compute
226: !
227: !               Y(:,1:M) = op(A)*X(:,1:M)
228: !
229: !  where op(A) is A or A' (the transpose of A). The A is the Ising 
230: !  matrix.
231: !
232: !  Arguments
233: !  =========
234: !
235: !  TRANS   (input) INTEGER
236: !          If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M) 
237: !          If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M) 
238: !           
239: !  N       (input) INTEGER
240: !          The order of the matrix A. N has to be an even number.
241: !
242: !  M       (input) INTEGER
243: !          The number of columns of X to multiply.
244: !
245: !  X       (input) DOUBLE PRECISION array, dimension ( LDX, M )
246: !          X contains the matrix (vectors) X.
247: !
248: !  LDX     (input) INTEGER
249: !          The leading dimension of array X, LDX >= max( 1, N )
250: !
251: !  Y       (output) DOUBLE PRECISION array, dimension (LDX, M )
252: !          contains the product of the matrix op(A) with X.
253: !
254: !  LDY     (input) INTEGER
255: !          The leading dimension of array Y, LDY >= max( 1, N )
256: !
257: !  ===================================================================
258: !
259: !
260: #include <petsc/finclude/petscsys.h>

262: !     .. Local Variables .. 
263:       PetscInt    I, K 
264:       PetscReal   ALPHA, BETA
265:       PetscReal   COSA, COSB, SINA
266:       PetscReal   SINB, TEMP, TEMP1 
267: !
268: !     .. Intrinsic functions ..
269:       INTRINSIC   COS, SIN 
270: !
271:       ALPHA = PETSC_PI/4
272:       BETA = PETSC_PI/4
273:       COSA = COS( ALPHA ) 
274:       SINA = SIN( ALPHA ) 
275:       COSB = COS( BETA )
276:       SINB = SIN( BETA ) 
277: !      
278:       IF ( TRANS.EQ.0 ) THEN 
279: !
280: !     Compute Y(:,1:M) = A*X(:,1:M)

282:          DO 30 K = 1, M
283: !
284:             Y( 1, K ) = COSB*X( 1, K ) - SINB*X( N, K ) 
285:             DO 10 I = 2, N-1, 2   
286:                Y( I, K )   =  COSB*X( I, K ) + SINB*X( I+1, K )
287:                Y( I+1, K ) = -SINB*X( I, K ) + COSB*X( I+1, K )  
288:    10       CONTINUE
289:             Y( N, K ) = SINB*X( 1, K ) + COSB*X( N, K ) 
290: !
291:             DO 20 I = 1, N, 2
292:                TEMP        =  COSA*Y( I, K ) + SINA*Y( I+1, K )
293:                Y( I+1, K ) = -SINA*Y( I, K ) + COSA*Y( I+1, K )  
294:                Y( I, K )   = TEMP 
295:    20       CONTINUE  
296: !
297:    30    CONTINUE 
298: !
299:       ELSE IF ( TRANS.EQ.1 ) THEN 
300: !
301: !        Compute Y(:1:M) = A'*X(:,1:M) 
302: !
303:          DO 60 K = 1, M 
304: !
305:             DO 40 I = 1, N, 2
306:                Y( I, K )   =  COSA*X( I, K ) - SINA*X( I+1, K )
307:                Y( I+1, K ) =  SINA*X( I, K ) + COSA*X( I+1, K )  
308:    40       CONTINUE  
309:             TEMP  = COSB*Y(1,K) + SINB*Y(N,K) 
310:             DO 50 I = 2, N-1, 2   
311:                TEMP1       =  COSB*Y( I, K ) - SINB*Y( I+1, K )
312:                Y( I+1, K ) =  SINB*Y( I, K ) + COSB*Y( I+1, K )  
313:                Y( I, K )   =  TEMP1
314:    50       CONTINUE
315:             Y( N, K ) = -SINB*Y( 1, K ) + COSB*Y( N, K ) 
316:             Y( 1, K ) = TEMP 
317: !
318:    60    CONTINUE
319: !
320:       END IF 
321: !
322:       RETURN
323: !  
324: !     END OF MVMISG
325:       END