Actual source code: ex20f90.F90
slepc-3.11.2 2019-07-30
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*/