Actual source code: ex20f90.F90

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2019, 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
 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:       alpha = 1.0
136:       call VecSet(x,alpha,ierr);CHKERRA(ierr)
137:       k = 1
138:       call NEPSetInitialSpace(nep,k,x,ierr);CHKERRA(ierr)

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

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

161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: !     Display solution and clean up
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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