Actual source code: ex10f.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> ./ex10f [-help] [-n <n>] [all SLEPc options]
11: !
12: ! Description: Illustrates the use of shell spectral transformations.
13: ! The problem to be solved is the same as ex1.c and corresponds to the
14: ! Laplacian operator in 1 dimension
15: !
16: ! The command line options are:
17: ! nm <n>, where <n> is the number of grid subdivisions = matrix dimension
18: ! ----------------------------------------------------------------------
20: ! Module contains data needed by shell ST
21: !
22: #include <slepc/finclude/slepceps.h>
23: module ex10fmodule
24: use slepceps
25: implicit none
27: KSP myksp
29: contains
30: ! -------------------------------------------------------------------
31: ! STApply_User - This routine demonstrates the use of a user-provided spectral
32: ! transformation. The transformation implemented in this code is just OP=A^-1.
33: !
34: ! Input Parameters:
35: ! st - spectral transformation context
36: ! x - input vector
37: !
38: ! Output Parameter:
39: ! y - output vector
40: !
41: subroutine STApply_User(st, x, y, ierr)
42: use slepceps
43: implicit none
45: ST :: st
46: Vec :: x, y
47: PetscErrorCode, intent(out) :: ierr
49: PetscCall(KSPSolve(myksp, x, y, ierr))
50: end subroutine
52: ! -------------------------------------------------------------------
53: ! STApplyTranspose_User - This is not required unless using a two-sided eigensolver
54: !
55: ! Input Parameters:
56: ! st - spectral transformation context
57: ! x - input vector
58: !
59: ! Output Parameter:
60: ! y - output vector
61: !
62: subroutine STApplyTranspose_User(st, x, y, ierr)
63: use slepceps
64: implicit none
66: ST :: st
67: Vec :: x, y
68: PetscErrorCode, intent(out) :: ierr
70: PetscCall(KSPSolveTranspose(myksp, x, y, ierr))
71: end subroutine
73: #if defined(PETSC_USE_COMPLEX)
74: ! -------------------------------------------------------------------
75: ! STApplyHermitianTranspose_User - This is not required unless using a two-sided eigensolver
76: ! in complex scalars
77: !
78: ! Input Parameters:
79: ! st - spectral transformation context
80: ! x - input vector
81: !
82: ! Output Parameter:
83: ! y - output vector
84: !
85: subroutine STApplyHermitianTranspose_User(st, x, y, ierr)
86: use slepceps
87: implicit none
89: ST :: st
90: Vec :: x, y, w
91: PetscErrorCode, intent(out) :: ierr
93: PetscCall(VecDuplicate(x, w, ierr))
94: PetscCall(VecCopy(x, w, ierr))
95: PetscCall(VecConjugate(w, ierr))
96: PetscCall(KSPSolveTranspose(myksp, w, y, ierr))
97: PetscCall(VecConjugate(y, ierr))
98: PetscCall(VecDestroy(w, ierr))
99: end subroutine
100: #endif
102: ! -------------------------------------------------------------------
103: ! STBackTransform_User - This routine demonstrates the use of a user-provided spectral
104: ! transformation
105: !
106: ! Input Parameters:
107: ! st - spectral transformation context
108: ! n - number of eigenvalues to transform
109: !
110: ! Output Parameters:
111: ! eigr - real part of eigenvalues
112: ! eigi - imaginary part of eigenvalues
113: !
114: subroutine STBackTransform_User(st, n, eigr, eigi, ierr)
115: use slepceps
116: implicit none
118: ST :: st
119: PetscInt :: n, j
120: PetscScalar :: eigr(*), eigi(*)
121: PetscErrorCode, intent(out) :: ierr
123: do j = 1, n
124: eigr(j) = 1.0/eigr(j)
125: end do
126: ierr = 0
127: end subroutine
129: end module ex10fmodule
131: ! ----------------------------------------------------------------------
133: program ex10f
134: use slepceps
135: use ex10fmodule
136: implicit none
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: ! Declarations
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: !
142: Mat :: A ! operator matrix
143: EPS :: eps ! eigenproblem solver context
144: ST :: st
145: EPSType :: tname
146: PetscInt :: n, i, Istart, Iend
147: PetscInt :: nev, row(1), col(3)
148: PetscScalar :: val(3)
149: PetscBool :: flg, isShell, terse
150: PetscMPIInt :: rank
151: PetscErrorCode :: ierr
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: ! Beginning of program
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
158: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
159: n = 30
160: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
162: if (rank == 0) then
163: write (*, '(/a,i6/)') '1-D Laplacian Eigenproblem (shell-enabled), n=', n
164: end if
166: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: ! Compute the operator matrix that defines the eigensystem, Ax=kx
168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
171: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
172: PetscCallA(MatSetFromOptions(A, ierr))
174: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
175: if (Istart == 0) then
176: row(1) = 0
177: col(1) = 0
178: col(2) = 1
179: val(1) = 2.0
180: val(2) = -1.0
181: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, row, 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
182: Istart = Istart + 1
183: end if
184: if (Iend == n) then
185: row(1) = n - 1
186: col(1) = n - 2
187: col(2) = n - 1
188: val(1) = -1.0
189: val(2) = 2.0
190: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, row, 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
191: Iend = Iend - 1
192: end if
193: val(1) = -1.0
194: val(2) = 2.0
195: val(3) = -1.0
196: do i = Istart, Iend - 1
197: row(1) = i
198: col(1) = i - 1
199: col(2) = i
200: col(3) = i + 1
201: PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, row, 3_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
202: end do
204: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
205: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
207: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: ! Create the eigensolver and set various options
209: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211: ! ** Create eigensolver context
212: PetscCallA(EPSCreate(PETSC_COMM_WORLD, eps, ierr))
214: ! ** Set operators. In this case, it is a standard eigenvalue problem
215: PetscCallA(EPSSetOperators(eps, A, PETSC_NULL_MAT, ierr))
216: PetscCallA(EPSSetProblemType(eps, EPS_NHEP, ierr))
218: ! ** Set solver parameters at runtime
219: PetscCallA(EPSSetFromOptions(eps, ierr))
221: ! ** Initialize shell spectral transformation if selected by user
222: PetscCallA(EPSGetST(eps, st, ierr))
223: PetscCallA(PetscObjectTypeCompare(st, STSHELL, isShell, ierr))
225: if (isShell) then
226: ! ** Change sorting criterion since this ST example computes values
227: ! ** closest to 0
228: PetscCallA(EPSSetWhichEigenpairs(eps, EPS_SMALLEST_REAL, ierr))
230: ! ** In Fortran, instead of a context for the user-defined spectral transform
231: ! ** we use a module containing any application-specific data, initialized here
232: PetscCallA(KSPCreate(PETSC_COMM_WORLD, myksp, ierr))
233: PetscCallA(KSPAppendOptionsPrefix(myksp, "st_", ierr))
235: ! ** (Required) Set the user-defined routine for applying the operator
236: PetscCallA(STShellSetApply(st, STApply_User, ierr))
238: ! ** (Optional) Set the user-defined routine for applying the transposed operator
239: PetscCallA(STShellSetApplyTranspose(st, STApplyTranspose_User, ierr))
241: #if defined(PETSC_USE_COMPLEX)
242: ! ** (Optional) Set the user-defined routine for applying the conjugate-transposed operator
243: PetscCallA(STShellSetApplyHermitianTranspose(st, STApplyHermitianTranspose_User, ierr))
244: #endif
246: ! ** (Optional) Set the user-defined routine for back-transformation
247: PetscCallA(STShellSetBackTransform(st, STBackTransform_User, ierr))
249: ! ** (Optional) Set a name for the transformation, used for STView()
250: PetscCallA(PetscObjectSetName(st, 'MyTransformation', ierr))
252: ! ** (Optional) Do any setup required for the new transformation
253: PetscCallA(KSPSetOperators(myksp, A, A, ierr))
254: PetscCallA(KSPSetFromOptions(myksp, ierr))
255: end if
257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: ! Solve the eigensystem
259: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
261: PetscCallA(EPSSolve(eps, ierr))
263: ! ** Optional: Get some information from the solver and display it
264: PetscCallA(EPSGetType(eps, tname, ierr))
265: if (rank == 0) then
266: write (*, '(a,a/)') ' Solution method: ', tname
267: end if
268: PetscCallA(EPSGetDimensions(eps, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
269: if (rank == 0) then
270: write (*, '(a,i2)') ' Number of requested eigenvalues:', nev
271: end if
273: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: ! Display solution and clean up
275: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: ! ** show detailed info unless -terse option is given by user
278: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
279: if (terse) then
280: PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_NULL_VIEWER, ierr))
281: else
282: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
283: PetscCallA(EPSConvergedReasonView(eps, PETSC_VIEWER_STDOUT_WORLD, ierr))
284: PetscCallA(EPSErrorView(eps, EPS_ERROR_RELATIVE, PETSC_VIEWER_STDOUT_WORLD, ierr))
285: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
286: end if
287: if (isShell) then
288: PetscCallA(KSPDestroy(myksp, ierr))
289: end if
290: PetscCallA(EPSDestroy(eps, ierr))
291: PetscCallA(MatDestroy(A, ierr))
292: PetscCallA(SlepcFinalize(ierr))
293: end program ex10f
295: !/*TEST
296: !
297: ! testset:
298: ! args: -eps_nev 5 -eps_non_hermitian -terse
299: ! output_file: output/ex10_1.out
300: ! requires: !single
301: ! test:
302: ! suffix: 1_sinvert
303: ! args: -st_type sinvert
304: ! test:
305: ! suffix: 1_shell
306: ! args: -st_type shell -eps_two_sided {{0 1}}
307: !
308: !TEST*/