Actual source code: ex27f.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> ./ex27f [-help] [-n <n>] [all SLEPc options]
11: !
12: ! Description: Simple NLEIGS example. Fortran90 equivalent of ex27.c
13: !
14: ! The command line options are:
15: ! -n <n>, where <n> = matrix dimension
16: !
17: ! ----------------------------------------------------------------------
18: ! Solve T(lambda)x=0 using NLEIGS solver
19: ! with T(lambda) = -D+sqrt(lambda)*I
20: ! where D is the Laplacian operator in 1 dimension
21: ! and with the interpolation interval [.01,16]
22: ! ----------------------------------------------------------------------
23: !
24: #include <slepc/finclude/slepcnep.h>
26: module ex27fmodule
27: use slepcnep
28: implicit none
30: contains
31: ! --------------------------------------------------------------
32: ! FormFunction - Computes Function matrix T(lambda)
33: !
34: subroutine FormFunction(nep, lambda, fun, B, ctx, ierr)
35: use slepcnep
36: implicit none
38: NEP :: nep
39: PetscScalar :: lambda, val(0:2), t
40: Mat :: fun, B
41: PetscInt :: ctx, i, n, col(0:2), Istart, Iend, Istart0, Iend0
42: PetscBool :: FirstBlock = PETSC_FALSE, LastBlock = PETSC_FALSE
43: PetscErrorCode, intent(out) :: ierr
45: ! ** Compute Function entries and insert into matrix
46: t = sqrt(lambda)
47: PetscCall(MatGetSize(fun, n, PETSC_NULL_INTEGER, ierr))
48: PetscCall(MatGetOwnershipRange(fun, Istart, Iend, ierr))
49: if (Istart == 0) FirstBlock = PETSC_TRUE
50: if (Iend == n) LastBlock = PETSC_TRUE
51: val(0) = 1.0
52: val(1) = t - 2.0
53: val(2) = 1.0
55: Istart0 = Istart
56: if (FirstBlock) Istart0 = Istart + 1
57: Iend0 = Iend
58: if (LastBlock) Iend0 = Iend - 1
60: do i = Istart0, Iend0 - 1
61: col(0) = i - 1
62: col(1) = i
63: col(2) = i + 1
64: PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
65: end do
67: if (LastBlock) then
68: i = n - 1
69: col(0) = n - 2
70: col(1) = n - 1
71: val(0) = 1.0
72: val(1) = t - 2.0
73: PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
74: end if
76: if (FirstBlock) then
77: i = 0
78: col(0) = 0
79: col(1) = 1
80: val(0) = t - 2.0
81: val(1) = 1.0
82: PetscCall(MatSetValues(fun, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, val, INSERT_VALUES, ierr))
83: end if
85: ! ** Assemble matrix
86: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY, ierr))
87: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY, ierr))
88: PetscCall(MatAssemblyBegin(fun, MAT_FINAL_ASSEMBLY, ierr))
89: PetscCall(MatAssemblyEnd(fun, MAT_FINAL_ASSEMBLY, ierr))
91: end subroutine FormFunction
93: ! --------------------------------------------------------------
94: ! FormJacobian - Computes Jacobian matrix T'(lambda)
95: !
96: subroutine FormJacobian(nep, lambda, jac, ctx, ierr)
97: use slepcnep
98: implicit none
100: NEP :: nep
101: PetscScalar :: lambda, t
102: Mat :: jac
103: PetscInt :: ctx
104: Vec :: d
105: PetscErrorCode, intent(out) :: ierr
107: PetscCall(MatCreateVecs(jac, d, PETSC_NULL_VEC, ierr))
108: t = 0.5/sqrt(lambda)
109: PetscCall(VecSet(d, t, ierr))
110: PetscCall(MatDiagonalSet(jac, d, INSERT_VALUES, ierr))
111: PetscCall(VecDestroy(d, ierr))
113: end subroutine FormJacobian
115: ! --------------------------------------------------------------
116: ! ComputeSingularities - This is a user-defined routine to compute maxnp
117: ! points (at most) in the complex plane where the function T(.) is not analytic.
118: !
119: ! In this case, we discretize the singularity region (-inf,0)~(-1e+6,-1e-5)
120: !
121: ! Input Parameters:
122: ! nep - nonlinear eigensolver context
123: ! maxnp - on input number of requested points in the discretization (can be set)
124: ! xi - computed values of the discretization
125: ! dummy - optional user-defined monitor context (unused here)
126: !
127: subroutine ComputeSingularities(nep, maxnp, xi, dummy, ierr)
128: use slepcnep
129: implicit none
131: NEP :: nep
132: PetscInt :: maxnp, dummy
133: PetscScalar :: xi(0:maxnp - 1)
134: PetscReal :: h
135: PetscInt :: i
136: PetscErrorCode, intent(out) :: ierr
138: h = 11.0/real(maxnp - 1, PETSC_REAL_KIND)
139: xi(0) = -1e-5
140: xi(maxnp - 1) = -1e+6
141: do i = 1, maxnp - 2
142: xi(i) = -10**(-5 + h*i)
143: end do
144: ierr = 0
146: end subroutine ComputeSingularities
148: end module ex27fmodule
150: program ex27f
151: use slepcnep
152: use ex27fmodule
153: implicit none
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: ! Declarations
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: NEP :: nep
160: Mat :: A(2), F, J
161: NEPType :: ntype
162: PetscInt :: n = 100, nev, Istart, Iend, i, col
163: PetscErrorCode :: ierr
164: PetscBool :: terse, flg, split = PETSC_TRUE
165: PetscReal :: ia, ib, ic, id
166: RG :: rg
167: FN :: fn(2)
168: PetscScalar :: coeffs, sigma
169: character(len=128) :: string
170: PetscScalar, parameter :: one = 1.0
172: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: ! Beginning of program
174: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
177: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-n", n, flg, ierr))
178: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, "-split", split, flg, ierr))
179: if (split) then
180: write (string, *) 'Square root eigenproblem, n=', n, ' (split-form)\n'
181: else
182: write (string, *) 'Square root eigenproblem, n=', n, '\n'
183: end if
184: PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(string), ierr))
186: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: ! Create nonlinear eigensolver context and set options
188: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: PetscCallA(NEPCreate(PETSC_COMM_WORLD, nep, ierr))
191: PetscCallA(NEPSetType(nep, NEPNLEIGS, ierr))
192: PetscCallA(NEPNLEIGSSetSingularitiesFunction(nep, ComputeSingularities, 0, ierr))
193: PetscCallA(NEPGetRG(nep, rg, ierr))
194: PetscCallA(RGSetType(rg, RGINTERVAL, ierr))
195: ia = 0.01
196: ib = 16.0
197: #if defined(PETSC_USE_COMPLEX)
198: ic = -0.001
199: id = 0.001
200: #else
201: ic = 0.0
202: id = 0.0
203: #endif
204: PetscCallA(RGIntervalSetEndpoints(rg, ia, ib, ic, id, ierr))
205: sigma = 1.1
206: PetscCallA(NEPSetTarget(nep, sigma, ierr))
208: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209: ! Define the nonlinear problem
210: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212: if (split) then
213: ! ** Create matrices for the split form
214: PetscCallA(MatCreate(PETSC_COMM_WORLD, A(1), ierr))
215: PetscCallA(MatSetSizes(A(1), PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
216: PetscCallA(MatSetFromOptions(A(1), ierr))
217: PetscCallA(MatGetOwnershipRange(A(1), Istart, Iend, ierr))
218: coeffs = -2.0
219: do i = Istart, Iend - 1
220: if (i > 0) then
221: col = i - 1
222: PetscCallA(MatSetValue(A(1), i, col, one, INSERT_VALUES, ierr))
223: end if
224: if (i < n - 1) then
225: col = i + 1
226: PetscCallA(MatSetValue(A(1), i, col, one, INSERT_VALUES, ierr))
227: end if
228: PetscCallA(MatSetValue(A(1), i, i, coeffs, INSERT_VALUES, ierr))
229: end do
230: PetscCallA(MatAssemblyBegin(A(1), MAT_FINAL_ASSEMBLY, ierr))
231: PetscCallA(MatAssemblyEnd(A(1), MAT_FINAL_ASSEMBLY, ierr))
233: PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, one, A(2), ierr))
235: ! ** Define functions for the split form
236: PetscCallA(FNCreate(PETSC_COMM_WORLD, fn(1), ierr))
237: PetscCallA(FNSetType(fn(1), FNRATIONAL, ierr))
238: PetscCallA(FNRationalSetNumerator(fn(1), 1_PETSC_INT_KIND, [one], ierr))
239: PetscCallA(FNCreate(PETSC_COMM_WORLD, fn(2), ierr))
240: PetscCallA(FNSetType(fn(2), FNSQRT, ierr))
241: PetscCallA(NEPSetSplitOperator(nep, 2_PETSC_INT_KIND, A, fn, SUBSET_NONZERO_PATTERN, ierr))
242: else
243: ! ** Callback form: create matrix and set Function evaluation routine
244: PetscCallA(MatCreate(PETSC_COMM_WORLD, F, ierr))
245: PetscCallA(MatSetSizes(F, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
246: PetscCallA(MatSetFromOptions(F, ierr))
247: PetscCallA(MatSeqAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
248: PetscCallA(MatMPIAIJSetPreallocation(F, 3_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
249: PetscCallA(NEPSetFunction(nep, F, F, FormFunction, PETSC_NULL_INTEGER, ierr))
251: PetscCallA(MatCreate(PETSC_COMM_WORLD, J, ierr))
252: PetscCallA(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
253: PetscCallA(MatSetFromOptions(J, ierr))
254: PetscCallA(MatSeqAIJSetPreallocation(J, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
255: PetscCallA(MatMPIAIJSetPreallocation(J, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, ierr))
256: PetscCallA(NEPSetJacobian(nep, J, FormJacobian, PETSC_NULL_INTEGER, ierr))
257: end if
259: PetscCallA(NEPSetFromOptions(nep, ierr))
261: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
262: ! Solve the eigensystem
263: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
265: PetscCallA(NEPSolve(nep, ierr))
266: PetscCallA(NEPGetType(nep, ntype, ierr))
267: write (string, *) 'Solution method: ', ntype, '\n'
268: PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(string), ierr))
269: PetscCallA(NEPGetDimensions(nep, nev, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
270: write (string, *) 'Number of requested eigenvalues:', nev, '\n'
271: PetscCallA(PetscPrintf(PETSC_COMM_WORLD, trim(string), ierr))
273: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: ! Display solution and clean up
275: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: ! ** show detailed info unless -terse option is given by user
278: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-terse', terse, ierr))
279: if (terse) then
280: PetscCallA(NEPErrorView(nep, NEP_ERROR_BACKWARD, PETSC_NULL_VIEWER, ierr))
281: else
282: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL, ierr))
283: PetscCallA(NEPConvergedReasonView(nep, PETSC_VIEWER_STDOUT_WORLD, ierr))
284: PetscCallA(NEPErrorView(nep, NEP_ERROR_BACKWARD, PETSC_VIEWER_STDOUT_WORLD, ierr))
285: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD, ierr))
286: end if
288: if (split) then
289: PetscCallA(MatDestroy(A(1), ierr))
290: PetscCallA(MatDestroy(A(2), ierr))
291: PetscCallA(FNDestroy(fn(1), ierr))
292: PetscCallA(FNDestroy(fn(2), ierr))
293: else
294: PetscCallA(MatDestroy(F, ierr))
295: PetscCallA(MatDestroy(J, ierr))
296: end if
297: PetscCallA(NEPDestroy(nep, ierr))
298: PetscCallA(SlepcFinalize(ierr))
300: end program ex27f
302: !/*TEST
303: !
304: ! test:
305: ! suffix: 1
306: ! args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
307: ! requires: !single
308: ! filter: sed -e "s/[+-]0\.0*i//g"
309: !
310: ! test:
311: ! suffix: 2
312: ! args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
313: ! requires: !single
314: ! filter: sed -e "s/[+-]0\.0*i//g"
315: !
316: !TEST*/