Actual source code: ex54f.F90
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: #include <slepc/finclude/slepcnep.h>
20: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
21: ! Modules needed to pass and get the context to/from the Mat
22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
24: module shell_ctx
25: use petscmat
26: implicit none
27: type :: MatCtx
28: PetscScalar :: lambda
29: end type MatCtx
30: end module shell_ctx
32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
33: ! Module used to implement the shell matrix operations
34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: module ex54fmodule
37: use slepcnep
38: implicit none
40: contains
41: ! --------------------------------------------------------------
42: ! MyNEPFunction - Computes Function matrix T(lambda)
43: !
44: subroutine MyNEPFunction(nep, lambda, T, P, ctx, ierr)
45: use slepcnep
46: use shell_ctx
47: implicit none
49: NEP :: nep
50: PetscScalar :: lambda
51: Mat :: T, P
52: PetscInt :: ctx
53: type(MatCtx), pointer :: ctxT
54: PetscErrorCode, intent(out) :: ierr
56: PetscCall(MatShellGetContext(T, ctxT, ierr))
57: ctxT%lambda = lambda
58: end subroutine MyNEPFunction
60: ! --------------------------------------------------------------
61: ! MyNEPJacobian - Computes Jacobian matrix T'(lambda)
62: !
63: subroutine MyNEPJacobian(nep, lambda, T, ctx, ierr)
64: use slepcnep
65: use shell_ctx
66: implicit none
68: NEP :: nep
69: PetscScalar :: lambda
70: Mat :: T
71: PetscInt :: ctx
72: type(MatCtx), pointer :: ctxT
73: PetscErrorCode, intent(out) :: ierr
75: PetscCall(MatShellGetContext(T, ctxT, ierr))
76: ctxT%lambda = lambda
77: end subroutine MyNEPJacobian
79: ! --------------------------------------------------------------
80: ! MatMult_A - Shell matrix operation, multiples y=A*x
81: ! Here A=(D-lambda*I) where D is a diagonal matrix
82: !
83: subroutine MatMult_A(A, x, y, ierr)
84: use shell_ctx
85: implicit none
87: Mat :: A
88: Vec :: x, y
89: PetscInt :: i, istart, iend
90: PetscScalar :: val
91: type(MatCtx), pointer :: ctxA
92: PetscErrorCode, intent(out) :: ierr
93: !
94: PetscCall(VecGetOwnershipRange(x, istart, iend, ierr))
95: do i = istart, iend - 1
96: val = 1.0/real(i + 1, PETSC_REAL_KIND)
97: PetscCall(VecSetValue(y, i, val, INSERT_VALUES, ierr))
98: end do
99: PetscCall(VecAssemblyBegin(y, ierr))
100: PetscCall(VecAssemblyEnd(y, ierr))
101: PetscCall(VecPointwiseMult(y, y, x, ierr))
102: PetscCall(MatShellGetContext(A, ctxA, ierr))
103: PetscCall(VecAXPY(y, -ctxA%lambda, x, ierr))
104: end subroutine MatMult_A
106: ! --------------------------------------------------------------
107: ! MatDuplicate_A - Shell matrix operation, duplicates A
108: !
109: subroutine MatDuplicate_A(A, opt, M, ierr)
110: use shell_ctx
111: implicit none
113: Mat :: A, M
114: MatDuplicateOption :: opt
115: PetscInt :: ml, nl
116: type(MatCtx), pointer :: ctxM, ctxA
117: PetscErrorCode, intent(out) :: ierr
119: PetscCall(MatGetLocalSize(A, ml, nl, ierr))
120: PetscCall(MatShellGetContext(A, ctxA, ierr))
121: allocate (ctxM)
122: ctxM%lambda = ctxA%lambda
123: PetscCall(MatCreateShell(PETSC_COMM_WORLD, ml, nl, PETSC_DETERMINE, PETSC_DETERMINE, ctxM, M, ierr))
124: PetscCall(MatShellSetOperation(M, MATOP_MULT, MatMult_A, ierr))
125: PetscCall(MatShellSetOperation(M, MATOP_DESTROY, MatDestroy_A, ierr))
126: end subroutine MatDuplicate_A
128: ! --------------------------------------------------------------
129: ! MatDestroy_A - Shell matrix operation, destroys A
130: !
131: subroutine MatDestroy_A(A, ierr)
132: use shell_ctx
133: implicit none
135: Mat :: A
136: type(MatCtx), pointer :: ctxA
137: PetscErrorCode, intent(out) :: ierr
139: PetscCall(MatShellGetContext(A, ctxA, ierr))
140: deallocate (ctxA)
141: end subroutine MatDestroy_A
143: ! --------------------------------------------------------------
144: ! MatMult_B - Shell matrix operation, multiples y=B*x
145: ! Here B=-I
146: !
147: subroutine MatMult_B(B, x, y, ierr)
148: use petscmat
149: implicit none
151: Mat :: B
152: Vec :: x
153: Vec :: y
154: PetscErrorCode, intent(out) :: ierr
155: PetscScalar, parameter :: mone = -1.0
157: PetscCall(VecCopy(x, y, ierr))
158: PetscCall(VecScale(y, mone, ierr))
159: end subroutine MatMult_B
161: end module ex54fmodule
163: !=================================================================================================
165: program ex54f
166: use slepcnep
167: use shell_ctx
168: use ex54fmodule
169: implicit none
171: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: ! Declarations
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: NEP :: nep
176: Mat :: A, B
177: PetscInt :: n = 400, nev = 3, nconv
178: PetscErrorCode :: ierr
179: PetscScalar :: sigma
180: PetscBool :: flg, terse
181: PetscMPIInt :: rank
182: type(MatCtx) :: ctxA, ctxB
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: ! Beginning of program
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
189: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
190: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
191: if (rank == 0) then
192: write (*, '(/a,i4)') 'Nonlinear eigenproblem with shell matrices, n =', n
193: end if
195: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: ! Create matrix-free operators for A and B corresponding to T and T'
197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: ctxA%lambda = 0.0
200: PetscCallA(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, ctxA, A, ierr))
201: PetscCallA(MatShellSetOperation(A, MATOP_MULT, MatMult_A, ierr))
202: PetscCallA(MatShellSetOperation(A, MATOP_DUPLICATE, MatDuplicate_A, ierr))
204: ctxB%lambda = 0.0 ! unused
205: PetscCallA(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, ctxB, B, ierr))
206: PetscCallA(MatShellSetOperation(B, MATOP_MULT, MatMult_B, ierr))
208: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209: ! Create the eigensolver and set various options
210: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212: PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))
213: PetscCallA(NEPSetFunction(nep, A, A, MyNEPFunction, PETSC_NULL_INTEGER, ierr))
214: PetscCallA(NEPSetJacobian(nep, B, MyNEPJacobian, PETSC_NULL_INTEGER, ierr))
215: PetscCallA(NEPSetDimensions(nep, nev, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
216: sigma = 1.05
217: PetscCallA(NEPSetTarget(nep, sigma, ierr))
218: PetscCallA(NEPSetFromOptions(nep, ierr))
220: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
221: ! Solve the eigensystem, display solution and clean up
222: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
224: PetscCallA(NEPSolve(nep, ierr))
226: PetscCallA(NEPGetConverged(nep, nconv, ierr))
227: if (rank == 0) then
228: write (*, '(a,i2/)') ' Number of converged approximate eigenpairs:', nconv
229: end if
230: !
231: ! ** show detailed info unless -terse option is given by user
232: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
233: if (terse) then
234: PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
235: else
236: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
237: PetscCallA(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_WORLD, ierr))
238: PetscCallA(NEPErrorView(nep, NEP_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
239: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
240: end if
242: PetscCallA(NEPDestroy(nep, ierr))
243: PetscCallA(MatDestroy(A, ierr))
244: PetscCallA(MatDestroy(B, ierr))
245: PetscCallA(SlepcFinalize(ierr))
247: end program ex54f
249: !/*TEST
250: !
251: ! testset:
252: ! args: -terse
253: ! output_file: output/ex54f_1.out
254: ! filter: grep -v approximate | sed -e "s/[+-]0\.0*i//g"
255: ! test:
256: ! suffix: 1_slp
257: ! args: -nep_type slp -nep_slp_ksp_type gmres -nep_slp_pc_type none
258: ! requires: double
259: ! test:
260: ! suffix: 1_nleigs
261: ! args: -nep_type nleigs -rg_interval_endpoints 0.2,1.1 -nep_nleigs_ksp_type gmres -nep_nleigs_pc_type none
262: ! requires: !complex
263: ! test:
264: ! suffix: 1_nleigs_complex
265: ! args: -nep_type nleigs -rg_interval_endpoints 0.2,1.1,-.1,.1 -nep_nleigs_ksp_type gmres -nep_nleigs_pc_type none
266: ! requires: complex
267: !
268: !TEST*/