Actual source code: ex6f.F90

slepc-3.22.1 2024-10-28
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-, 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: Eigensystem from the Ising model for ferromagnetic materials.
 13: !  Information about the model can be found at the following
 14: !  site https://math.nist.gov/MatrixMarket/data/NEP
 15: !
 16: !  The command line options are:
 17: !    -m <m>, where <m> is the number of 2x2 blocks, i.e. matrix size N=2*m
 18: !
 19: ! ----------------------------------------------------------------------
 20: !
 21:       program main
 22: #include <slepc/finclude/slepceps.h>
 23:       use slepceps
 24:       implicit none

 26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 27: !     Declarations
 28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 29: !
 30: !  Variables:
 31: !     A     operator matrix
 32: !     eps   eigenproblem solver context

 34:       Mat            A
 35:       EPS            eps
 36:       EPSType        tname
 37:       PetscReal      tol
 38:       PetscInt       N, m
 39:       PetscInt       nev, maxit, its
 40:       PetscMPIInt    sz, rank
 41:       PetscErrorCode ierr
 42:       PetscBool      flg, terse

 44: !     This is the routine to use for matrix-free approach
 45: !
 46:       external MatIsing_Mult

 48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49: !     Beginning of program
 50: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 52:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
 53: #if defined(PETSC_USE_COMPLEX)
 54:       SETERRA(PETSC_COMM_SELF,PETSC_ERR_SUP,'This example requires real numbers')
 55: #endif
 56:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,sz,ierr))
 57:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 58:       if (sz .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only!'); endif
 59:       m = 30
 60:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr))
 61:       N = 2*m

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

 67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68: !     Register the matrix-vector subroutine for the operator that defines
 69: !     the eigensystem, Ax=kx
 70: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 72:       PetscCallA(MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,PETSC_NULL_INTEGER,A,ierr))
 73:       PetscCallA(MatShellSetOperation(A,MATOP_MULT,MatIsing_Mult,ierr))

 75: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76: !     Create the eigensolver and display info
 77: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 79: !     ** Create eigensolver context
 80:       PetscCallA(EPSCreate(PETSC_COMM_WORLD,eps,ierr))

 82: !     ** Set operators. In this case, it is a standard eigenvalue problem
 83:       PetscCallA(EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr))
 84:       PetscCallA(EPSSetProblemType(eps,EPS_NHEP,ierr))

 86: !     ** Set solver parameters at runtime
 87:       PetscCallA(EPSSetFromOptions(eps,ierr))

 89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90: !     Solve the eigensystem
 91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 93:       PetscCallA(EPSSolve(eps,ierr))
 94:       PetscCallA(EPSGetIterationNumber(eps,its,ierr))
 95:       if (rank .eq. 0) then
 96:         write(*,'(A,I4)') ' Number of iterations of the method: ',its
 97:       endif

 99: !     ** Optional: Get some information from the solver and display it
100:       PetscCallA(EPSGetType(eps,tname,ierr))
101:       if (rank .eq. 0) then
102:         write(*,'(A,A)') ' Solution method: ', tname
103:       endif
104:       PetscCallA(EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
105:       if (rank .eq. 0) then
106:         write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
107:       endif
108:       PetscCallA(EPSGetTolerances(eps,tol,maxit,ierr))
109:       if (rank .eq. 0) then
110:         write(*,'(A,1PE11.4,A,I6)') ' Stopping condition: tol=',tol,', maxit=', maxit
111:       endif

113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: !     Display solution and clean up
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

117: !     ** show detailed info unless -terse option is given by user
118:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
119:       if (terse) then
120:         PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
121:       else
122:         PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
123:         PetscCallA(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr))
124:         PetscCallA(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
125:         PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
126:       endif
127:       PetscCallA(EPSDestroy(eps,ierr))
128:       PetscCallA(MatDestroy(A,ierr))
129:       PetscCallA(SlepcFinalize(ierr))
130:       end

132: ! -------------------------------------------------------------------
133: !
134: !   MatIsing_Mult - user provided matrix-vector multiply
135: !
136: !   Input Parameters:
137: !   A - matrix
138: !   x - input vector
139: !
140: !   Output Parameter:
141: !   y - output vector
142: !
143:       subroutine MatIsing_Mult(A,x,y,ierr)
144: #include <petsc/finclude/petscmat.h>
145:       use petscmat
146:       implicit none

148:       Mat            A
149:       Vec            x,y
150:       PetscInt       trans,one,N
151:       PetscScalar    x_array(1),y_array(1)
152:       PetscOffset    i_x,i_y
153:       PetscErrorCode ierr

155: !     The actual routine for the matrix-vector product
156:       external mvmisg

158:       PetscCall(MatGetSize(A,N,PETSC_NULL_INTEGER,ierr))
159:       PetscCall(VecGetArrayRead(x,x_array,i_x,ierr))
160:       PetscCall(VecGetArray(y,y_array,i_y,ierr))

162:       trans = 0
163:       one = 1
164:       call mvmisg(trans,N,one,x_array(i_x+1),N,y_array(i_y+1),N)

166:       PetscCall(VecRestoreArrayRead(x,x_array,i_x,ierr))
167:       PetscCall(VecRestoreArray(y,y_array,i_y,ierr))

169:       end

171: ! -------------------------------------------------------------------
172: !     The actual routine for the matrix-vector product
173: !     See https://math.nist.gov/MatrixMarket/data/NEP/mvmisg/mvmisg.html

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

225: !     .. Local Variables ..
226:       PetscInt    I, K
227:       PetscReal   ALPHA, BETA
228:       PetscReal   COSA, COSB, SINA, SINB
229:       PetscScalar TEMP, TEMP1
230: !
231: !     .. Intrinsic functions ..
232:       INTRINSIC   COS, SIN
233: !
234:       ALPHA = PETSC_PI/4
235:       BETA = PETSC_PI/4
236:       COSA = COS(ALPHA)
237:       SINA = SIN(ALPHA)
238:       COSB = COS(BETA)
239:       SINB = SIN(BETA)
240: !
241:       IF (TRANS.EQ.0) THEN
242: !
243: !     Compute Y(:,1:M) = A*X(:,1:M)

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

290: !/*TEST
291: !
292: !   test:
293: !      suffix: 1
294: !      args: -eps_max_it 1000 -eps_ncv 12 -eps_tol 1e-5 -eps_nev 4 -eps_largest_imaginary -terse
295: !      requires: !complex
296: !
297: !TEST*/