Actual source code: ex20f90.F90

slepc-3.17.1 2022-04-11
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:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 74:       if (ierr .ne. 0) then
 75:         print*,'SlepcInitialize failed'
 76:         stop
 77:       endif
 78:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
 79:       n = 128
 80:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(ierr)
 81:       if (rank .eq. 0) then
 82:         write(*,'(/A,I4)') 'Nonlinear Eigenproblem, n =',n
 83:       endif

 85:       ctx%h = 1.0/real(n)
 86:       ctx%kappa = 1.0

 88:       three = 3
 89:       one = 1

 91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92: !     Create matrix data structure to hold the Function and the Jacobian
 93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 95:       call MatCreate(PETSC_COMM_WORLD,F,ierr);CHKERRA(ierr)
 96:       call MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
 97:       call MatSetFromOptions(F,ierr);CHKERRA(ierr)
 98:       call MatSeqAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
 99:       call MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
100:       call MatSetUp(F,ierr);CHKERRA(ierr)

102:       call MatCreate(PETSC_COMM_WORLD,J,ierr);CHKERRA(ierr)
103:       call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
104:       call MatSetFromOptions(J,ierr);CHKERRA(ierr)
105:       call MatSeqAIJSetPreallocation(J,three,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
106:       call MatMPIAIJSetPreallocation(J,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
107:       call MatSetUp(J,ierr);CHKERRA(ierr)

109: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: !     Create the eigensolver and set various options
111: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

113: !     ** Create eigensolver context
114:       call NEPCreate(PETSC_COMM_WORLD,nep,ierr);CHKERRA(ierr)

116: !     ** Set routines for evaluation of Function and Jacobian
117:       call NEPSetFunction(nep,F,F,FormFunction,ctx,ierr);CHKERRA(ierr)
118:       call NEPSetJacobian(nep,J,FormJacobian,ctx,ierr);CHKERRA(ierr)

120: !     ** Customize nonlinear solver
121:       tol = 1e-9
122:       call NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
123:       k = 1
124:       call NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)

126: !     ** Set solver parameters at runtime
127:       call NEPSetFromOptions(nep,ierr);CHKERRA(ierr)

129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: !     Solve the eigensystem
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

133: !     ** Evaluate initial guess
134:       call MatCreateVecs(F,x,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
135:       call VecDuplicate(x,v(1),ierr);CHKERRA(ierr)
136:       alpha = 1.0
137:       call VecSet(v(1),alpha,ierr);CHKERRA(ierr)
138:       k = 1
139:       call NEPSetInitialSpace(nep,k,v,ierr);CHKERRA(ierr)
140:       call VecDestroy(v(1),ierr);CHKERRA(ierr)

142: !     ** Call the solver
143:       call NEPSolve(nep,ierr);CHKERRA(ierr)
144:       call NEPGetIterationNumber(nep,its,ierr);CHKERRA(ierr)
145:       if (rank .eq. 0) then
146:         write(*,'(A,I3)') ' Number of NEP iterations =',its
147:       endif

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

163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: !     Display solution and clean up
165: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

167:       call NEPGetConverged(nep,nconv,ierr);CHKERRA(ierr)
168:       if (rank .eq. 0) then
169:         write(*,'(A,I2/)') ' Number of converged approximate eigenpairs:',nconv
170:       endif

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

182: !         ** Compute residual norm and error
183:           call NEPComputeError(nep,i,NEP_ERROR_RELATIVE,norm,ierr);CHKERRA(ierr)
184:           if (rank .eq. 0) then
185:             write(*,'(1P,E15.4,E18.4)') PetscRealPart(lambda), norm
186:           endif
187:         enddo
188:         if (rank .eq. 0) then
189:           write(*,*)
190:         endif
191:       endif

193:       call NEPDestroy(nep,ierr);CHKERRA(ierr)
194:       call MatDestroy(F,ierr);CHKERRA(ierr)
195:       call MatDestroy(J,ierr);CHKERRA(ierr)
196:       call VecDestroy(x,ierr);CHKERRA(ierr)
197:       call SlepcFinalize(ierr)
198:       end

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

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

213: !     ** Compute Function entries and insert into matrix
214:       call MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
215:       call MatGetOwnershipRange(fun,Istart,Iend,ierr);CHKERRQ(ierr)
216:       h = ctx%h
217:       c = ctx%kappa/(lambda-ctx%kappa)
218:       d = n
219:       one = 1
220:       two = 2
221:       three = 3

223: !     ** Boundary points
224:       if (Istart .eq. 0) then
225:         i = 0
226:         j(1) = 0
227:         j(2) = 1
228:         A(1) = 2.0*(d-lambda*h/3.0)
229:         A(2) = -d-lambda*h/6.0
230:         call MatSetValues(fun,one,i,two,j,A,INSERT_VALUES,ierr);CHKERRQ(ierr)
231:         Istart = Istart + 1
232:       endif

234:       if (Iend .eq. n) then
235:         i = n-1
236:         j(1) = n-2
237:         j(2) = n-1
238:         A(1) = -d-lambda*h/6.0
239:         A(2) = d-lambda*h/3.0+lambda*c
240:         call MatSetValues(fun,one,i,two,j,A,INSERT_VALUES,ierr);CHKERRQ(ierr)
241:         Iend = Iend - 1
242:       endif

244: !     ** Interior grid points
245:       do i=Istart,Iend-1
246:         j(1) = i-1
247:         j(2) = i
248:         j(3) = i+1
249:         A(1) = -d-lambda*h/6.0
250:         A(2) = 2.0*(d-lambda*h/3.0)
251:         A(3) = -d-lambda*h/6.0
252:         call MatSetValues(fun,one,i,three,j,A,INSERT_VALUES,ierr);CHKERRQ(ierr)
253:       enddo

255: !     ** Assemble matrix
256:       call MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
257:       call MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
258:       return
259:       end

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

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

274: !     ** Compute Jacobian entries and insert into matrix
275:       call MatGetSize(jac,n,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr)
276:       call MatGetOwnershipRange(jac,Istart,Iend,ierr);CHKERRQ(ierr)
277:       h = ctx%h
278:       c = ctx%kappa/(lambda-ctx%kappa)
279:       one = 1
280:       two = 2
281:       three = 3

283: !     ** Boundary points
284:       if (Istart .eq. 0) then
285:         i = 0
286:         j(1) = 0
287:         j(2) = 1
288:         A(1) = -2.0*h/3.0
289:         A(2) = -h/6.0
290:         call MatSetValues(jac,one,i,two,j,A,INSERT_VALUES,ierr);CHKERRQ(ierr)
291:         Istart = Istart + 1
292:       endif

294:       if (Iend .eq. n) then
295:         i = n-1
296:         j(1) = n-2
297:         j(2) = n-1
298:         A(1) = -h/6.0
299:         A(2) = -h/3.0-c*c
300:         call MatSetValues(jac,one,i,two,j,A,INSERT_VALUES,ierr);CHKERRQ(ierr)
301:         Iend = Iend - 1
302:       endif

304: !     ** Interior grid points
305:       do i=Istart,Iend-1
306:         j(1) = i-1
307:         j(2) = i
308:         j(3) = i+1
309:         A(1) = -h/6.0
310:         A(2) = -2.0*h/3.0
311:         A(3) = -h/6.0
312:         call MatSetValues(jac,one,i,three,j,A,INSERT_VALUES,ierr);CHKERRQ(ierr)
313:       enddo

315: !     ** Assemble matrix
316:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
317:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
318:       return
319:       end

321: !/*TEST
322: !
323: !   test:
324: !      suffix: 1
325: !      args: -nep_target 4
326: !      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 = /"
327: !      requires: !single
328: !
329: !TEST*/