Actual source code: ex22f90.F90

slepc-3.17.2 2022-08-09
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> ./ex22f90 [-n <n>] [-tau <tau>] [SLEPc opts]
 11: !
 12: !  Description: Delay differential equation. Fortran90 equivalent of ex22.c
 13: !
 14: !  The command line options are:
 15: !    -n <n>, where <n> = number of grid subdivisions
 16: !    -tau <tau>, where <tau> = delay parameter
 17: !
 18: ! ----------------------------------------------------------------------
 19: !  Solve parabolic partial differential equation with time delay tau
 20: !
 21: !           u_t = u_xx + aa*u(t) + bb*u(t-tau)
 22: !           u(0,t) = u(pi,t) = 0
 23: !
 24: !  with aa = 20 and bb(x) = -4.1+x*(1-exp(x-pi)).
 25: !
 26: !  Discretization leads to a DDE of dimension n
 27: !
 28: !           -u' = A*u(t) + B*u(t-tau)
 29: !
 30: !  which results in the nonlinear eigenproblem
 31: !
 32: !           (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
 33: ! ----------------------------------------------------------------------
 34: !
 35:       program main
 36: #include <slepc/finclude/slepcnep.h>
 37:       use slepcnep
 38:       implicit none

 40: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 41: !     Declarations
 42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43: !
 44: !  Variables:
 45: !     nep       nonlinear eigensolver context
 46: !     Id,A,B    problem matrices
 47: !     f1,f2,f3  functions to define the nonlinear operator

 49:       Mat            Id, A, B, mats(3)
 50:       FN             f1, f2, f3, funs(3)
 51:       NEP            nep
 52:       NEPType        tname
 53:       PetscScalar    one, bb, coeffs(2), scal
 54:       PetscReal      tau, h, aa, xi, tol
 55:       PetscInt       n, i, k, nev, Istart, Iend
 56:       PetscMPIInt    rank
 57:       PetscErrorCode ierr
 58:       PetscBool      flg, terse

 60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61: !     Beginning of program
 62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 64:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 65:       if (ierr .ne. 0) then
 66:         print*,'SlepcInitialize failed'
 67:         stop
 68:       endif
 69:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr);CHKERRA(ierr)
 70:       n = 128
 71:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr);CHKERRA(ierr)
 72:       tau = 0.001
 73:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-tau',tau,flg,ierr);CHKERRA(ierr)
 74:       if (rank .eq. 0) then
 75:         write(*,100) n, tau
 76:       endif
 77:  100  format (/'Delay Eigenproblem, n =',I4,', tau =',F6.3)

 79:       one = 1.0
 80:       aa = 20.0
 81:       h = PETSC_PI/real(n+1)

 83: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84: !     Create problem matrices
 85: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 87: !     ** Id is the identity matrix
 88:       call MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,one,Id,ierr);CHKERRA(ierr)
 89:       call MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)

 91: !     ** A = 1/h^2*tridiag(1,-2,1) + aa*I
 92:       call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
 93:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
 94:       call MatSetFromOptions(A,ierr);CHKERRA(ierr)
 95:       call MatSetUp(A,ierr);CHKERRA(ierr)
 96:       call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
 97:       coeffs(1) = 1.0/(h*h)
 98:       coeffs(2) = -2.0/(h*h)+aa
 99:       do i=Istart,Iend-1
100:         if (i .gt. 0) then
101:           call MatSetValue(A,i,i-1,coeffs(1),INSERT_VALUES,ierr);CHKERRA(ierr)
102:         endif
103:         if (i .lt. n-1) then
104:           call MatSetValue(A,i,i+1,coeffs(1),INSERT_VALUES,ierr);CHKERRA(ierr)
105:         endif
106:         call MatSetValue(A,i,i,coeffs(2),INSERT_VALUES,ierr);CHKERRA(ierr)
107:       end do
108:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
109:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
110:       call MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)

112: !     ** B = diag(bb(xi))
113:       call MatCreate(PETSC_COMM_WORLD,B,ierr);CHKERRA(ierr)
114:       call MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
115:       call MatSetFromOptions(B,ierr);CHKERRA(ierr)
116:       call MatSetUp(B,ierr);CHKERRA(ierr)
117:       call MatGetOwnershipRange(B,Istart,Iend,ierr);CHKERRA(ierr)
118:       do i=Istart,Iend-1
119:         xi = (i+1)*h
120:         bb  = -4.1+xi*(1.0-exp(xi-PETSC_PI))
121:         call MatSetValue(B,i,i,bb,INSERT_VALUES,ierr);CHKERRA(ierr)
122:       end do
123:       call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
124:       call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
125:       call MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)

127: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: !     Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

131:       call FNCreate(PETSC_COMM_WORLD,f1,ierr);CHKERRA(ierr)
132:       call FNSetType(f1,FNRATIONAL,ierr);CHKERRA(ierr)
133:       k = 2
134:       coeffs(1) = -1.0
135:       coeffs(2) = 0.0
136:       call FNRationalSetNumerator(f1,k,coeffs,ierr);CHKERRA(ierr)

138:       call FNCreate(PETSC_COMM_WORLD,f2,ierr);CHKERRA(ierr)
139:       call FNSetType(f2,FNRATIONAL,ierr);CHKERRA(ierr)
140:       k = 1
141:       coeffs(1) = 1.0
142:       call FNRationalSetNumerator(f2,k,coeffs,ierr);CHKERRA(ierr)

144:       call FNCreate(PETSC_COMM_WORLD,f3,ierr);CHKERRA(ierr)
145:       call FNSetType(f3,FNEXP,ierr);CHKERRA(ierr)
146:       scal = -tau
147:       call FNSetScale(f3,scal,one,ierr);CHKERRA(ierr)

149: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: !     Create the eigensolver and set various options
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

153: !     ** Create eigensolver context
154:       call NEPCreate(PETSC_COMM_WORLD,nep,ierr);CHKERRA(ierr)

156: !     ** Set the split operator. Note that A is passed first so that
157: !     ** SUBSET_NONZERO_PATTERN can be used
158:       k = 3
159:       mats(1) = A
160:       mats(2) = Id
161:       mats(3) = B
162:       funs(1) = f2
163:       funs(2) = f1
164:       funs(3) = f3
165:       call NEPSetSplitOperator(nep,k,mats,funs,SUBSET_NONZERO_PATTERN,ierr);CHKERRA(ierr)
166:       call NEPSetProblemType(nep,NEP_GENERAL,ierr);CHKERRA(ierr)

168: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: !     Customize nonlinear solver; set runtime options
170: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

172:       tol = 1e-9
173:       call NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
174:       k = 1
175:       call NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
176:       k = 0
177:       call NEPRIISetLagPreconditioner(nep,k,ierr);CHKERRA(ierr)

179: !     ** Set solver parameters at runtime
180:       call NEPSetFromOptions(nep,ierr);CHKERRA(ierr)

182: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: !     Solve the eigensystem
184: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

186:       call NEPSolve(nep,ierr);CHKERRA(ierr)

188: !     ** Optional: Get some information from the solver and display it
189:       call NEPGetType(nep,tname,ierr);CHKERRA(ierr)
190:       if (rank .eq. 0) then
191:         write(*,120) tname
192:       endif
193:  120  format (' Solution method: ',A)
194:       call NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
195:       if (rank .eq. 0) then
196:         write(*,130) nev
197:       endif
198:  130  format (' Number of requested eigenvalues:',I4)

200: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: !     Display solution and clean up
202: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

204: !     ** show detailed info unless -terse option is given by user
205:       call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
206:       if (terse) then
207:         call NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
208:       else
209:         call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
210:         call NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
211:         call NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
212:         call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
213:       endif
214:       call NEPDestroy(nep,ierr);CHKERRA(ierr)
215:       call MatDestroy(Id,ierr);CHKERRA(ierr)
216:       call MatDestroy(A,ierr);CHKERRA(ierr)
217:       call MatDestroy(B,ierr);CHKERRA(ierr)
218:       call FNDestroy(f1,ierr);CHKERRA(ierr)
219:       call FNDestroy(f2,ierr);CHKERRA(ierr)
220:       call FNDestroy(f3,ierr);CHKERRA(ierr)
221:       call SlepcFinalize(ierr)
222:       end

224: !/*TEST
225: !
226: !   test:
227: !      suffix: 1
228: !      args: -terse
229: !      requires: !single
230: !      filter: sed -e "s/[+-]0\.0*i//g"
231: !
232: !TEST*/