Actual source code: ex27f90.F90
slepc-3.11.2 2019-07-30
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2019, 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: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
54: if (ierr .ne. 0) then
55: print*,'SlepcInitialize failed'
56: stop
57: end if
58: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-n",n,flg,ierr);CHKERRA(ierr)
59: call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-split",split,flg,ierr);CHKERRA(ierr)
60: if (split) then
61: write(string,*) 'Square root eigenproblem, n=',n,' (split-form)\n'
62: else
63: write(string,*) 'Square root eigenproblem, n=',n,'\n'
64: end if
65: call PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr);CHKERRA(ierr)
66: done = 1.0
67: one = 1
68: two = 2
69: three = 3
70:
71: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: ! Create nonlinear eigensolver context and set options
73: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: call NEPCreate(PETSC_COMM_WORLD,nep,ierr);CHKERRA(ierr)
76: call NEPSetType(nep,NEPNLEIGS,ierr);CHKERRA(ierr)
77: call NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,0,ierr);CHKERRA(ierr)
78: call NEPGetRG(nep,rg,ierr);CHKERRA(ierr)
79: call RGSetType(rg,RGINTERVAL,ierr);CHKERRA(ierr)
80: ia = 0.01
81: ib = 16.0
82: #if defined(PETSC_USE_COMPLEX)
83: ic = -0.001
84: id = 0.001
85: #else
86: ic = 0.0
87: id = 0.0
88: #endif
89: call RGIntervalSetEndpoints(rg,ia,ib,ic,id,ierr);CHKERRA(ierr)
90: sigma = 1.1
91: call NEPSetTarget(nep,sigma,ierr);CHKERRA(ierr)
93: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: ! Define the nonlinear problem
95: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: if (split) then
98: ! ** Create matrices for the split form
99: call MatCreate(PETSC_COMM_WORLD,A(1),ierr);CHKERRA(ierr)
100: call MatSetSizes(A(1),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
101: call MatSetFromOptions(A(1),ierr);CHKERRA(ierr)
102: call MatSetUp(A(1),ierr);CHKERRA(ierr)
103: call MatGetOwnershipRange(A(1),Istart,Iend,ierr);CHKERRA(ierr)
104: coeffs = -2.0
105: do i=Istart,Iend-1
106: if (i.gt.0) then
107: col = i-1
108: call MatSetValue(A(1),i,col,done,INSERT_VALUES,ierr);CHKERRA(ierr)
109: end if
110: if (i.lt.n-1) then
111: col = i+1
112: call MatSetValue(A(1),i,col,done,INSERT_VALUES,ierr);CHKERRA(ierr)
113: end if
114: call MatSetValue(A(1),i,i,coeffs,INSERT_VALUES,ierr);CHKERRA(ierr)
115: end do
116: call MatAssemblyBegin(A(1),MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
117: call MatAssemblyEnd(A(1),MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
119: call MatCreate(PETSC_COMM_WORLD,A(2),ierr);CHKERRA(ierr)
120: call MatSetSizes(A(2),PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
121: call MatSetFromOptions(A(2),ierr);CHKERRA(ierr)
122: call MatSetUp(A(2),ierr);CHKERRA(ierr)
123: call MatAssemblyBegin(A(2),MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
124: call MatAssemblyEnd(A(2),MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
125: call MatShift(A(2),done,ierr);CHKERRA(ierr)
127: ! ** Define funcions for the split form
128: call FNCreate(PETSC_COMM_WORLD,fn(1),ierr);CHKERRA(ierr)
129: call FNSetType(fn(1),FNRATIONAL,ierr);CHKERRA(ierr)
130: call FNRationalSetNumerator(fn(1),one,done,ierr);CHKERRA(ierr)
131: call FNCreate(PETSC_COMM_WORLD,fn(2),ierr);CHKERRA(ierr)
132: call FNSetType(fn(2),FNSQRT,ierr);CHKERRA(ierr)
133: call NEPSetSplitOperator(nep,two,A,fn,SUBSET_NONZERO_PATTERN,ierr);CHKERRA(ierr)
134: else
135: ! ** Callback form: create matrix and set Function evaluation routine
136: call MatCreate(PETSC_COMM_WORLD,F,ierr);CHKERRA(ierr)
137: call MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
138: call MatSetFromOptions(F,ierr);CHKERRA(ierr)
139: call MatSeqAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
140: call MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
141: Call MatSetUp(F,ierr);CHKERRA(ierr)
142: call NEPSetFunction(nep,F,F,FormFunction,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
144: call MatCreate(PETSC_COMM_WORLD,J,ierr);CHKERRA(ierr)
145: call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
146: call MatSetFromOptions(J,ierr);CHKERRA(ierr)
147: call MatSeqAIJSetPreallocation(J,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
148: call MatMPIAIJSetPreallocation(J,one,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
149: call MatSetUp(J,ierr);CHKERRA(ierr)
150: call NEPSetJacobian(nep,J,FormJacobian,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
151: end if
153: call NEPSetFromOptions(nep,ierr);CHKERRA(ierr)
155: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: ! Solve the eigensystem
157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: call NEPSolve(nep,ierr);CHKERRA(ierr)
160: call NEPGetType(nep,ntype,ierr);CHKERRA(ierr)
161: write(string,*) 'Solution method: ',ntype,'\n'
162: call PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr);CHKERRA(ierr)
163: call NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
164: write(string,*) 'Number of requested eigenvalues:',nev,'\n'
165: call PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr);CHKERRA(ierr)
167: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168: ! Display solution and clean up
169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: ! ** show detailed info unless -terse option is given by user
172: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
173: if (terse) then
174: call NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
175: else
176: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
177: call NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
178: call NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
179: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
180: end if
181:
182: if (split) then
183: call MatDestroy(A(1),ierr);CHKERRA(ierr)
184: call MatDestroy(A(2),ierr);CHKERRA(ierr)
185: call FNDestroy(fn(1),ierr);CHKERRA(ierr)
186: call FNDestroy(fn(2),ierr);CHKERRA(ierr)
187: else
188: call MatDestroy(F,ierr);CHKERRA(ierr)
189: call MatDestroy(J,ierr);CHKERRA(ierr)
190: end if
191: call NEPDestroy(nep,ierr)
192: call SlepcFinalize(ierr)
194: END PROGRAM main
196: ! --------------------------------------------------------------
197: !
198: ! FormFunction - Computes Function matrix T(lambda)
199: !
200: SUBROUTINE FormFunction(nep,lambda,fun,B,ctx,ierr)
201: #include <slepc/finclude/slepcnep.h>
202: use slepcnep
203: implicit none
205: NEP :: nep
206: PetscScalar :: lambda,val(0:2),t
207: Mat :: fun,B
208: PetscInt :: ctx,i,n,col(0:2),Istart,Iend,Istart0,Iend0,one,two,three
209: PetscErrorCode :: ierr
210: PetscBool :: FirstBlock=PETSC_FALSE, LastBlock=PETSC_FALSE
212: one = 1
213: two = 2
214: three = 3
216: ! ** Compute Function entries and insert into matrix
217: t = sqrt(lambda)
218: call MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
219: call MatGetOwnershipRange(fun,Istart,Iend,ierr);CHKERRA(ierr)
220: if (Istart.eq.0) FirstBlock=PETSC_TRUE;
221: if (Iend.eq.n) LastBlock=PETSC_TRUE;
222: val(0)=1.0; val(1)=t-2.0; val(2)=1.0;
224: Istart0 = Istart
225: if (FirstBlock) Istart0 = Istart+1
226: Iend0 = Iend
227: if (LastBlock) Iend0 = Iend-1
229: do i=Istart0,Iend0-1
230: col(0) = i-1
231: col(1) = i
232: col(2) = i+1
233: call MatSetValues(fun,one,i,three,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
234: end do
236: if (LastBlock) then
237: i = n-1
238: col(0) = n-2
239: col(1) = n-1
240: val(0) = 1.0
241: val(1) = t-2.0
242: call MatSetValues(fun,one,i,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
243: end if
245: if (FirstBlock) then
246: i = 0
247: col(0) = 0
248: col(1) = 1
249: val(0) = t-2.0
250: val(1) = 1.0
251: call MatSetValues(fun,one,i,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
252: end if
254: ! ** Assemble matrix
255: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
256: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
257: call MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
258: call MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr)
260: END SUBROUTINE FormFunction
262: ! --------------------------------------------------------------
263: !
264: ! FormJacobian - Computes Jacobian matrix T'(lambda)
265: !
266: SUBROUTINE FormJacobian(nep,lambda,jac,ctx,ierr)
267: #include <slepc/finclude/slepcnep.h>
268: USE slepcnep
269: implicit none
271: NEP :: nep
272: PetscScalar :: lambda,t
273: Mat :: jac
274: PetscInt :: ctx
275: PetscErrorCode :: ierr
276: Vec :: d
278: call MatCreateVecs(jac,d,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
279: t = 0.5/sqrt(lambda)
280: call VecSet(d,t,ierr);CHKERRA(ierr)
281: call MatDiagonalSet(jac,d,INSERT_VALUES,ierr);CHKERRA(ierr)
282: calL VecDestroy(d,ierr);CHKERRA(ierr)
284: END SUBROUTINE FormJacobian
286: ! --------------------------------------------------------------
287: !
288: ! ComputeSingularities - This is a user-defined routine to compute maxnp
289: ! points (at most) in the complex plane where the function T(.) is not analytic.
290: !
291: ! In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
292: !
293: ! Input Parameters:
294: ! nep - nonlinear eigensolver context
295: ! maxnp - on input number of requested points in the discretization (can be set)
296: ! xi - computed values of the discretization
297: ! dummy - optional user-defined monitor context (unused here)
298: !
299: SUBROUTINE ComputeSingularities(nep,maxnp,xi,dummy,ierr)
300: #include <slepc/finclude/slepcnep.h>
301: use slepcnep
302: implicit none
304: NEP :: nep
305: PetscInt :: maxnp, dummy
306: PetscScalar :: xi(0:maxnp-1)
307: PetscErrorCode :: ierr
308: PetscReal :: h
309: PetscInt :: i
311: h = 11.0/real(maxnp-1)
312: xi(0) = -1e-5
313: xi(maxnp-1) = -1e+6
314: do i=1,maxnp-2
315: xi(i) = -10**(-5+h*i)
316: end do
318: END SUBROUTINE ComputeSingularities
320: !/*TEST
321: !
322: ! test:
323: ! suffix: 1
324: ! args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
325: ! requires: !single
326: !
327: ! test:
328: ! suffix: 2
329: ! args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
330: ! requires: !single
331: !
332: !TEST*/