Actual source code: ex20f90.F90
slepc-3.17.1 2022-04-11
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*/