Actual source code: ex6f90.F90

slepc-3.8.3 2018-04-03
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2017, 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:       SETERRA(PETSC_COMM_SELF,1,'This example requires real numbers')
 60: #endif
 61:       call MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr);CHKERRA(ierr)
 62:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
 63:       if (sz .ne. 1) then
 64:         SETERRA(PETSC_COMM_SELF,1,'This is a uniprocessor example only!')
 65:       endif
 66:       m = 30
 67:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr);CHKERRA(ierr)
 68:       N = 2*m

 70:       if (rank .eq. 0) then
 71:         write(*,'(/A,I6,A/)') 'Ising Model Eigenproblem, m=',m,', (N=2*m)'
 72:       endif

 74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75: !     Register the matrix-vector subroutine for the operator that defines
 76: !     the eigensystem, Ax=kx
 77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 79:       call MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_INTEGER,A,ierr);CHKERRA(ierr)
 80:       call MatShellSetOperation(A,MATOP_MULT,MatIsing_Mult,ierr);CHKERRA(ierr)

 82: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83: !     Create the eigensolver and display info
 84: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 86: !     ** Create eigensolver context
 87:       call EPSCreate(PETSC_COMM_WORLD,eps,ierr);CHKERRA(ierr)

 89: !     ** Set operators. In this case, it is a standard eigenvalue problem
 90:       call EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr);CHKERRA(ierr)
 91:       call EPSSetProblemType(eps,EPS_NHEP,ierr);CHKERRA(ierr)

 93: !     ** Set solver parameters at runtime
 94:       call EPSSetFromOptions(eps,ierr);CHKERRA(ierr)

 96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97: !     Solve the eigensystem
 98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

100:       call EPSSolve(eps,ierr);CHKERRA(ierr)
101:       call EPSGetIterationNumber(eps,its,ierr);CHKERRA(ierr)
102:       if (rank .eq. 0) then
103:         write(*,'(A,I4)') ' Number of iterations of the method: ',its
104:       endif

106: !     ** Optional: Get some information from the solver and display it
107:       call EPSGetType(eps,tname,ierr);CHKERRA(ierr)
108:       if (rank .eq. 0) then
109:         write(*,'(A,A)') ' Solution method: ', tname
110:       endif
111:       call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
112:       if (rank .eq. 0) then
113:         write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
114:       endif
115:       call EPSGetTolerances(eps,tol,maxit,ierr);CHKERRA(ierr)
116:       if (rank .eq. 0) then
117:         write(*,'(A,1PE10.4,A,I6)') ' Stopping condition: tol=',tol,', maxit=', maxit
118:       endif

120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: !     Display solution and clean up
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

124: !     ** show detailed info unless -terse option is given by user
125:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
126:       if (terse) then
127:         call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
128:       else
129:         call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
130:         call EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
131:         call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
132:         call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
133:       endif
134:       call EPSDestroy(eps,ierr);CHKERRA(ierr)
135:       call MatDestroy(A,ierr);CHKERRA(ierr)
136:       call SlepcFinalize(ierr)
137:       end

139: ! -------------------------------------------------------------------
140: !
141: !   MatIsing_Mult - user provided matrix-vector multiply
142: !
143: !   Input Parameters:
144: !   A - matrix
145: !   x - input vector
146: !
147: !   Output Parameter:
148: !   y - output vector
149: !
150:       subroutine MatIsing_Mult(A,x,y,ierr)
151: #include <petsc/finclude/petscmat.h>
152:       use petscmat
153:       implicit none

155:       Mat            A
156:       Vec            x,y
157:       PetscInt       trans,one,N
158:       PetscScalar    x_array(1),y_array(1)
159:       PetscOffset    i_x,i_y
160:       PetscErrorCode ierr

162: !     The actual routine for the matrix-vector product
163:       external mvmisg

165:       call MatGetSize(A,N,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
166:       call VecGetArrayRead(x,x_array,i_x,ierr);CHKERRQ(ierr)
167:       call VecGetArray(y,y_array,i_y,ierr);CHKERRQ(ierr)

169:       trans = 0
170:       one = 1
171:       call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)

173:       call VecRestoreArrayRead(x,x_array,i_x,ierr);CHKERRQ(ierr)
174:       call VecRestoreArray(y,y_array,i_y,ierr);CHKERRQ(ierr)

176:       return
177:       end

179: ! -------------------------------------------------------------------
180: !     The actual routine for the matrix-vector product
181: !     See http://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html

183:       SUBROUTINE MVMISG( TRANS, N, M, X, LDX, Y, LDY )
184: #include <petsc/finclude/petscsys.h>
185:       use petscsys      
186: !     ..
187: !     .. Scalar Arguments ..
188:       PetscInt     LDY, LDX, M, N, TRANS
189: !     ..
190: !     .. Array Arguments ..
191:       PetscScalar  Y( LDY, * ), X( LDX, * )
192: !     ..
193: !
194: !  Purpose
195: !  =======
196: !
197: !  Compute
198: !
199: !               Y(:,1:M) = op(A)*X(:,1:M)
200: !
201: !  where op(A) is A or A' (the transpose of A). The A is the Ising
202: !  matrix.
203: !
204: !  Arguments
205: !  =========
206: !
207: !  TRANS   (input) INTEGER
208: !          If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
209: !          If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
210: !
211: !  N       (input) INTEGER
212: !          The order of the matrix A. N has to be an even number.
213: !
214: !  M       (input) INTEGER
215: !          The number of columns of X to multiply.
216: !
217: !  X       (input) DOUBLE PRECISION array, dimension ( LDX, M )
218: !          X contains the matrix (vectors) X.
219: !
220: !  LDX     (input) INTEGER
221: !          The leading dimension of array X, LDX >= max( 1, N )
222: !
223: !  Y       (output) DOUBLE PRECISION array, dimension (LDX, M )
224: !          contains the product of the matrix op(A) with X.
225: !
226: !  LDY     (input) INTEGER
227: !          The leading dimension of array Y, LDY >= max( 1, N )
228: !
229: !  ===================================================================
230: !
231: !

233: !     .. Local Variables ..
234:       PetscInt    I, K
235:       PetscReal   ALPHA, BETA
236:       PetscReal   COSA, COSB, SINA
237:       PetscReal   SINB, TEMP, TEMP1
238: !
239: !     .. Intrinsic functions ..
240:       INTRINSIC   COS, SIN
241: !
242:       ALPHA = PETSC_PI/4
243:       BETA = PETSC_PI/4
244:       COSA = COS( ALPHA )
245:       SINA = SIN( ALPHA )
246:       COSB = COS( BETA )
247:       SINB = SIN( BETA )
248: !
249:       IF ( TRANS.EQ.0 ) THEN
250: !
251: !     Compute Y(:,1:M) = A*X(:,1:M)

253:          DO 30 K = 1, M
254: !
255:             Y( 1, K ) = COSB*X( 1, K ) - SINB*X( N, K )
256:             DO 10 I = 2, N-1, 2
257:                Y( I, K )   =  COSB*X( I, K ) + SINB*X( I+1, K )
258:                Y( I+1, K ) = -SINB*X( I, K ) + COSB*X( I+1, K )
259:    10       CONTINUE
260:             Y( N, K ) = SINB*X( 1, K ) + COSB*X( N, K )
261: !
262:             DO 20 I = 1, N, 2
263:                TEMP        =  COSA*Y( I, K ) + SINA*Y( I+1, K )
264:                Y( I+1, K ) = -SINA*Y( I, K ) + COSA*Y( I+1, K )
265:                Y( I, K )   = TEMP
266:    20       CONTINUE
267: !
268:    30    CONTINUE
269: !
270:       ELSE IF ( TRANS.EQ.1 ) THEN
271: !
272: !        Compute Y(:1:M) = A'*X(:,1:M)
273: !
274:          DO 60 K = 1, M
275: !
276:             DO 40 I = 1, N, 2
277:                Y( I, K )   =  COSA*X( I, K ) - SINA*X( I+1, K )
278:                Y( I+1, K ) =  SINA*X( I, K ) + COSA*X( I+1, K )
279:    40       CONTINUE
280:             TEMP  = COSB*Y(1,K) + SINB*Y(N,K)
281:             DO 50 I = 2, N-1, 2
282:                TEMP1       =  COSB*Y( I, K ) - SINB*Y( I+1, K )
283:                Y( I+1, K ) =  SINB*Y( I, K ) + COSB*Y( I+1, K )
284:                Y( I, K )   =  TEMP1
285:    50       CONTINUE
286:             Y( N, K ) = -SINB*Y( 1, K ) + COSB*Y( N, K )
287:             Y( 1, K ) = TEMP
288: !
289:    60    CONTINUE
290: !
291:       END IF
292: !
293:       RETURN
294: !
295: !     END OF MVMISG
296:       END