Actual source code: ex20f.F90
slepc-main 2024-11-09
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> ./ex20f [-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_ARRAY,ierr))
95: PetscCallA(MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER_ARRAY,one,PETSC_NULL_INTEGER_ARRAY,ierr))
97: PetscCallA(MatCreate(PETSC_COMM_WORLD,J,ierr))
98: PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
99: PetscCallA(MatSetFromOptions(J,ierr))
100: PetscCallA(MatSeqAIJSetPreallocation(J,three,PETSC_NULL_INTEGER_ARRAY,ierr))
101: PetscCallA(MatMPIAIJSetPreallocation(J,three,PETSC_NULL_INTEGER_ARRAY,one,PETSC_NULL_INTEGER_ARRAY,ierr))
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: ! Create the eigensolver and set various options
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: ! ** Create eigensolver context
108: PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
110: ! ** Set routines for evaluation of Function and Jacobian
111: PetscCallA(NEPSetFunction(nep,F,F,FormFunction,ctx,ierr))
112: PetscCallA(NEPSetJacobian(nep,J,FormJacobian,ctx,ierr))
114: ! ** Customize nonlinear solver
115: tol = 1e-9
116: PetscCallA(NEPSetTolerances(nep,tol,PETSC_CURRENT_INTEGER,ierr))
117: k = 1
118: PetscCallA(NEPSetDimensions(nep,k,PETSC_DETERMINE_INTEGER,PETSC_DETERMINE_INTEGER,ierr))
120: ! ** Set solver parameters at runtime
121: PetscCallA(NEPSetFromOptions(nep,ierr))
123: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! Solve the eigensystem
125: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: ! ** Evaluate initial guess
128: PetscCallA(MatCreateVecs(F,x,PETSC_NULL_VEC,ierr))
129: PetscCallA(VecDuplicate(x,v(1),ierr))
130: alpha = 1.0
131: PetscCallA(VecSet(v(1),alpha,ierr))
132: k = 1
133: PetscCallA(NEPSetInitialSpace(nep,k,v,ierr))
134: PetscCallA(VecDestroy(v(1),ierr))
136: ! ** Call the solver
137: PetscCallA(NEPSolve(nep,ierr))
138: PetscCallA(NEPGetIterationNumber(nep,its,ierr))
139: if (rank .eq. 0) then
140: write(*,'(A,I3)') ' Number of NEP iterations =',its
141: endif
143: ! ** Optional: Get some information from the solver and display it
144: PetscCallA(NEPGetType(nep,tname,ierr))
145: if (rank .eq. 0) then
146: write(*,'(A,A10)') ' Solution method: ',tname
147: endif
148: PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
149: if (rank .eq. 0) then
150: write(*,'(A,I4)') ' Number of requested eigenvalues:',nev
151: endif
152: PetscCallA(NEPGetTolerances(nep,tol,maxit,ierr))
153: if (rank .eq. 0) then
154: write(*,'(A,F12.9,A,I5)') ' Stopping condition: tol=',tol,', maxit=',maxit
155: endif
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: ! Display solution and clean up
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: PetscCallA(NEPGetConverged(nep,nconv,ierr))
162: if (rank .eq. 0) then
163: write(*,'(A,I2/)') ' Number of converged approximate eigenpairs:',nconv
164: endif
166: ! ** Display eigenvalues and relative errors
167: if (nconv .gt. 0) then
168: if (rank .eq. 0) then
169: write(*,*) ' k ||T(k)x||'
170: write(*,*) '----------------- ------------------'
171: endif
172: do i=0,nconv-1
173: ! ** Get converged eigenpairs: (in this example they are always real)
174: PetscCallA(NEPGetEigenpair(nep,i,lambda,PETSC_NULL_SCALAR,x,PETSC_NULL_VEC,ierr))
176: ! ** Compute residual norm and error
177: PetscCallA(NEPComputeError(nep,i,NEP_ERROR_RELATIVE,norm,ierr))
178: if (rank .eq. 0) then
179: write(*,'(1P,E15.4,E18.4)') PetscRealPart(lambda), norm
180: endif
181: enddo
182: if (rank .eq. 0) then
183: write(*,*)
184: endif
185: endif
187: PetscCallA(NEPDestroy(nep,ierr))
188: PetscCallA(MatDestroy(F,ierr))
189: PetscCallA(MatDestroy(J,ierr))
190: PetscCallA(VecDestroy(x,ierr))
191: PetscCallA(SlepcFinalize(ierr))
192: end
194: ! --------------- Evaluate Function matrix T(lambda) ----------------
196: subroutine FormFunction(nep,lambda,fun,B,ctx,ierr)
197: use UserModule
198: implicit none
199: NEP nep
200: PetscScalar lambda, A(3), c, d
201: Mat fun,B
202: type(User) ctx
203: PetscReal h
204: PetscInt i, n, j(3), Istart, Iend, one, two, three
205: PetscErrorCode ierr
207: ! ** Compute Function entries and insert into matrix
208: PetscCall(MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr))
209: PetscCall(MatGetOwnershipRange(fun,Istart,Iend,ierr))
210: h = ctx%h
211: c = ctx%kappa/(lambda-ctx%kappa)
212: d = n
213: one = 1
214: two = 2
215: three = 3
217: ! ** Boundary points
218: if (Istart .eq. 0) then
219: i = 0
220: j(1) = 0
221: j(2) = 1
222: A(1) = 2.0*(d-lambda*h/3.0)
223: A(2) = -d-lambda*h/6.0
224: PetscCall(MatSetValues(fun,one,[i],two,j,A,INSERT_VALUES,ierr))
225: Istart = Istart + 1
226: endif
228: if (Iend .eq. n) then
229: i = n-1
230: j(1) = n-2
231: j(2) = n-1
232: A(1) = -d-lambda*h/6.0
233: A(2) = d-lambda*h/3.0+lambda*c
234: PetscCall(MatSetValues(fun,one,[i],two,j,A,INSERT_VALUES,ierr))
235: Iend = Iend - 1
236: endif
238: ! ** Interior grid points
239: do i=Istart,Iend-1
240: j(1) = i-1
241: j(2) = i
242: j(3) = i+1
243: A(1) = -d-lambda*h/6.0
244: A(2) = 2.0*(d-lambda*h/3.0)
245: A(3) = -d-lambda*h/6.0
246: PetscCall(MatSetValues(fun,one,[i],three,j,A,INSERT_VALUES,ierr))
247: enddo
249: ! ** Assemble matrix
250: PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr))
251: PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr))
253: end
255: ! --------------- Evaluate Jacobian matrix T'(lambda) ---------------
257: subroutine FormJacobian(nep,lambda,jac,ctx,ierr)
258: use UserModule
259: implicit none
260: NEP nep
261: PetscScalar lambda, A(3), c
262: Mat jac
263: type(User) ctx
264: PetscReal h
265: PetscInt i, n, j(3), Istart, Iend, one, two, three
266: PetscErrorCode ierr
268: ! ** Compute Jacobian entries and insert into matrix
269: PetscCall(MatGetSize(jac,n,PETSC_NULL_INTEGER,ierr))
270: PetscCall(MatGetOwnershipRange(jac,Istart,Iend,ierr))
271: h = ctx%h
272: c = ctx%kappa/(lambda-ctx%kappa)
273: one = 1
274: two = 2
275: three = 3
277: ! ** Boundary points
278: if (Istart .eq. 0) then
279: i = 0
280: j(1) = 0
281: j(2) = 1
282: A(1) = -2.0*h/3.0
283: A(2) = -h/6.0
284: PetscCall(MatSetValues(jac,one,[i],two,j,A,INSERT_VALUES,ierr))
285: Istart = Istart + 1
286: endif
288: if (Iend .eq. n) then
289: i = n-1
290: j(1) = n-2
291: j(2) = n-1
292: A(1) = -h/6.0
293: A(2) = -h/3.0-c*c
294: PetscCall(MatSetValues(jac,one,[i],two,j,A,INSERT_VALUES,ierr))
295: Iend = Iend - 1
296: endif
298: ! ** Interior grid points
299: do i=Istart,Iend-1
300: j(1) = i-1
301: j(2) = i
302: j(3) = i+1
303: A(1) = -h/6.0
304: A(2) = -2.0*h/3.0
305: A(3) = -h/6.0
306: PetscCall(MatSetValues(jac,one,[i],three,j,A,INSERT_VALUES,ierr))
307: enddo
309: ! ** Assemble matrix
310: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr))
311: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr))
313: end
315: !/*TEST
316: !
317: ! test:
318: ! suffix: 1
319: ! args: -nep_target 4
320: ! 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 = /"
321: ! requires: !single
322: !
323: !TEST*/