Actual source code: ex20f.F90
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: #include <slepc/finclude/slepcnep.h>
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: ! User-defined module with application context and callback functions
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: module ex20fmodule
30: use slepcnep
31: type User
32: PetscScalar :: kappa
33: PetscReal :: h
34: end type User
36: contains
37: ! --------------- Evaluate Function matrix T(lambda) ----------------
39: subroutine FormFunction(nep, lambda, fun, B, ctx, ierr)
40: implicit none
41: NEP :: nep
42: PetscScalar :: lambda, A(3), c, d
43: Mat :: fun, B
44: type(User) :: ctx
45: PetscReal :: h
46: PetscInt :: i, n, j(3), Istart, Iend
47: PetscErrorCode, intent(out) :: ierr
49: ! ** Compute Function entries and insert into matrix
50: PetscCall(MatGetSize(fun, n, PETSC_NULL_INTEGER, ierr))
51: PetscCall(MatGetOwnershipRange(fun, Istart, Iend, ierr))
52: h = ctx%h
53: c = ctx%kappa/(lambda - ctx%kappa)
54: d = n
56: ! ** Boundary points
57: if (Istart == 0) then
58: i = 0
59: j(1) = 0
60: j(2) = 1
61: A(1) = 2.0*(d - lambda*h/3.0)
62: A(2) = -d - lambda*h/6.0
63: PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
64: Istart = Istart + 1
65: end if
67: if (Iend == n) then
68: i = n - 1
69: j(1) = n - 2
70: j(2) = n - 1
71: A(1) = -d - lambda*h/6.0
72: A(2) = d - lambda*h/3.0 + lambda*c
73: PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
74: Iend = Iend - 1
75: end if
77: ! ** Interior grid points
78: do i = Istart, Iend - 1
79: j(1) = i - 1
80: j(2) = i
81: j(3) = i + 1
82: A(1) = -d - lambda*h/6.0
83: A(2) = 2.0*(d - lambda*h/3.0)
84: A(3) = -d - lambda*h/6.0
85: PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
86: end do
88: ! ** Assemble matrix
89: PetscCall(MatAssemblyBegin(fun, MAT_FINAL_ASSEMBLY, ierr))
90: PetscCall(MatAssemblyEnd(fun, MAT_FINAL_ASSEMBLY, ierr))
92: end subroutine
94: ! --------------- Evaluate Jacobian matrix T'(lambda) ---------------
96: subroutine FormJacobian(nep, lambda, jac, ctx, ierr)
97: implicit none
98: NEP :: nep
99: PetscScalar :: lambda, A(3), c
100: Mat :: jac
101: type(User) :: ctx
102: PetscReal :: h
103: PetscInt :: i, n, j(3), Istart, Iend
104: PetscErrorCode, intent(out) :: ierr
106: ! ** Compute Jacobian entries and insert into matrix
107: PetscCall(MatGetSize(jac, n, PETSC_NULL_INTEGER, ierr))
108: PetscCall(MatGetOwnershipRange(jac, Istart, Iend, ierr))
109: h = ctx%h
110: c = ctx%kappa/(lambda - ctx%kappa)
112: ! ** Boundary points
113: if (Istart == 0) then
114: i = 0
115: j(1) = 0
116: j(2) = 1
117: A(1) = -2.0*h/3.0
118: A(2) = -h/6.0
119: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
120: Istart = Istart + 1
121: end if
123: if (Iend == n) then
124: i = n - 1
125: j(1) = n - 2
126: j(2) = n - 1
127: A(1) = -h/6.0
128: A(2) = -h/3.0 - c*c
129: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
130: Iend = Iend - 1
131: end if
133: ! ** Interior grid points
134: do i = Istart, Iend - 1
135: j(1) = i - 1
136: j(2) = i
137: j(3) = i + 1
138: A(1) = -h/6.0
139: A(2) = -2.0*h/3.0
140: A(3) = -h/6.0
141: PetscCall(MatSetValues(jac, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, j, A, INSERT_VALUES, ierr))
142: end do
144: ! ** Assemble matrix
145: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY, ierr))
146: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY, ierr))
148: end subroutine
150: end module ex20fmodule
152: program ex20f
153: use ex20fmodule
154: implicit none
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: ! Declarations
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: NEP :: nep ! nonlinear eigensolver
161: Vec :: x, v ! eigenvector, auxiliary vector
162: PetscScalar :: lambda ! eigenvalue
163: Mat :: F, J ! Function and Jacobian matrices
164: type(User) :: ctx ! user-defined context
165: NEPType :: tname
166: PetscInt :: n, i, nev, its, maxit, nconv
167: PetscReal :: tol, norm
168: PetscMPIInt :: rank
169: PetscBool :: flg
170: PetscErrorCode :: ierr
171: PetscScalar, parameter :: one = 1.0
173: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: ! Beginning of program
175: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
178: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
179: n = 128
180: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
181: if (rank == 0) then
182: write (*, '(/a,i4)') 'Nonlinear Eigenproblem, n =', n
183: end if
185: ctx%h = 1.0/real(n, PETSC_REAL_KIND)
186: ctx%kappa = 1.0
188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: ! Create matrix data structure to hold the Function and the Jacobian
190: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: PetscCallA(MatCreate(PETSC_COMM_WORLD, F, ierr))
193: PetscCallA(MatSetSizes(F, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
194: PetscCallA(MatSetFromOptions(F, ierr))
195: PetscCallA(MatSeqAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
196: PetscCallA(MatMPIAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
198: PetscCallA(MatCreate(PETSC_COMM_WORLD, J, ierr))
199: PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
200: PetscCallA(MatSetFromOptions(J, ierr))
201: PetscCallA(MatSeqAIJSetPreallocation(J, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
202: PetscCallA(MatMPIAIJSetPreallocation(J, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
204: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: ! Create the eigensolver and set various options
206: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: ! ** Create eigensolver context
209: PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))
211: ! ** Set routines for evaluation of Function and Jacobian
212: PetscCallA(NEPSetFunction(nep, F, F, FormFunction, ctx, ierr))
213: PetscCallA(NEPSetJacobian(nep, J, FormJacobian, ctx, ierr))
215: ! ** Customize nonlinear solver
216: tol = 1e-9
217: PetscCallA(NEPSetTolerances(nep, tol, PETSC_CURRENT_INTEGER, ierr))
218: PetscCallA(NEPSetDimensions(nep, 1_PETSC_INT_KIND, PETSC_DETERMINE_INTEGER, PETSC_DETERMINE_INTEGER, ierr))
220: ! ** Set solver parameters at runtime
221: PetscCallA(NEPSetFromOptions(nep, ierr))
223: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
224: ! Solve the eigensystem
225: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227: ! ** Evaluate initial guess
228: PetscCallA(MatCreateVecs(F, x, PETSC_NULL_VEC, ierr))
229: PetscCallA(VecDuplicate(x, v, ierr))
230: PetscCallA(VecSet(v, one, ierr))
231: PetscCallA(NEPSetInitialSpace(nep, 1_PETSC_INT_KIND, [v], ierr))
232: PetscCallA(VecDestroy(v, ierr))
234: ! ** Call the solver
235: PetscCallA(NEPSolve(nep, ierr))
236: PetscCallA(NEPGetIterationNumber(nep, its, ierr))
237: if (rank == 0) then
238: write (*, '(a,i3)') ' Number of NEP iterations =', its
239: end if
241: ! ** Optional: Get some information from the solver and display it
242: PetscCallA(NEPGetType(nep, tname, ierr))
243: if (rank == 0) then
244: write (*, '(a,a10)') ' Solution method: ', tname
245: end if
246: PetscCallA(NEPGetDimensions(nep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
247: if (rank == 0) then
248: write (*, '(a,i4)') ' Number of requested eigenvalues:', nev
249: end if
250: PetscCallA(NEPGetTolerances(nep, tol, maxit, ierr))
251: if (rank == 0) then
252: write (*, '(a,f12.9,a,i5)') ' Stopping condition: tol=', tol, ', maxit=', maxit
253: end if
255: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
256: ! Display solution and clean up
257: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259: PetscCallA(NEPGetConverged(nep, nconv, ierr))
260: if (rank == 0) then
261: write (*, '(a,i2/)') ' Number of converged approximate eigenpairs:', nconv
262: end if
264: ! ** Display eigenvalues and relative errors
265: if (nconv > 0) then
266: if (rank == 0) then
267: write (*, *) ' k ||T(k)x||'
268: write (*, *) '----------------- ------------------'
269: end if
270: do i = 0, nconv - 1
271: ! ** Get converged eigenpairs: (in this example they are always real)
272: PetscCallA(NEPGetEigenpair(nep, i, lambda, PETSC_NULL_SCALAR, x, PETSC_NULL_VEC, ierr))
274: ! ** Compute residual norm and error
275: PetscCallA(NEPComputeError(nep, i, NEP_ERROR_RELATIVE, norm, ierr))
276: if (rank == 0) then
277: write (*, '(1p,e15.4,e18.4)') PetscRealPart(lambda), norm
278: end if
279: end do
280: if (rank == 0) then
281: write (*, *)
282: end if
283: end if
285: PetscCallA(NEPDestroy(nep, ierr))
286: PetscCallA(MatDestroy(F, ierr))
287: PetscCallA(MatDestroy(J, ierr))
288: PetscCallA(VecDestroy(x, ierr))
289: PetscCallA(SlepcFinalize(ierr))
290: end program ex20f
292: !/*TEST
293: !
294: ! test:
295: ! suffix: 1
296: ! args: -nep_target 4
297: ! 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 = /"
298: ! requires: !single
299: !
300: !TEST*/