Actual source code: ex20f90.F90
slepc-3.18.2 2023-01-26
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*/