Actual source code: ex54f.F90

slepc-main 2024-11-22
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> ./ex54f [-help] [-n <n>] [all SLEPc options]
 11: !
 12: !  Description: Illustrates use of shell matrices in callback interface from Fortran.
 13: !  Similar to ex21.c. This one solves a simple diagonal linear eigenproblem as a NEP.
 14: !
 15: !  The command line options are:
 16: !    -n <n>, where <n> = matrix dimension

 18: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 19: !    Modules needed to pass and get the context to/from the Mat
 20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 22:    module shell_ctx
 23: #include <petsc/finclude/petscmat.h>
 24:       use petscmat
 25:       implicit none
 26:       type :: MatCtx
 27:          PetscScalar :: lambda
 28:       end type MatCtx
 29:    end module shell_ctx

 31:    module shell_ctx_interfaces
 32:       use shell_ctx
 33:       implicit none

 35:       interface MatCreateShell
 36:         subroutine MatCreateShell(comm,mloc,nloc,m,n,ctx,mat,ierr)
 37:           use shell_ctx
 38:           MPI_Comm       :: comm
 39:           PetscInt       :: mloc,nloc,m,n
 40:           type(MatCtx)   :: ctx
 41:           Mat            :: mat
 42:           PetscErrorCode :: ierr
 43:         end subroutine MatCreateShell
 44:       end interface MatCreateShell

 46:       interface MatShellSetContext
 47:         subroutine MatShellSetContext(mat,ctx,ierr)
 48:           use shell_ctx
 49:           Mat            :: mat
 50:           type(MatCtx)   :: ctx
 51:           PetscErrorCode :: ierr
 52:         end subroutine MatShellSetContext
 53:       end interface MatShellSetContext

 55:       interface MatShellGetContext
 56:         subroutine MatShellGetContext(mat,ctx,ierr)
 57:           use shell_ctx
 58:           Mat                   :: mat
 59:           type(MatCtx), pointer :: ctx
 60:           PetscErrorCode        :: ierr
 61:         end subroutine matShellGetContext
 62:       end interface MatShellGetContext

 64:    end module shell_ctx_interfaces

 66: !=================================================================================================

 68:    program main
 69: #include <slepc/finclude/slepcnep.h>
 70:       use slepcnep
 71:       use shell_ctx_interfaces
 72:       implicit none

 74: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75: !     Declarations
 76: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 78:       NEP            :: nep
 79:       Mat            :: A,B
 80:       PetscInt       :: n=400,nev=3,nconv
 81:       PetscErrorCode :: ierr
 82:       PetscScalar    :: sigma
 83:       PetscBool      :: flg,terse
 84:       PetscMPIInt    :: rank
 85:       type(MatCtx)   :: ctxA,ctxB

 87:       external MyNEPFunction,MyNEPJacobian,MatMult_A,MatDuplicate_A,MatMult_B

 89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90: !     Beginning of program
 91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 93:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
 94:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
 95:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 96:       if (rank .eq. 0) then
 97:          write(*,'(/A,I4)') 'Nonlinear eigenproblem with shell matrices, n =',n
 98:       endif

100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: !     Create matrix-free operators for A and B corresponding to T and T'
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

104:       ctxA%lambda = 0.0
105:       PetscCallA(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,ctxA,A,ierr))
106:       PetscCallA(MatShellSetOperation(A,MATOP_MULT,MatMult_A,ierr))
107:       PetscCallA(MatShellSetOperation(A,MATOP_DUPLICATE,MatDuplicate_A,ierr))

109:       ctxB%lambda = 0.0   ! unused
110:       PetscCallA(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,ctxB,B,ierr))
111:       PetscCallA(MatShellSetOperation(B,MATOP_MULT,MatMult_B,ierr))

113: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: !     Create the eigensolver and set various options
115: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

117:       PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
118:       PetscCallA(NEPSetFunction(nep,A,A,MyNEPFunction,PETSC_NULL_INTEGER,ierr))
119:       PetscCallA(NEPSetJacobian(nep,B,MyNEPJacobian,PETSC_NULL_INTEGER,ierr))
120:       PetscCallA(NEPSetDimensions(nep,nev,PETSC_DETERMINE_INTEGER,PETSC_DETERMINE_INTEGER,ierr))
121:       sigma = 1.05
122:       PetscCallA(NEPSetTarget(nep,sigma,ierr))
123:       PetscCallA(NEPSetFromOptions(nep,ierr))

125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: !     Solve the eigensystem, display solution and clean up
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

129:       PetscCallA(NEPSolve(nep, ierr))

131:       PetscCallA(NEPGetConverged(nep,nconv,ierr))
132:       if (rank .eq. 0) then
133:          write(*,'(A,I2/)') ' Number of converged approximate eigenpairs:',nconv
134:       endif
135: !
136: !     ** show detailed info unless -terse option is given by user
137:       PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
138:       if (terse) then
139:          PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr))
140:       else
141:          PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
142:          PetscCallA(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr))
143:          PetscCallA(NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr))
144:          PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
145:       endif

147:       PetscCallA(NEPDestroy(nep,ierr))
148:       PetscCallA(MatDestroy(A,ierr))
149:       PetscCallA(MatDestroy(B,ierr))
150:       PetscCallA(SlepcFinalize(ierr))

152:    end program main

154: ! --------------------------------------------------------------
155: !
156: !  MyNEPFunction - Computes Function matrix  T(lambda)
157: !
158:    subroutine MyNEPFunction(nep,lambda,T,P,ctx,ierr)
159:       use slepcnep
160:       use shell_ctx_interfaces
161:       implicit none

163:       NEP                   :: nep
164:       PetscScalar           :: lambda
165:       Mat                   :: T,P
166:       PetscInt              :: ctx
167:       PetscErrorCode        :: ierr
168:       type(MatCtx), pointer :: ctxT

170:       PetscCall(MatShellGetContext(T,ctxT,ierr))
171:       ctxT%lambda = lambda

173:    end subroutine MyNEPFunction

175: ! --------------------------------------------------------------
176: !
177: !  MyNEPJacobian - Computes Jacobian matrix  T'(lambda)
178: !
179:    subroutine MyNEPJacobian(nep,lambda,T,ctx,ierr)
180:       use slepcnep
181:       use shell_ctx_interfaces
182:       implicit none

184:       NEP                   :: nep
185:       PetscScalar           :: lambda
186:       Mat                   :: T
187:       PetscInt              :: ctx
188:       PetscErrorCode        :: ierr
189:       type(MatCtx), pointer :: ctxT

191:       PetscCall(MatShellGetContext(T,ctxT,ierr))
192:       ctxT%lambda = lambda

194:    end subroutine MyNEPJacobian

196: ! --------------------------------------------------------------
197: !
198: !  MatMult_A - Shell matrix operation, multiples y=A*x
199: !  Here A=(D-lambda*I) where D is a diagonal matrix
200: !
201:    subroutine MatMult_A(A,x,y,ierr)
202:       use shell_ctx_interfaces
203:       implicit none

205:       Mat                   :: A
206:       Vec                   :: x,y
207:       PetscErrorCode        :: ierr
208:       PetscInt              :: i,istart,iend
209:       PetscScalar           :: val
210:       type(MatCtx), pointer :: ctxA
211: !
212:       PetscCall(VecGetOwnershipRange(x,istart,iend,ierr))
213:       do i=istart,iend-1
214:          val = i+1
215:          val = 1.0/val
216:          PetscCall(VecSetValue(y,i,val,INSERT_VALUES,ierr))
217:       end do
218:       PetscCall(VecAssemblyBegin(y,ierr))
219:       PetscCall(VecAssemblyEnd(y,ierr))
220:       PetscCall(VecPointwiseMult(y,y,x,ierr))
221:       PetscCall(MatShellGetContext(A,ctxA,ierr))
222:       PetscCall(VecAXPY(y,-ctxA%lambda,x,ierr))

224:    end subroutine MatMult_A

226: ! --------------------------------------------------------------
227: !
228: !  MatDuplicate_A - Shell matrix operation, duplicates A
229: !
230:    subroutine MatDuplicate_A(A,opt,M,ierr)
231:       use shell_ctx_interfaces
232:       implicit none

234:       Mat                   :: A,M
235:       MatDuplicateOption    :: opt
236:       PetscErrorCode        :: ierr
237:       PetscInt              :: ml,nl
238:       type(MatCtx), pointer :: ctxM,ctxA

240:       external MatMult_A,MatDestroy_A

242:       PetscCall(MatGetLocalSize(A,ml,nl,ierr));
243:       PetscCall(MatShellGetContext(A,ctxA,ierr))
244:       allocate(ctxM)
245:       ctxM%lambda = ctxA%lambda
246:       PetscCall(MatCreateShell(PETSC_COMM_WORLD,ml,nl,PETSC_DETERMINE,PETSC_DETERMINE,ctxM,M,ierr))
247:       PetscCall(MatShellSetOperation(M,MATOP_MULT,MatMult_A,ierr))
248:       PetscCall(MatShellSetOperation(M,MATOP_DESTROY,MatDestroy_A,ierr))

250:    end subroutine MatDuplicate_A

252: ! --------------------------------------------------------------
253: !
254: !  MatDestroy_A - Shell matrix operation, destroys A
255: !
256:    subroutine MatDestroy_A(A,ierr)
257:       use shell_ctx_interfaces
258:       implicit none

260:       Mat                   :: A
261:       PetscErrorCode        :: ierr
262:       type(MatCtx), pointer :: ctxA

264:       PetscCall(MatShellGetContext(A,ctxA,ierr))
265:       deallocate(ctxA)

267:    end subroutine MatDestroy_A

269: ! --------------------------------------------------------------
270: !
271: !  MatMult_B - Shell matrix operation, multiples y=B*x
272: !  Here B=-I
273: !
274:    subroutine MatMult_B(B,x,y,ierr)
275:       use petscmat
276:       implicit none

278:       Mat            :: B
279:       Vec            :: x
280:       Vec            :: y
281:       PetscErrorCode :: ierr
282:       PetscScalar    :: mone

284:       PetscCall(VecCopy(x,y,ierr))
285:       mone = -1.0
286:       PetscCall(VecScale(y,mone,ierr))

288:    end subroutine MatMult_B

290: !/*TEST
291: !
292: !   testset:
293: !      args: -terse
294: !      output_file: output/ex54f_1.out
295: !      filter: grep -v approximate | sed -e "s/[+-]0\.0*i//g"
296: !      test:
297: !         suffix: 1_slp
298: !         args: -nep_type slp -nep_slp_ksp_type gmres -nep_slp_pc_type none
299: !         requires: double
300: !      test:
301: !         suffix: 1_nleigs
302: !         args: -nep_type nleigs -rg_interval_endpoints 0.2,1.1 -nep_nleigs_ksp_type gmres -nep_nleigs_pc_type none
303: !         requires: !complex
304: !      test:
305: !         suffix: 1_nleigs_complex
306: !         args: -nep_type nleigs -rg_interval_endpoints 0.2,1.1,-.1,.1 -nep_nleigs_ksp_type gmres -nep_nleigs_pc_type none
307: !         requires: complex
308: !
309: !TEST*/