Actual source code: ex27f.F90
slepc-3.22.1 2024-10-28
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: 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(MatGetOwnershipRange(A(1),Istart,Iend,ierr))
99: coeffs = -2.0
100: do i=Istart,Iend-1
101: if (i.gt.0) then
102: col = i-1
103: PetscCallA(MatSetValue(A(1),i,col,done,INSERT_VALUES,ierr))
104: end if
105: if (i.lt.n-1) then
106: col = i+1
107: PetscCallA(MatSetValue(A(1),i,col,done,INSERT_VALUES,ierr))
108: end if
109: PetscCallA(MatSetValue(A(1),i,i,coeffs,INSERT_VALUES,ierr))
110: end do
111: PetscCallA(MatAssemblyBegin(A(1),MAT_FINAL_ASSEMBLY,ierr))
112: PetscCallA(MatAssemblyEnd(A(1),MAT_FINAL_ASSEMBLY,ierr))
114: PetscCallA(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,done,A(2),ierr))
116: ! ** Define functions for the split form
117: PetscCallA(FNCreate(PETSC_COMM_WORLD,fn(1),ierr))
118: PetscCallA(FNSetType(fn(1),FNRATIONAL,ierr))
119: PetscCallA(FNRationalSetNumerator(fn(1),one,[done],ierr))
120: PetscCallA(FNCreate(PETSC_COMM_WORLD,fn(2),ierr))
121: PetscCallA(FNSetType(fn(2),FNSQRT,ierr))
122: PetscCallA(NEPSetSplitOperator(nep,two,A,fn,SUBSET_NONZERO_PATTERN,ierr))
123: else
124: ! ** Callback form: create matrix and set Function evaluation routine
125: PetscCallA(MatCreate(PETSC_COMM_WORLD,F,ierr))
126: PetscCallA(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
127: PetscCallA(MatSetFromOptions(F,ierr))
128: PetscCallA(MatSeqAIJSetPreallocation(F,three,PETSC_NULL_INTEGER_ARRAY,ierr))
129: PetscCallA(MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER_ARRAY,one,PETSC_NULL_INTEGER_ARRAY,ierr))
130: PetscCallA(NEPSetFunction(nep,F,F,FormFunction,PETSC_NULL_INTEGER,ierr))
132: PetscCallA(MatCreate(PETSC_COMM_WORLD,J,ierr))
133: PetscCallA(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr))
134: PetscCallA(MatSetFromOptions(J,ierr))
135: PetscCallA(MatSeqAIJSetPreallocation(J,one,PETSC_NULL_INTEGER_ARRAY,ierr))
136: PetscCallA(MatMPIAIJSetPreallocation(J,one,PETSC_NULL_INTEGER_ARRAY,one,PETSC_NULL_INTEGER_ARRAY,ierr))
137: PetscCallA(NEPSetJacobian(nep,J,FormJacobian,PETSC_NULL_INTEGER,ierr))
138: end if
140: PetscCallA(NEPSetFromOptions(nep,ierr))
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: ! Solve the eigensystem
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: PetscCallA(NEPSolve(nep,ierr))
147: PetscCallA(NEPGetType(nep,ntype,ierr))
148: write(string,*) 'Solution method: ',ntype,'\n'
149: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr))
150: PetscCallA(NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
151: write(string,*) 'Number of requested eigenvalues:',nev,'\n'
152: PetscCallA(PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr))
154: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: ! Display solution and clean up
156: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: ! ** show detailed info unless -terse option is given by user
159: PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr))
160: if (terse) then
161: PetscCallA(NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr))
162: else
163: PetscCallA(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr))
164: PetscCallA(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr))
165: PetscCallA(NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr))
166: PetscCallA(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr))
167: end if
169: if (split) then
170: PetscCallA(MatDestroy(A(1),ierr))
171: PetscCallA(MatDestroy(A(2),ierr))
172: PetscCallA(FNDestroy(fn(1),ierr))
173: PetscCallA(FNDestroy(fn(2),ierr))
174: else
175: PetscCallA(MatDestroy(F,ierr))
176: PetscCallA(MatDestroy(J,ierr))
177: end if
178: PetscCallA(NEPDestroy(nep,ierr))
179: PetscCallA(SlepcFinalize(ierr))
181: END PROGRAM main
183: ! --------------------------------------------------------------
184: !
185: ! FormFunction - Computes Function matrix T(lambda)
186: !
187: SUBROUTINE FormFunction(nep,lambda,fun,B,ctx,ierr)
188: #include <slepc/finclude/slepcnep.h>
189: use slepcnep
190: implicit none
192: NEP :: nep
193: PetscScalar :: lambda,val(0:2),t
194: Mat :: fun,B
195: PetscInt :: ctx,i,n,col(0:2),Istart,Iend,Istart0,Iend0,one,two,three
196: PetscErrorCode :: ierr
197: PetscBool :: FirstBlock=PETSC_FALSE, LastBlock=PETSC_FALSE
199: one = 1
200: two = 2
201: three = 3
203: ! ** Compute Function entries and insert into matrix
204: t = sqrt(lambda)
205: PetscCall(MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr))
206: PetscCall(MatGetOwnershipRange(fun,Istart,Iend,ierr))
207: if (Istart.eq.0) FirstBlock=PETSC_TRUE;
208: if (Iend.eq.n) LastBlock=PETSC_TRUE;
209: val(0)=1.0; val(1)=t-2.0; val(2)=1.0;
211: Istart0 = Istart
212: if (FirstBlock) Istart0 = Istart+1
213: Iend0 = Iend
214: if (LastBlock) Iend0 = Iend-1
216: do i=Istart0,Iend0-1
217: col(0) = i-1
218: col(1) = i
219: col(2) = i+1
220: PetscCall(MatSetValues(fun,one,[i],three,col,val,INSERT_VALUES,ierr))
221: end do
223: if (LastBlock) then
224: i = n-1
225: col(0) = n-2
226: col(1) = n-1
227: val(0) = 1.0
228: val(1) = t-2.0
229: PetscCall(MatSetValues(fun,one,[i],two,col,val,INSERT_VALUES,ierr))
230: end if
232: if (FirstBlock) then
233: i = 0
234: col(0) = 0
235: col(1) = 1
236: val(0) = t-2.0
237: val(1) = 1.0
238: PetscCall(MatSetValues(fun,one,[i],two,col,val,INSERT_VALUES,ierr))
239: end if
241: ! ** Assemble matrix
242: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
243: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
244: PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr))
245: PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr))
247: END SUBROUTINE FormFunction
249: ! --------------------------------------------------------------
250: !
251: ! FormJacobian - Computes Jacobian matrix T'(lambda)
252: !
253: SUBROUTINE FormJacobian(nep,lambda,jac,ctx,ierr)
254: #include <slepc/finclude/slepcnep.h>
255: USE slepcnep
256: implicit none
258: NEP :: nep
259: PetscScalar :: lambda,t
260: Mat :: jac
261: PetscInt :: ctx
262: PetscErrorCode :: ierr
263: Vec :: d
265: PetscCall(MatCreateVecs(jac,d,PETSC_NULL_VEC,ierr))
266: t = 0.5/sqrt(lambda)
267: PetscCall(VecSet(d,t,ierr))
268: PetscCall(MatDiagonalSet(jac,d,INSERT_VALUES,ierr))
269: PetscCall(VecDestroy(d,ierr))
271: END SUBROUTINE FormJacobian
273: ! --------------------------------------------------------------
274: !
275: ! ComputeSingularities - This is a user-defined routine to compute maxnp
276: ! points (at most) in the complex plane where the function T(.) is not analytic.
277: !
278: ! In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
279: !
280: ! Input Parameters:
281: ! nep - nonlinear eigensolver context
282: ! maxnp - on input number of requested points in the discretization (can be set)
283: ! xi - computed values of the discretization
284: ! dummy - optional user-defined monitor context (unused here)
285: !
286: SUBROUTINE ComputeSingularities(nep,maxnp,xi,dummy,ierr)
287: #include <slepc/finclude/slepcnep.h>
288: use slepcnep
289: implicit none
291: NEP :: nep
292: PetscInt :: maxnp, dummy
293: PetscScalar :: xi(0:maxnp-1)
294: PetscErrorCode :: ierr
295: PetscReal :: h
296: PetscInt :: i
298: h = 11.0/real(maxnp-1)
299: xi(0) = -1e-5
300: xi(maxnp-1) = -1e+6
301: do i=1,maxnp-2
302: xi(i) = -10**(-5+h*i)
303: end do
304: ierr = 0
306: END SUBROUTINE ComputeSingularities
308: !/*TEST
309: !
310: ! test:
311: ! suffix: 1
312: ! args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
313: ! requires: !single
314: ! filter: sed -e "s/[+-]0\.0*i//g"
315: !
316: ! test:
317: ! suffix: 2
318: ! args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
319: ! requires: !single
320: ! filter: sed -e "s/[+-]0\.0*i//g"
321: !
322: !TEST*/