Actual source code: ex10f90.F90
slepc-3.20.2 2024-03-15
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 mymoduleex10f90
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 mymoduleex10f90
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, 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)
90: call MatSetUp(A,ierr);CHKERRA(ierr)
92: call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
93: if (Istart .eq. 0) then
94: row(1) = 0
95: col(1) = 0
96: col(2) = 1
97: val(1) = 2.0
98: val(2) = -1.0
99: call MatSetValues(A,one,row,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
100: Istart = Istart+1
101: endif
102: if (Iend .eq. n) then
103: row(1) = n-1
104: col(1) = n-2
105: col(2) = n-1
106: val(1) = -1.0
107: val(2) = 2.0
108: call MatSetValues(A,one,row,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
109: Iend = Iend-1
110: endif
111: val(1) = -1.0
112: val(2) = 2.0
113: val(3) = -1.0
114: do i=Istart,Iend-1
115: row(1) = i
116: col(1) = i-1
117: col(2) = i
118: col(3) = i+1
119: call MatSetValues(A,one,row,three,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
120: enddo
122: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
123: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: ! Create the eigensolver and set various options
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: ! ** Create eigensolver context
130: call EPSCreate(PETSC_COMM_WORLD,eps,ierr);CHKERRA(ierr)
132: ! ** Set operators. In this case, it is a standard eigenvalue problem
133: call EPSSetOperators(eps,A,PETSC_NULL_MAT,ierr);CHKERRA(ierr)
134: call EPSSetProblemType(eps,EPS_NHEP,ierr);CHKERRA(ierr)
136: ! ** Set solver parameters at runtime
137: call EPSSetFromOptions(eps,ierr);CHKERRA(ierr)
139: ! ** Initialize shell spectral transformation if selected by user
140: call EPSGetST(eps,st,ierr);CHKERRA(ierr)
141: call PetscObjectTypeCompare(st,STSHELL,isShell,ierr);CHKERRA(ierr)
143: if (isShell) then
144: ! ** Change sorting criterion since this ST example computes values
145: ! ** closest to 0
146: call EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL,ierr);CHKERRA(ierr)
148: ! ** In Fortran, instead of a context for the user-defined spectral transform
149: ! ** we use a module containing any application-specific data, initialized here
150: call KSPCreate(PETSC_COMM_WORLD,myksp,ierr);CHKERRA(ierr)
151: call KSPAppendOptionsPrefix(myksp,"st_",ierr);CHKERRA(ierr)
153: ! ** (Required) Set the user-defined routine for applying the operator
154: call STShellSetApply(st,STApply_User,ierr);CHKERRA(ierr)
156: ! ** (Optional) Set the user-defined routine for applying the transposed operator
157: call STShellSetApplyTranspose(st,STApplyTranspose_User,ierr);CHKERRA(ierr)
159: ! ** (Optional) Set the user-defined routine for back-transformation
160: call STShellSetBackTransform(st,STBackTransform_User,ierr);CHKERRA(ierr)
162: ! ** (Optional) Set a name for the transformation, used for STView()
163: call PetscObjectSetName(st,'MyTransformation',ierr);CHKERRA(ierr)
165: ! ** (Optional) Do any setup required for the new transformation
166: call KSPSetOperators(myksp,A,A,ierr);CHKERRA(ierr)
167: call KSPSetFromOptions(myksp,ierr);CHKERRA(ierr)
168: endif
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: ! Solve the eigensystem
172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: call EPSSolve(eps,ierr);CHKERRA(ierr)
176: ! ** Optional: Get some information from the solver and display it
177: call EPSGetType(eps,tname,ierr);CHKERRA(ierr)
178: if (rank .eq. 0) then
179: write(*,'(A,A,/)') ' Solution method: ', tname
180: endif
181: call EPSGetDimensions(eps,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
182: if (rank .eq. 0) then
183: write(*,'(A,I2)') ' Number of requested eigenvalues:',nev
184: endif
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: ! Display solution and clean up
188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: ! ** show detailed info unless -terse option is given by user
191: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
192: if (terse) then
193: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
194: else
195: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
196: call EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
197: call EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
198: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
199: endif
200: if (isShell) then
201: call KSPDestroy(myksp,ierr);CHKERRA(ierr)
202: endif
203: call EPSDestroy(eps,ierr);CHKERRA(ierr)
204: call MatDestroy(A,ierr);CHKERRA(ierr)
205: call SlepcFinalize(ierr)
206: end
208: ! -------------------------------------------------------------------
209: !
210: ! STApply_User - This routine demonstrates the use of a user-provided spectral
211: ! transformation. The transformation implemented in this code is just OP=A^-1.
212: !
213: ! Input Parameters:
214: ! st - spectral transformation context
215: ! x - input vector
216: !
217: ! Output Parameter:
218: ! y - output vector
219: !
220: subroutine STApply_User(st,x,y,ierr)
221: #include <slepc/finclude/slepceps.h>
222: use slepceps
223: use mymoduleex10f90
224: implicit none
226: ST st
227: Vec x,y
228: PetscErrorCode ierr
230: call KSPSolve(myksp,x,y,ierr);CHKERRQ(ierr)
232: return
233: end
235: ! -------------------------------------------------------------------
236: !
237: ! STApplyTranspose_User - This is not required unless using a two-sided eigensolver
238: !
239: ! Input Parameters:
240: ! st - spectral transformation context
241: ! x - input vector
242: !
243: ! Output Parameter:
244: ! y - output vector
245: !
246: subroutine STApplyTranspose_User(st,x,y,ierr)
247: #include <slepc/finclude/slepceps.h>
248: use slepceps
249: use mymoduleex10f90
250: implicit none
252: ST st
253: Vec x,y
254: PetscErrorCode ierr
256: call KSPSolveTranspose(myksp,x,y,ierr);CHKERRQ(ierr)
258: return
259: end
261: ! -------------------------------------------------------------------
262: !
263: ! STBackTransform_User - This routine demonstrates the use of a user-provided spectral
264: ! transformation
265: !
266: ! Input Parameters:
267: ! st - spectral transformation context
268: ! n - number of eigenvalues to transform
269: !
270: ! Output Parameters:
271: ! eigr - real part of eigenvalues
272: ! eigi - imaginary part of eigenvalues
273: !
274: subroutine STBackTransform_User(st,n,eigr,eigi,ierr)
275: #include <slepc/finclude/slepceps.h>
276: use slepceps
277: use mymoduleex10f90
278: implicit none
280: ST st
281: PetscInt n, j
282: PetscScalar eigr(*), eigi(*)
283: PetscErrorCode ierr
285: do j=1,n
286: eigr(j) = 1.0 / eigr(j)
287: enddo
288: ierr = 0
290: return
291: end
293: !/*TEST
294: !
295: ! testset:
296: ! args: -eps_nev 5 -eps_non_hermitian -terse
297: ! output_file: output/ex10_1.out
298: ! requires: !single
299: ! test:
300: ! suffix: 1_sinvert
301: ! args: -st_type sinvert
302: ! test:
303: ! suffix: 1_shell
304: ! args: -st_type shell
305: !
306: !TEST*/