Actual source code: ex20f90.F90

slepc-3.20.0 2023-09-29
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> ./ex20f90 [-n <n>] [SLEPc opts]
 11: !
 12: !  Description: Simple 1-D nonlinear eigenproblem. Fortran90 equivalent of ex20.c
 13: !
 14: !  The command line options are:
 15: !    -n <n>, where <n> = number of grid subdivisions
 16: !
 17: ! ----------------------------------------------------------------------
 18: !  Solve 1-D PDE
 19: !           -u'' = lambda*u
 20: !  on [0,1] subject to
 21: !           u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
 22: ! ----------------------------------------------------------------------
 23: !

 25: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 26: !     User-defined application context
 27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 28:       module UserModule
 29: #include <slepc/finclude/slepcnep.h>
 30:       use slepcnep
 31:       type User
 32:         PetscScalar kappa
 33:         PetscReal   h
 34:       end type User
 35:       end module

 37:       program main
 38: #include <slepc/finclude/slepcnep.h>
 39:       use UserModule
 40:       implicit none

 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !     Declarations
 44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45: !
 46: !  Variables:
 47: !     nep       nonlinear eigensolver context
 48: !     x         eigenvector
 49: !     lambda    eigenvalue
 50: !     F,J       Function and Jacobian matrices
 51: !     ctx       user-defined context

 53:       NEP            nep
 54:       Vec            x, v(1)
 55:       PetscScalar    lambda
 56:       Mat            F, J
 57:       type(User)     ctx
 58:       NEPType        tname
 59:       PetscInt       n, i, k, nev, its, maxit, nconv, three, one
 60:       PetscReal      tol, norm
 61:       PetscScalar    alpha
 62:       PetscMPIInt    rank
 63:       PetscBool      flg
 64:       PetscErrorCode ierr
 65: !  Note: Any user-defined Fortran routines (such as FormJacobian)
 66: !  MUST be declared as external.
 67:       external       FormFunction, FormJacobian

 69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 70: !     Beginning of program
 71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 73:       PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
 74:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
 75:       n = 128
 76:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
 77:       if (rank .eq. 0) then
 78:         write(*,'(/A,I4)') 'Nonlinear Eigenproblem, n =',n
 79:       endif

 81:       ctx%h = 1.0/real(n)
 82:       ctx%kappa = 1.0

 84:       three = 3
 85:       one = 1

 87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88: !     Create matrix data structure to hold the Function and the Jacobian
 89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 91:       PetscCallA(MatCreate(PETSC_COMM_WORLD,F,ierr))
 92:       PetscCallA(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
 93:       PetscCallA(MatSetFromOptions(F,ierr))
 94:       PetscCallA(MatSeqAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,ierr))
 95:       PetscCallA(MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr))
 96:       PetscCallA(MatSetUp(F,ierr))

 98:       PetscCallA(MatCreate(PETSC_COMM_WORLD,J,ierr))
 99:       PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
100:       PetscCallA(MatSetFromOptions(J,ierr))
101:       PetscCallA(MatSeqAIJSetPreallocation(J,three,PETSC_NULL_INTEGER,ierr))
102:       PetscCallA(MatMPIAIJSetPreallocation(J,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr))
103:       PetscCallA(MatSetUp(J,ierr))

105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: !     Create the eigensolver and set various options
107: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

109: !     ** Create eigensolver context
110:       PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))

112: !     ** Set routines for evaluation of Function and Jacobian
113:       PetscCallA(NEPSetFunction(nep,F,F,FormFunction,ctx,ierr))
114:       PetscCallA(NEPSetJacobian(nep,J,FormJacobian,ctx,ierr))

116: !     ** Customize nonlinear solver
117:       tol = 1e-9
118:       PetscCallA(NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr))
119:       k = 1
120:       PetscCallA(NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr))

122: !     ** Set solver parameters at runtime
123:       PetscCallA(NEPSetFromOptions(nep,ierr))

125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: !     Solve the eigensystem
127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

129: !     ** Evaluate initial guess
130:       PetscCallA(MatCreateVecs(F,x,PETSC_NULL_VEC,ierr))
131:       PetscCallA(VecDuplicate(x,v(1),ierr))
132:       alpha = 1.0
133:       PetscCallA(VecSet(v(1),alpha,ierr))
134:       k = 1
135:       PetscCallA(NEPSetInitialSpace(nep,k,v,ierr))
136:       PetscCallA(VecDestroy(v(1),ierr))

138: !     ** Call the solver
139:       PetscCallA(NEPSolve(nep,ierr))
140:       PetscCallA(NEPGetIterationNumber(nep,its,ierr))
141:       if (rank .eq. 0) then
142:         write(*,'(A,I3)') ' Number of NEP iterations =',its
143:       endif

145: !     ** Optional: Get some information from the solver and display it
146:       PetscCallA(NEPGetType(nep,tname,ierr))
147:       if (rank .eq. 0) then
148:         write(*,'(A,A10)') ' Solution method: ',tname
149:       endif
150:       PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
151:       if (rank .eq. 0) then
152:         write(*,'(A,I4)') ' Number of requested eigenvalues:',nev
153:       endif
154:       PetscCallA(NEPGetTolerances(nep,tol,maxit,ierr))
155:       if (rank .eq. 0) then
156:         write(*,'(A,F12.9,A,I5)') ' Stopping condition: tol=',tol,', maxit=',maxit
157:       endif

159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: !     Display solution and clean up
161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

163:       PetscCallA(NEPGetConverged(nep,nconv,ierr))
164:       if (rank .eq. 0) then
165:         write(*,'(A,I2/)') ' Number of converged approximate eigenpairs:',nconv
166:       endif

168: !     ** Display eigenvalues and relative errors
169:       if (nconv .gt. 0) then
170:         if (rank .eq. 0) then
171:           write(*,*) '        k              ||T(k)x||'
172:           write(*,*) '----------------- ------------------'
173:         endif
174:         do i=0,nconv-1
175: !         ** Get converged eigenpairs: (in this example they are always real)
176:           PetscCallA(NEPGetEigenpair(nep,i,lambda,PETSC_NULL_SCALAR,x,PETSC_NULL_VEC,ierr))

178: !         ** Compute residual norm and error
179:           PetscCallA(NEPComputeError(nep,i,NEP_ERROR_RELATIVE,norm,ierr))
180:           if (rank .eq. 0) then
181:             write(*,'(1P,E15.4,E18.4)') PetscRealPart(lambda), norm
182:           endif
183:         enddo
184:         if (rank .eq. 0) then
185:           write(*,*)
186:         endif
187:       endif

189:       PetscCallA(NEPDestroy(nep,ierr))
190:       PetscCallA(MatDestroy(F,ierr))
191:       PetscCallA(MatDestroy(J,ierr))
192:       PetscCallA(VecDestroy(x,ierr))
193:       PetscCallA(SlepcFinalize(ierr))
194:       end

196: ! ---------------  Evaluate Function matrix  T(lambda)  ----------------

198:       subroutine FormFunction(nep,lambda,fun,B,ctx,ierr)
199:       use UserModule
200:       implicit none
201:       NEP            nep
202:       PetscScalar    lambda, A(3), c, d
203:       Mat            fun,B
204:       type(User)     ctx
205:       PetscReal      h
206:       PetscInt       i, n, j(3), Istart, Iend, one, two, three
207:       PetscErrorCode ierr

209: !     ** Compute Function entries and insert into matrix
210:       PetscCall(MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr))
211:       PetscCall(MatGetOwnershipRange(fun,Istart,Iend,ierr))
212:       h = ctx%h
213:       c = ctx%kappa/(lambda-ctx%kappa)
214:       d = n
215:       one = 1
216:       two = 2
217:       three = 3

219: !     ** Boundary points
220:       if (Istart .eq. 0) then
221:         i = 0
222:         j(1) = 0
223:         j(2) = 1
224:         A(1) = 2.0*(d-lambda*h/3.0)
225:         A(2) = -d-lambda*h/6.0
226:         PetscCall(MatSetValues(fun,one,i,two,j,A,INSERT_VALUES,ierr))
227:         Istart = Istart + 1
228:       endif

230:       if (Iend .eq. n) then
231:         i = n-1
232:         j(1) = n-2
233:         j(2) = n-1
234:         A(1) = -d-lambda*h/6.0
235:         A(2) = d-lambda*h/3.0+lambda*c
236:         PetscCall(MatSetValues(fun,one,i,two,j,A,INSERT_VALUES,ierr))
237:         Iend = Iend - 1
238:       endif

240: !     ** Interior grid points
241:       do i=Istart,Iend-1
242:         j(1) = i-1
243:         j(2) = i
244:         j(3) = i+1
245:         A(1) = -d-lambda*h/6.0
246:         A(2) = 2.0*(d-lambda*h/3.0)
247:         A(3) = -d-lambda*h/6.0
248:         PetscCall(MatSetValues(fun,one,i,three,j,A,INSERT_VALUES,ierr))
249:       enddo

251: !     ** Assemble matrix
252:       PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr))
253:       PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr))
254:       return
255:       end

257: ! ---------------  Evaluate Jacobian matrix  T'(lambda)  ---------------

259:       subroutine FormJacobian(nep,lambda,jac,ctx,ierr)
260:       use UserModule
261:       implicit none
262:       NEP            nep
263:       PetscScalar    lambda, A(3), c
264:       Mat            jac
265:       type(User)     ctx
266:       PetscReal      h
267:       PetscInt       i, n, j(3), Istart, Iend, one, two, three
268:       PetscErrorCode ierr

270: !     ** Compute Jacobian entries and insert into matrix
271:       PetscCall(MatGetSize(jac,n,PETSC_NULL_INTEGER,ierr))
272:       PetscCall(MatGetOwnershipRange(jac,Istart,Iend,ierr))
273:       h = ctx%h
274:       c = ctx%kappa/(lambda-ctx%kappa)
275:       one = 1
276:       two = 2
277:       three = 3

279: !     ** Boundary points
280:       if (Istart .eq. 0) then
281:         i = 0
282:         j(1) = 0
283:         j(2) = 1
284:         A(1) = -2.0*h/3.0
285:         A(2) = -h/6.0
286:         PetscCall(MatSetValues(jac,one,i,two,j,A,INSERT_VALUES,ierr))
287:         Istart = Istart + 1
288:       endif

290:       if (Iend .eq. n) then
291:         i = n-1
292:         j(1) = n-2
293:         j(2) = n-1
294:         A(1) = -h/6.0
295:         A(2) = -h/3.0-c*c
296:         PetscCall(MatSetValues(jac,one,i,two,j,A,INSERT_VALUES,ierr))
297:         Iend = Iend - 1
298:       endif

300: !     ** Interior grid points
301:       do i=Istart,Iend-1
302:         j(1) = i-1
303:         j(2) = i
304:         j(3) = i+1
305:         A(1) = -h/6.0
306:         A(2) = -2.0*h/3.0
307:         A(3) = -h/6.0
308:         PetscCall(MatSetValues(jac,one,i,three,j,A,INSERT_VALUES,ierr))
309:       enddo

311: !     ** Assemble matrix
312:       PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
313:       PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
314:       return
315:       end

317: !/*TEST
318: !
319: !   test:
320: !      suffix: 1
321: !      args: -nep_target 4
322: !      filter: sed -e "s/[0-9]\.[0-9]*E-[0-9]*/removed/g" -e "s/ Number of NEP iterations = [ 0-9]*/ Number of NEP iterations = /"
323: !      requires: !single
324: !
325: !TEST*/