Actual source code: ex6f90.F90

slepc-3.20.1 2023-11-27
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:       return
170:       end

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

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

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

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

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