Actual source code: ex10f.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> ./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: !
19: ! Note: this example illustrates old error checking with CHKERRA instead
20: ! of PetscCallA()
21: ! ----------------------------------------------------------------------
22: !
23: ! Module contains data needed by shell ST
24: !
25: module mymoduleex10f
26: #include <slepc/finclude/slepceps.h>
27: use slepceps
28: implicit none
30: KSP myksp
31: end module
33: ! ----------------------------------------------------------------------
35: program main
36: #include <slepc/finclude/slepceps.h>
37: use slepceps
38: use mymoduleex10f
39: implicit none
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: ! Declarations
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44: !
45: ! Variables:
46: ! A operator matrix
47: ! eps eigenproblem solver context
49: Mat A
50: EPS eps
51: ST st
52: EPSType tname
53: PetscInt n, i, Istart, Iend, one, two, three
54: PetscInt nev, row(1), col(3)
55: PetscScalar val(3)
56: PetscBool flg, isShell, terse
57: PetscMPIInt rank
58: PetscErrorCode ierr
60: ! Note: Any user-defined Fortran routines MUST be declared as external.
61: external STApply_User, STApplyTranspose_User, STApplyHermitianTranspose_User, STBackTransform_User
63: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: ! Beginning of program
65: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: one = 1
68: two = 2
69: three = 3
70: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
71: if (ierr .ne. 0) then
72: print*,'SlepcInitialize failed'
73: stop
74: endif
75: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRMPIA(ierr)
76: n = 30
77: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(ierr)
79: if (rank .eq. 0) then
80: write(*,'(/A,I6/)') '1-D Laplacian Eigenproblem (shell-enabled), n=',n
81: endif
83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: ! Compute the operator matrix that defines the eigensystem, Ax=kx
85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
88: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
89: call MatSetFromOptions(A,ierr);CHKERRA(ierr)
91: call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
92: if (Istart .eq. 0) then
93: row(1) = 0
94: col(1) = 0
95: col(2) = 1
96: val(1) = 2.0
97: val(2) = -1.0
98: call MatSetValues(A,one,row,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
99: Istart = Istart+1
100: endif
101: if (Iend .eq. n) then
102: row(1) = n-1
103: col(1) = n-2
104: col(2) = n-1
105: val(1) = -1.0
106: val(2) = 2.0
107: call MatSetValues(A,one,row,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
108: Iend = Iend-1
109: endif
110: val(1) = -1.0
111: val(2) = 2.0
112: val(3) = -1.0
113: do i=Istart,Iend-1
114: row(1) = i
115: col(1) = i-1
116: col(2) = i
117: col(3) = i+1
118: call MatSetValues(A,one,row,three,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
119: enddo
121: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
122: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Create the eigensolver and set various options
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: ! ** Create eigensolver context
129: call EPSCreate(PETSC_COMM_WORLD,eps,ierr);CHKERRA(ierr)
131: ! ** Set operators. In this case, it is a standard eigenvalue problem
132: call EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr);CHKERRA(ierr)
133: call EPSSetProblemType(eps,EPS_NHEP,ierr);CHKERRA(ierr)
135: ! ** Set solver parameters at runtime
136: call EPSSetFromOptions(eps,ierr);CHKERRA(ierr)
138: ! ** Initialize shell spectral transformation if selected by user
139: call EPSGetST(eps,st,ierr);CHKERRA(ierr)
140: call PetscObjectTypeCompare(st,STSHELL,isShell,ierr);CHKERRA(ierr)
142: if (isShell) then
143: ! ** Change sorting criterion since this ST example computes values
144: ! ** closest to 0
145: call EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL,ierr);CHKERRA(ierr)
147: ! ** In Fortran, instead of a context for the user-defined spectral transform
148: ! ** we use a module containing any application-specific data, initialized here
149: call KSPCreate(PETSC_COMM_WORLD,myksp,ierr);CHKERRA(ierr)
150: call KSPAppendOptionsPrefix(myksp,"st_",ierr);CHKERRA(ierr)
152: ! ** (Required) Set the user-defined routine for applying the operator
153: call STShellSetApply(st,STApply_User,ierr);CHKERRA(ierr)
155: ! ** (Optional) Set the user-defined routine for applying the transposed operator
156: call STShellSetApplyTranspose(st,STApplyTranspose_User,ierr);CHKERRA(ierr)
158: #if defined(PETSC_USE_COMPLEX)
159: ! ** (Optional) Set the user-defined routine for applying the conjugate-transposed operator
160: call STShellSetApplyHermitianTranspose(st,STApplyHermitianTranspose_User,ierr);CHKERRA(ierr)
161: #endif
163: ! ** (Optional) Set the user-defined routine for back-transformation
164: call STShellSetBackTransform(st,STBackTransform_User,ierr);CHKERRA(ierr)
166: ! ** (Optional) Set a name for the transformation, used for STView()
167: call PetscObjectSetName(st,'MyTransformation',ierr);CHKERRA(ierr)
169: ! ** (Optional) Do any setup required for the new transformation
170: call KSPSetOperators(myksp,A,A,ierr);CHKERRA(ierr)
171: call KSPSetFromOptions(myksp,ierr);CHKERRA(ierr)
172: endif
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: ! Solve the eigensystem
176: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: call EPSSolve(eps,ierr);CHKERRA(ierr)
180: ! ** Optional: Get some information from the solver and display it
181: call EPSGetType(eps,tname,ierr);CHKERRA(ierr)
182: if (rank .eq. 0) then
183: write(*,'(A,A,/)') ' Solution method: ', tname
184: endif
185: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
186: if (rank .eq. 0) then
187: write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
188: endif
190: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191: ! Display solution and clean up
192: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: ! ** show detailed info unless -terse option is given by user
195: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
196: if (terse) then
197: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
198: else
199: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
200: call EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
201: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
202: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
203: endif
204: if (isShell) then
205: call KSPDestroy(myksp,ierr);CHKERRA(ierr)
206: endif
207: call EPSDestroy(eps,ierr);CHKERRA(ierr)
208: call MatDestroy(A,ierr);CHKERRA(ierr)
209: call SlepcFinalize(ierr)
210: end
212: ! -------------------------------------------------------------------
213: !
214: ! STApply_User - This routine demonstrates the use of a user-provided spectral
215: ! transformation. The transformation implemented in this code is just OP=A^-1.
216: !
217: ! Input Parameters:
218: ! st - spectral transformation context
219: ! x - input vector
220: !
221: ! Output Parameter:
222: ! y - output vector
223: !
224: subroutine STApply_User(st,x,y,ierr)
225: #include <slepc/finclude/slepceps.h>
226: use slepceps
227: use mymoduleex10f
228: implicit none
230: ST st
231: Vec x,y
232: PetscErrorCode ierr
234: call KSPSolve(myksp,x,y,ierr);CHKERRQ(ierr)
236: end
238: ! -------------------------------------------------------------------
239: !
240: ! STApplyTranspose_User - This is not required unless using a two-sided eigensolver
241: !
242: ! Input Parameters:
243: ! st - spectral transformation context
244: ! x - input vector
245: !
246: ! Output Parameter:
247: ! y - output vector
248: !
249: subroutine STApplyTranspose_User(st,x,y,ierr)
250: #include <slepc/finclude/slepceps.h>
251: use slepceps
252: use mymoduleex10f
253: implicit none
255: ST st
256: Vec x,y
257: PetscErrorCode ierr
259: call KSPSolveTranspose(myksp,x,y,ierr);CHKERRQ(ierr)
261: end
263: #if defined(PETSC_USE_COMPLEX)
264: ! -------------------------------------------------------------------
265: !
266: ! STApplyHermitianTranspose_User - This is not required unless using a two-sided eigensolver
267: ! in complex scalars
268: !
269: ! Input Parameters:
270: ! st - spectral transformation context
271: ! x - input vector
272: !
273: ! Output Parameter:
274: ! y - output vector
275: !
276: subroutine STApplyHermitianTranspose_User(st,x,y,ierr)
277: #include <slepc/finclude/slepceps.h>
278: use slepceps
279: use mymoduleex10f
280: implicit none
282: ST st
283: Vec x,y,w
284: PetscErrorCode ierr
286: call VecDuplicate(x,w,ierr);CHKERRQ(ierr)
287: call VecCopy(x,w,ierr);CHKERRQ(ierr)
288: call VecConjugate(w,ierr);CHKERRQ(ierr)
289: call KSPSolveTranspose(myksp,w,y,ierr);CHKERRQ(ierr)
290: call VecConjugate(y,ierr);CHKERRQ(ierr)
291: call VecDestroy(w,ierr);CHKERRQ(ierr)
293: end
294: #endif
296: ! -------------------------------------------------------------------
297: !
298: ! STBackTransform_User - This routine demonstrates the use of a user-provided spectral
299: ! transformation
300: !
301: ! Input Parameters:
302: ! st - spectral transformation context
303: ! n - number of eigenvalues to transform
304: !
305: ! Output Parameters:
306: ! eigr - real part of eigenvalues
307: ! eigi - imaginary part of eigenvalues
308: !
309: subroutine STBackTransform_User(st,n,eigr,eigi,ierr)
310: #include <slepc/finclude/slepceps.h>
311: use slepceps
312: use mymoduleex10f
313: implicit none
315: ST st
316: PetscInt n, j
317: PetscScalar eigr(*), eigi(*)
318: PetscErrorCode ierr
320: do j=1,n
321: eigr(j) = 1.0 / eigr(j)
322: enddo
323: ierr = 0
325: end
327: !/*TEST
328: !
329: ! testset:
330: ! args: -eps_nev 5 -eps_non_hermitian -terse
331: ! output_file: output/ex10_1.out
332: ! requires: !single
333: ! test:
334: ! suffix: 1_sinvert
335: ! args: -st_type sinvert
336: ! test:
337: ! suffix: 1_shell
338: ! args: -st_type shell
339: !
340: !TEST*/