Actual source code: ex27f.F90

slepc-main 2024-11-15
Report Typos and Errors
  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*/