Actual source code: ex10f.F90

slepc-3.22.1 2024-10-28
Report Typos and Errors
  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*/