Actual source code: ex27f90.F90
slepc-3.18.2 2023-01-26
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> ./ex27f90 [-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: PROGRAM main
25: #include <slepc/finclude/slepcnep.h>
26: USE slepcnep
27: implicit none
29: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: ! Declarations
31: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
33: NEP :: nep
34: Mat :: A(2),F,J
35: NEPType :: ntype
36: PetscInt :: n=100,nev,Istart,Iend,i,col,one,two,three
37: PetscErrorCode :: ierr
38: PetscBool :: terse,flg,split=PETSC_TRUE
39: PetscReal :: ia,ib,ic,id
40: RG :: rg
41: FN :: fn(2)
42: PetscScalar :: coeffs,sigma,done
43: CHARACTER(LEN=128) :: string
45: ! NOTE: Any user-defined Fortran routines (such as ComputeSingularities)
46: ! MUST be declared as external.
47: external ComputeSingularities, FormFunction, FormJacobian
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: ! Beginning of program
51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
54: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-n",n,flg,ierr))
55: PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-split",split,flg,ierr))
56: if (split) then
57: write(string,*) 'Square root eigenproblem, n=',n,' (split-form)\n'
58: else
59: write(string,*) 'Square root eigenproblem, n=',n,'\n'
60: end if
61: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr))
62: done = 1.0
63: one = 1
64: two = 2
65: three = 3
67: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: ! Create nonlinear eigensolver context and set options
69: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: PetscCallA(NEPCreate(PETSC_COMM_WORLD,nep,ierr))
72: PetscCallA(NEPSetType(nep,NEPNLEIGS,ierr))
73: PetscCallA(NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,0,ierr))
74: PetscCallA(NEPGetRG(nep,rg,ierr))
75: PetscCallA(RGSetType(rg,RGINTERVAL,ierr))
76: ia = 0.01
77: ib = 16.0
78: #if defined(PETSC_USE_COMPLEX)
79: ic = -0.001
80: id = 0.001
81: #else
82: ic = 0.0
83: id = 0.0
84: #endif
85: PetscCallA(RGIntervalSetEndpoints(rg,ia,ib,ic,id,ierr))
86: sigma = 1.1
87: PetscCallA(NEPSetTarget(nep,sigma,ierr))
89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: ! Define the nonlinear problem
91: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: if (split) then
94: ! ** Create matrices for the split form
95: PetscCallA(MatCreate(PETSC_COMM_WORLD,A(1),ierr))
96: PetscCallA(MatSetSizes(A(1),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
97: PetscCallA(MatSetFromOptions(A(1),ierr))
98: PetscCallA(MatSetUp(A(1),ierr))
99: PetscCallA(MatGetOwnershipRange(A(1),Istart,Iend,ierr))
100: coeffs = -2.0
101: do i=Istart,Iend-1
102: if (i.gt.0) then
103: col = i-1
104: PetscCallA(MatSetValue(A(1),i,col,done,INSERT_VALUES,ierr))
105: end if
106: if (i.lt.n-1) then
107: col = i+1
108: PetscCallA(MatSetValue(A(1),i,col,done,INSERT_VALUES,ierr))
109: end if
110: PetscCallA(MatSetValue(A(1),i,i,coeffs,INSERT_VALUES,ierr))
111: end do
112: PetscCallA(MatAssemblyBegin(A(1),MAT_FINAL_ASSEMBLY,ierr))
113: PetscCallA(MatAssemblyEnd(A(1),MAT_FINAL_ASSEMBLY,ierr))
115: PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,done,A(2),ierr))
117: ! ** Define functions for the split form
118: PetscCallA(FNCreate(PETSC_COMM_WORLD,fn(1),ierr))
119: PetscCallA(FNSetType(fn(1),FNRATIONAL,ierr))
120: PetscCallA(FNRationalSetNumerator(fn(1),one,done,ierr))
121: PetscCallA(FNCreate(PETSC_COMM_WORLD,fn(2),ierr))
122: PetscCallA(FNSetType(fn(2),FNSQRT,ierr))
123: PetscCallA(NEPSetSplitOperator(nep,two,A,fn,SUBSET_NONZERO_PATTERN,ierr))
124: else
125: ! ** Callback form: create matrix and set Function evaluation routine
126: PetscCallA(MatCreate(PETSC_COMM_WORLD,F,ierr))
127: PetscCallA(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
128: PetscCallA(MatSetFromOptions(F,ierr))
129: PetscCallA(MatSeqAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,ierr))
130: PetscCallA(MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr))
131: PetscCallA(MatSetUp(F,ierr))
132: PetscCallA(NEPSetFunction(nep,F,F,FormFunction,PETSC_NULL_INTEGER,ierr))
134: PetscCallA(MatCreate(PETSC_COMM_WORLD,J,ierr))
135: PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
136: PetscCallA(MatSetFromOptions(J,ierr))
137: PetscCallA(MatSeqAIJSetPreallocation(J,one,PETSC_NULL_INTEGER,ierr))
138: PetscCallA(MatMPIAIJSetPreallocation(J,one,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr))
139: PetscCallA(MatSetUp(J,ierr))
140: PetscCallA(NEPSetJacobian(nep,J,FormJacobian,PETSC_NULL_INTEGER,ierr))
141: end if
143: PetscCallA(NEPSetFromOptions(nep,ierr))
145: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: ! Solve the eigensystem
147: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: PetscCallA(NEPSolve(nep,ierr))
150: PetscCallA(NEPGetType(nep,ntype,ierr))
151: write(string,*) 'Solution method: ',ntype,'\n'
152: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr))
153: PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
154: write(string,*) 'Number of requested eigenvalues:',nev,'\n'
155: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr))
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: ! Display solution and clean up
159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: ! ** show detailed info unless -terse option is given by user
162: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
163: if (terse) then
164: PetscCallA(NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr))
165: else
166: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
167: PetscCallA(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr))
168: PetscCallA(NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr))
169: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
170: end if
172: if (split) then
173: PetscCallA(MatDestroy(A(1),ierr))
174: PetscCallA(MatDestroy(A(2),ierr))
175: PetscCallA(FNDestroy(fn(1),ierr))
176: PetscCallA(FNDestroy(fn(2),ierr))
177: else
178: PetscCallA(MatDestroy(F,ierr))
179: PetscCallA(MatDestroy(J,ierr))
180: end if
181: PetscCallA(NEPDestroy(nep,ierr))
182: PetscCallA(SlepcFinalize(ierr))
184: END PROGRAM main
186: ! --------------------------------------------------------------
187: !
188: ! FormFunction - Computes Function matrix T(lambda)
189: !
190: SUBROUTINE FormFunction(nep,lambda,fun,B,ctx,ierr)
191: #include <slepc/finclude/slepcnep.h>
192: use slepcnep
193: implicit none
195: NEP :: nep
196: PetscScalar :: lambda,val(0:2),t
197: Mat :: fun,B
198: PetscInt :: ctx,i,n,col(0:2),Istart,Iend,Istart0,Iend0,one,two,three
199: PetscErrorCode :: ierr
200: PetscBool :: FirstBlock=PETSC_FALSE, LastBlock=PETSC_FALSE
202: one = 1
203: two = 2
204: three = 3
206: ! ** Compute Function entries and insert into matrix
207: t = sqrt(lambda)
208: PetscCall(MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr))
209: PetscCall(MatGetOwnershipRange(fun,Istart,Iend,ierr))
210: if (Istart.eq.0) FirstBlock=PETSC_TRUE;
211: if (Iend.eq.n) LastBlock=PETSC_TRUE;
212: val(0)=1.0; val(1)=t-2.0; val(2)=1.0;
214: Istart0 = Istart
215: if (FirstBlock) Istart0 = Istart+1
216: Iend0 = Iend
217: if (LastBlock) Iend0 = Iend-1
219: do i=Istart0,Iend0-1
220: col(0) = i-1
221: col(1) = i
222: col(2) = i+1
223: PetscCall(MatSetValues(fun,one,i,three,col,val,INSERT_VALUES,ierr))
224: end do
226: if (LastBlock) then
227: i = n-1
228: col(0) = n-2
229: col(1) = n-1
230: val(0) = 1.0
231: val(1) = t-2.0
232: PetscCall(MatSetValues(fun,one,i,two,col,val,INSERT_VALUES,ierr))
233: end if
235: if (FirstBlock) then
236: i = 0
237: col(0) = 0
238: col(1) = 1
239: val(0) = t-2.0
240: val(1) = 1.0
241: PetscCall(MatSetValues(fun,one,i,two,col,val,INSERT_VALUES,ierr))
242: end if
244: ! ** Assemble matrix
245: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
246: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
247: PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr))
248: PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr))
250: END SUBROUTINE FormFunction
252: ! --------------------------------------------------------------
253: !
254: ! FormJacobian - Computes Jacobian matrix T'(lambda)
255: !
256: SUBROUTINE FormJacobian(nep,lambda,jac,ctx,ierr)
257: #include <slepc/finclude/slepcnep.h>
258: USE slepcnep
259: implicit none
261: NEP :: nep
262: PetscScalar :: lambda,t
263: Mat :: jac
264: PetscInt :: ctx
265: PetscErrorCode :: ierr
266: Vec :: d
268: PetscCall(MatCreateVecs(jac,d,PETSC_NULL_VEC,ierr))
269: t = 0.5/sqrt(lambda)
270: PetscCall(VecSet(d,t,ierr))
271: PetscCall(MatDiagonalSet(jac,d,INSERT_VALUES,ierr))
272: PetscCall(VecDestroy(d,ierr))
274: END SUBROUTINE FormJacobian
276: ! --------------------------------------------------------------
277: !
278: ! ComputeSingularities - This is a user-defined routine to compute maxnp
279: ! points (at most) in the complex plane where the function T(.) is not analytic.
280: !
281: ! In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
282: !
283: ! Input Parameters:
284: ! nep - nonlinear eigensolver context
285: ! maxnp - on input number of requested points in the discretization (can be set)
286: ! xi - computed values of the discretization
287: ! dummy - optional user-defined monitor context (unused here)
288: !
289: SUBROUTINE ComputeSingularities(nep,maxnp,xi,dummy,ierr)
290: #include <slepc/finclude/slepcnep.h>
291: use slepcnep
292: implicit none
294: NEP :: nep
295: PetscInt :: maxnp, dummy
296: PetscScalar :: xi(0:maxnp-1)
297: PetscErrorCode :: ierr
298: PetscReal :: h
299: PetscInt :: i
301: h = 11.0/real(maxnp-1)
302: xi(0) = -1e-5
303: xi(maxnp-1) = -1e+6
304: do i=1,maxnp-2
305: xi(i) = -10**(-5+h*i)
306: end do
307: ierr = 0
309: END SUBROUTINE ComputeSingularities
311: !/*TEST
312: !
313: ! test:
314: ! suffix: 1
315: ! args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
316: ! requires: !single
317: ! filter: sed -e "s/[+-]0\.0*i//g"
318: !
319: ! test:
320: ! suffix: 2
321: ! args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
322: ! requires: !single
323: ! filter: sed -e "s/[+-]0\.0*i//g"
324: !
325: !TEST*/