Actual source code: ex22f90.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> ./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 MatCreate(PETSC_COMM_WORLD,Id,ierr);CHKERRA(ierr)
89: call MatSetSizes(Id,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
90: call MatSetFromOptions(Id,ierr);CHKERRA(ierr)
91: call MatSetUp(Id,ierr);CHKERRA(ierr)
92: call MatGetOwnershipRange(Id,Istart,Iend,ierr);CHKERRA(ierr)
93: do i=Istart,Iend-1
94: call MatSetValue(Id,i,i,one,INSERT_VALUES,ierr);CHKERRA(ierr)
95: end do
96: call MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
97: call MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
98: call MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)
100: ! ** A = 1/h^2*tridiag(1,-2,1) + aa*I
101: call MatCreate(PETSC_COMM_WORLD,A,ierr);CHKERRA(ierr)
102: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
103: call MatSetFromOptions(A,ierr);CHKERRA(ierr)
104: call MatSetUp(A,ierr);CHKERRA(ierr)
105: call MatGetOwnershipRange(A,Istart,Iend,ierr);CHKERRA(ierr)
106: coeffs(1) = 1.0/(h*h)
107: coeffs(2) = -2.0/(h*h)+aa
108: do i=Istart,Iend-1
109: if (i .gt. 0) then
110: call MatSetValue(A,i,i-1,coeffs(1),INSERT_VALUES,ierr);CHKERRA(ierr)
111: endif
112: if (i .lt. n-1) then
113: call MatSetValue(A,i,i+1,coeffs(1),INSERT_VALUES,ierr);CHKERRA(ierr)
114: endif
115: call MatSetValue(A,i,i,coeffs(2),INSERT_VALUES,ierr);CHKERRA(ierr)
116: end do
117: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
118: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
119: call MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)
121: ! ** B = diag(bb(xi))
122: call MatCreate(PETSC_COMM_WORLD,B,ierr);CHKERRA(ierr)
123: call MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n,ierr);CHKERRA(ierr)
124: call MatSetFromOptions(B,ierr);CHKERRA(ierr)
125: call MatSetUp(B,ierr);CHKERRA(ierr)
126: call MatGetOwnershipRange(B,Istart,Iend,ierr);CHKERRA(ierr)
127: do i=Istart,Iend-1
128: xi = (i+1)*h
129: bb = -4.1+xi*(1.0-exp(xi-PETSC_PI))
130: call MatSetValue(B,i,i,bb,INSERT_VALUES,ierr);CHKERRA(ierr)
131: end do
132: call MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
133: call MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
134: call MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE,ierr);CHKERRA(ierr)
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: ! Create problem functions, f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: call FNCreate(PETSC_COMM_WORLD,f1,ierr);CHKERRA(ierr)
141: call FNSetType(f1,FNRATIONAL,ierr);CHKERRA(ierr)
142: k = 2
143: coeffs(1) = -1.0
144: coeffs(2) = 0.0
145: call FNRationalSetNumerator(f1,k,coeffs,ierr);CHKERRA(ierr)
147: call FNCreate(PETSC_COMM_WORLD,f2,ierr);CHKERRA(ierr)
148: call FNSetType(f2,FNRATIONAL,ierr);CHKERRA(ierr)
149: k = 1
150: coeffs(1) = 1.0
151: call FNRationalSetNumerator(f2,k,coeffs,ierr);CHKERRA(ierr)
153: call FNCreate(PETSC_COMM_WORLD,f3,ierr);CHKERRA(ierr)
154: call FNSetType(f3,FNEXP,ierr);CHKERRA(ierr)
155: scal = -tau
156: call FNSetScale(f3,scal,one,ierr);CHKERRA(ierr)
158: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: ! Create the eigensolver and set various options
160: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: ! ** Create eigensolver context
163: call NEPCreate(PETSC_COMM_WORLD,nep,ierr);CHKERRA(ierr)
165: ! ** Set the split operator. Note that A is passed first so that
166: ! ** SUBSET_NONZERO_PATTERN can be used
167: k = 3
168: mats(1) = A
169: mats(2) = Id
170: mats(3) = B
171: funs(1) = f2
172: funs(2) = f1
173: funs(3) = f3
174: call NEPSetSplitOperator(nep,k,mats,funs,SUBSET_NONZERO_PATTERN,ierr);CHKERRA(ierr)
175: call NEPSetProblemType(nep,NEP_GENERAL,ierr);CHKERRA(ierr)
177: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: ! Customize nonlinear solver; set runtime options
179: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: tol = 1e-9
182: call NEPSetTolerances(nep,tol,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
183: k = 1
184: call NEPSetDimensions(nep,k,PETSC_DEFAULT_INTEGER,PETSC_DEFAULT_INTEGER,ierr);CHKERRA(ierr)
185: k = 0
186: call NEPRIISetLagPreconditioner(nep,k,ierr);CHKERRA(ierr)
188: ! ** Set solver parameters at runtime
189: call NEPSetFromOptions(nep,ierr);CHKERRA(ierr)
191: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: ! Solve the eigensystem
193: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: call NEPSolve(nep,ierr);CHKERRA(ierr)
197: ! ** Optional: Get some information from the solver and display it
198: call NEPGetType(nep,tname,ierr);CHKERRA(ierr)
199: if (rank .eq. 0) then
200: write(*,120) tname
201: endif
202: 120 format (' Solution method: ',A)
203: call NEPGetDimensions(nep,nev,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRA(ierr)
204: if (rank .eq. 0) then
205: write(*,130) nev
206: endif
207: 130 format (' Number of requested eigenvalues:',I4)
209: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: ! Display solution and clean up
211: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: ! ** show detailed info unless -terse option is given by user
214: call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-terse',terse,ierr);CHKERRA(ierr)
215: if (terse) then
216: call NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_NULL_VIEWER,ierr);CHKERRA(ierr)
217: else
218: call PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL,ierr);CHKERRA(ierr)
219: call NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
220: call NEPErrorView(nep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
221: call PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRA(ierr)
222: endif
223: call NEPDestroy(nep,ierr);CHKERRA(ierr)
224: call MatDestroy(Id,ierr);CHKERRA(ierr)
225: call MatDestroy(A,ierr);CHKERRA(ierr)
226: call MatDestroy(B,ierr);CHKERRA(ierr)
227: call FNDestroy(f1,ierr);CHKERRA(ierr)
228: call FNDestroy(f2,ierr);CHKERRA(ierr)
229: call FNDestroy(f3,ierr);CHKERRA(ierr)
230: call SlepcFinalize(ierr)
231: end
233: !/*TEST
234: !
235: ! test:
236: ! suffix: 1
237: ! args: -terse
238: ! requires: !single
239: !
240: !TEST*/