Actual source code: ex54f.F90
slepc-3.21.1 2024-04-26
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_DEFAULT_INTEGER,PETSC_DEFAULT_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*/