Actual source code: ex10f90.F90

slepc-3.20.2 2024-03-15
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 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*/