Actual source code: ex6f90.F90

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

175:       return
176:       end

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

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

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

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

297: !/*TEST
298: !
299: !   test:
300: !      suffix: 1
301: !      args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
302: !      requires: !complex
303: !
304: !TEST*/