Actual source code: ex27f90.F90

slepc-3.17.0 2022-03-31
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> ./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

 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 MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,done,A(2),ierr);CHKERRA(ierr)

121:      ! ** Define functions for the split form
122:      call FNCreate(PETSC_COMM_WORLD,fn(1),ierr);CHKERRA(ierr)
123:      call FNSetType(fn(1),FNRATIONAL,ierr);CHKERRA(ierr)
124:      call FNRationalSetNumerator(fn(1),one,done,ierr);CHKERRA(ierr)
125:      call FNCreate(PETSC_COMM_WORLD,fn(2),ierr);CHKERRA(ierr)
126:      call FNSetType(fn(2),FNSQRT,ierr);CHKERRA(ierr)
127:      call NEPSetSplitOperator(nep,two,A,fn,SUBSET_NONZERO_PATTERN,ierr);CHKERRA(ierr)
128:   else
129:     ! ** Callback form: create matrix and set Function evaluation routine
130:     call MatCreate(PETSC_COMM_WORLD,F,ierr);CHKERRA(ierr)
131:     call MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
132:     call MatSetFromOptions(F,ierr);CHKERRA(ierr)
133:     call MatSeqAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
134:     call MatMPIAIJSetPreallocation(F,three,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
135:     Call MatSetUp(F,ierr);CHKERRA(ierr)
136:     call NEPSetFunction(nep,F,F,FormFunction,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)

138:     call MatCreate(PETSC_COMM_WORLD,J,ierr);CHKERRA(ierr)
139:     call MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
140:     call MatSetFromOptions(J,ierr);CHKERRA(ierr)
141:     call MatSeqAIJSetPreallocation(J,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
142:     call MatMPIAIJSetPreallocation(J,one,PETSC_NULL_INTEGER,one,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
143:     call MatSetUp(J,ierr);CHKERRA(ierr)
144:     call NEPSetJacobian(nep,J,FormJacobian,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
145:   end if

147:   call NEPSetFromOptions(nep,ierr);CHKERRA(ierr)

149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !     Solve the eigensystem
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

153:   call NEPSolve(nep,ierr);CHKERRA(ierr)
154:   call NEPGetType(nep,ntype,ierr);CHKERRA(ierr)
155:   write(string,*) 'Solution method: ',ntype,'\n'
156:   call PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr);CHKERRA(ierr)
157:   call NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
158:   write(string,*) 'Number of requested eigenvalues:',nev,'\n'
159:   call PetscPrintf(PETSC_COMM_WORLD,trim(string),ierr);CHKERRA(ierr)

161: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: !     Display solution and clean up
163: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

165:   ! ** show detailed info unless -terse option is given by user
166:   call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
167:   if (terse) then
168:     call NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
169:   else
170:     call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
171:     call NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
172:     call NEPErrorView(nep,NEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
173:     call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
174:   end if

176:   if (split) then
177:     call MatDestroy(A(1),ierr);CHKERRA(ierr)
178:     call MatDestroy(A(2),ierr);CHKERRA(ierr)
179:     call FNDestroy(fn(1),ierr);CHKERRA(ierr)
180:     call FNDestroy(fn(2),ierr);CHKERRA(ierr)
181:   else
182:     call MatDestroy(F,ierr);CHKERRA(ierr)
183:     call MatDestroy(J,ierr);CHKERRA(ierr)
184:   end if
185:   call NEPDestroy(nep,ierr)
186:   call SlepcFinalize(ierr)

188: END PROGRAM main

190: ! --------------------------------------------------------------
191: !
192: !   FormFunction - Computes Function matrix  T(lambda)
193: !
194: SUBROUTINE FormFunction(nep,lambda,fun,B,ctx,ierr)
195: #include <slepc/finclude/slepcnep.h>
196:   use slepcnep
197:   implicit none

199:   NEP            :: nep
200:   PetscScalar    :: lambda,val(0:2),t
201:   Mat            :: fun,B
202:   PetscInt       :: ctx,i,n,col(0:2),Istart,Iend,Istart0,Iend0,one,two,three
203:   PetscErrorCode :: ierr
204:   PetscBool      :: FirstBlock=PETSC_FALSE, LastBlock=PETSC_FALSE

206:   one   = 1
207:   two   = 2
208:   three = 3

210:   ! ** Compute Function entries and insert into matrix
211:   t = sqrt(lambda)
212:   call MatGetSize(fun,n,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
213:   call MatGetOwnershipRange(fun,Istart,Iend,ierr);CHKERRA(ierr)
214:   if (Istart.eq.0) FirstBlock=PETSC_TRUE;
215:   if (Iend.eq.n) LastBlock=PETSC_TRUE;
216:   val(0)=1.0; val(1)=t-2.0; val(2)=1.0;

218:   Istart0 = Istart
219:   if (FirstBlock) Istart0 = Istart+1
220:   Iend0 = Iend
221:   if (LastBlock) Iend0 = Iend-1

223:   do i=Istart0,Iend0-1
224:      col(0) = i-1
225:      col(1) = i
226:      col(2) = i+1
227:      call MatSetValues(fun,one,i,three,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
228:   end do

230:   if (LastBlock) then
231:      i = n-1
232:      col(0) = n-2
233:      col(1) = n-1
234:      val(0) = 1.0
235:      val(1) = t-2.0
236:      call MatSetValues(fun,one,i,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
237:   end if

239:   if (FirstBlock) then
240:      i = 0
241:      col(0) = 0
242:      col(1) = 1
243:      val(0) = t-2.0
244:      val(1) = 1.0
245:      call MatSetValues(fun,one,i,two,col,val,INSERT_VALUES,ierr);CHKERRA(ierr)
246:   end if

248:   ! ** Assemble matrix
249:   call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
250:   call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
251:   call MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)
252:   call MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY,ierr);PetscCall(ierr)

254: END SUBROUTINE FormFunction

256: ! --------------------------------------------------------------
257: !
258: !   FormJacobian - Computes Jacobian matrix  T'(lambda)
259: !
260: SUBROUTINE FormJacobian(nep,lambda,jac,ctx,ierr)
261: #include <slepc/finclude/slepcnep.h>
262:   USE slepcnep
263:   implicit none

265:   NEP            :: nep
266:   PetscScalar    :: lambda,t
267:   Mat            :: jac
268:   PetscInt       :: ctx
269:   PetscErrorCode :: ierr
270:   Vec            :: d

272:   call MatCreateVecs(jac,d,PETSC_NULL_VEC,ierr);CHKERRA(ierr)
273:   t = 0.5/sqrt(lambda)
274:   call VecSet(d,t,ierr);CHKERRA(ierr)
275:   call MatDiagonalSet(jac,d,INSERT_VALUES,ierr);CHKERRA(ierr)
276:   calL VecDestroy(d,ierr);CHKERRA(ierr)

278: END SUBROUTINE FormJacobian

280: ! --------------------------------------------------------------
281: !
282: !  ComputeSingularities - This is a user-defined routine to compute maxnp
283: !  points (at most) in the complex plane where the function T(.) is not analytic.
284: !
285: !  In this case, we discretize the singularity region (-inf,0)~(-10e+6,-10e-6)
286: !
287: !  Input Parameters:
288: !    nep   - nonlinear eigensolver context
289: !    maxnp - on input number of requested points in the discretization (can be set)
290: !    xi    - computed values of the discretization
291: !    dummy - optional user-defined monitor context (unused here)
292: !
293: SUBROUTINE ComputeSingularities(nep,maxnp,xi,dummy,ierr)
294: #include <slepc/finclude/slepcnep.h>
295:   use slepcnep
296:   implicit none

298:   NEP            :: nep
299:   PetscInt       :: maxnp, dummy
300:   PetscScalar    :: xi(0:maxnp-1)
301:   PetscErrorCode :: ierr
302:   PetscReal      :: h
303:   PetscInt       :: i

305:   h = 11.0/real(maxnp-1)
306:   xi(0) = -1e-5
307:   xi(maxnp-1) = -1e+6
308:   do i=1,maxnp-2
309:      xi(i) = -10**(-5+h*i)
310:   end do

312: END SUBROUTINE ComputeSingularities

314: !/*TEST
315: !
316: !   test:
317: !      suffix: 1
318: !      args: -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
319: !      requires: !single
320: !      filter: sed -e "s/[+-]0\.0*i//g"
321: !
322: !   test:
323: !      suffix: 2
324: !      args: -split 0 -nep_nev 3 -nep_nleigs_interpolation_degree 90 -terse
325: !      requires: !single
326: !      filter: sed -e "s/[+-]0\.0*i//g"
327: !
328: !TEST*/