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: 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
210: val(1)=t-2.0
211: val(2)=1.0
213: Istart0 = Istart
214: if (FirstBlock) Istart0 = Istart+1
215: Iend0 = Iend
216: if (LastBlock) Iend0 = Iend-1
218: do i=Istart0,Iend0-1
219: col(0) = i-1
220: col(1) = i
221: col(2) = i+1
222: PetscCall(MatSetValues(fun,one,[i],three,col,val,INSERT_VALUES,ierr))
223: end do
225: if (LastBlock) then
226: i = n-1
227: col(0) = n-2
228: col(1) = n-1
229: val(0) = 1.0
230: val(1) = t-2.0
231: PetscCall(MatSetValues(fun,one,[i],two,col,val,INSERT_VALUES,ierr))
232: end if
234: if (FirstBlock) then
235: i = 0
236: col(0) = 0
237: col(1) = 1
238: val(0) = t-2.0
239: val(1) = 1.0
240: PetscCall(MatSetValues(fun,one,[i],two,col,val,INSERT_VALUES,ierr))
241: end if
243: ! ** Assemble matrix
244: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr))
245: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr))
246: PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr))
247: PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr))
249: END SUBROUTINE FormFunction
251: ! --------------------------------------------------------------
252: !
253: ! FormJacobian - Computes Jacobian matrix T'(lambda)
254: !
255: SUBROUTINE FormJacobian(nep,lambda,jac,ctx,ierr)
256: #include <slepc/finclude/slepcnep.h>
257: USE slepcnep
258: implicit none
260: NEP :: nep
261: PetscScalar :: lambda,t
262: Mat :: jac
263: PetscInt :: ctx
264: PetscErrorCode :: ierr
265: Vec :: d
267: PetscCall(MatCreateVecs(jac,d,PETSC_NULL_VEC,ierr))
268: t = 0.5/sqrt(lambda)
269: PetscCall(VecSet(d,t,ierr))
270: PetscCall(MatDiagonalSet(jac,d,INSERT_VALUES,ierr))
271: PetscCall(VecDestroy(d,ierr))
273: END SUBROUTINE FormJacobian
275: ! --------------------------------------------------------------
276: !
277: ! ComputeSingularities - This is a user-defined routine to compute maxnp
278: ! points (at most) in the complex plane where the function T(.) is not analytic.
279: !
280: ! In this case, we discretize the singularity region (-inf,0)~(-1e+6,-1e-5)
281: !
282: ! Input Parameters:
283: ! nep - nonlinear eigensolver context
284: ! maxnp - on input number of requested points in the discretization (can be set)
285: ! xi - computed values of the discretization
286: ! dummy - optional user-defined monitor context (unused here)
287: !
288: SUBROUTINE ComputeSingularities(nep,maxnp,xi,dummy,ierr)
289: #include <slepc/finclude/slepcnep.h>
290: use slepcnep
291: implicit none
293: NEP :: nep
294: PetscInt :: maxnp, dummy
295: PetscScalar :: xi(0:maxnp-1)
296: PetscErrorCode :: ierr
297: PetscReal :: h
298: PetscInt :: i
300: h = 11.0/real(maxnp-1)
301: xi(0) = -1e-5
302: xi(maxnp-1) = -1e+6
303: do i=1,maxnp-2
304: xi(i) = -10**(-5+h*i)
305: end do
306: ierr = 0
308: END SUBROUTINE ComputeSingularities
310: !/*TEST
311: !
312: ! test:
313: ! suffix: 1
314: ! args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
315: ! requires: !single
316: ! filter: sed -e "s/[+-]0\.0*i//g"
317: !
318: ! test:
319: ! suffix: 2
320: ! args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
321: ! requires: !single
322: ! filter: sed -e "s/[+-]0\.0*i//g"
323: !
324: !TEST*/