Actual source code: ex22.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2013, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Delay differential equation.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = number of grid subdivisions.\n"
25: " -tau <tau>, where <tau> is the delay parameter.\n\n";
27: /*
28: Solve parabolic partial differential equation with time delay tau
30: u_t = u_xx + a*u(t) + b*u(t-tau)
31: u(0,t) = u(pi,t) = 0
33: with a = 20 and b(x) = -4.1+x*(1-exp(x-pi)).
35: Discretization leads to a DDE of dimension n
37: -u' = A*u(t) + B*u(t-tau)
39: which results in the nonlinear eigenproblem
41: (-lambda*I + A + exp(-tau*lambda)*B)*u = 0
42: */
44: #include <slepcnep.h>
48: int main(int argc,char **argv)
49: {
50: NEP nep; /* nonlinear eigensolver context */
51: PetscScalar lambda; /* eigenvalue */
52: Mat Id,A,B; /* problem matrices */
53: FN f1,f2,f3; /* functions to define the nonlinear operator */
54: Mat mats[3];
55: FN funs[3];
56: NEPType type;
57: PetscScalar value[3],coeffs[2],b;
58: PetscInt n=128,nev,Istart,Iend,col[3],i,its,nconv;
59: PetscReal tau=0.001,h,a=20,xi,re,im,norm;
60: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
63: SlepcInitialize(&argc,&argv,(char*)0,help);
64: PetscOptionsGetInt(NULL,"-n",&n,NULL);
65: PetscOptionsGetReal(NULL,"-tau",&tau,NULL);
66: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%D, tau=%G\n\n",n,tau);
67: h = PETSC_PI/(PetscReal)(n+1);
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Create nonlinear eigensolver context
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: NEPCreate(PETSC_COMM_WORLD,&nep);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create problem matrices and coefficient functions. Pass them to NEP
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: /*
80: Identity matrix
81: */
82: MatCreate(PETSC_COMM_WORLD,&Id);
83: MatSetSizes(Id,PETSC_DECIDE,PETSC_DECIDE,n,n);
84: MatSetFromOptions(Id);
85: MatSetUp(Id);
86: MatAssemblyBegin(Id,MAT_FINAL_ASSEMBLY);
87: MatAssemblyEnd(Id,MAT_FINAL_ASSEMBLY);
88: MatShift(Id,1.0);
89: MatSetOption(Id,MAT_HERMITIAN,PETSC_TRUE);
91: /*
92: A = 1/h^2*tridiag(1,-2,1) + a*I
93: */
94: MatCreate(PETSC_COMM_WORLD,&A);
95: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
96: MatSetFromOptions(A);
97: MatSetUp(A);
98: MatGetOwnershipRange(A,&Istart,&Iend);
99: if (Istart==0) FirstBlock=PETSC_TRUE;
100: if (Iend==n) LastBlock=PETSC_TRUE;
101: value[0]=1.0/(h*h); value[1]=-2.0/(h*h)+a; value[2]=1.0/(h*h);
102: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
103: col[0]=i-1; col[1]=i; col[2]=i+1;
104: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
105: }
106: if (LastBlock) {
107: i=n-1; col[0]=n-2; col[1]=n-1;
108: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
109: }
110: if (FirstBlock) {
111: i=0; col[0]=0; col[1]=1; value[0]=-2.0/(h*h)+a; value[1]=1.0/(h*h);
112: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
113: }
114: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
116: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
118: /*
119: B = diag(b(xi))
120: */
121: MatCreate(PETSC_COMM_WORLD,&B);
122: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
123: MatSetFromOptions(B);
124: MatSetUp(B);
125: MatGetOwnershipRange(B,&Istart,&Iend);
126: for (i=Istart;i<Iend;i++) {
127: xi = (i+1)*h;
128: b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
129: MatSetValues(B,1,&i,1,&i,&b,INSERT_VALUES);
130: }
131: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
132: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
133: MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);
135: /*
136: Functions: f1=-lambda, f2=1.0, f3=exp(-tau*lambda)
137: */
138: FNCreate(PETSC_COMM_WORLD,&f1);
139: FNSetType(f1,FNRATIONAL);
140: coeffs[0] = -1.0; coeffs[1] = 0.0;
141: FNSetParameters(f1,2,coeffs,0,NULL);
143: FNCreate(PETSC_COMM_WORLD,&f2);
144: FNSetType(f2,FNRATIONAL);
145: coeffs[0] = 1.0;
146: FNSetParameters(f2,1,coeffs,0,NULL);
148: FNCreate(PETSC_COMM_WORLD,&f3);
149: FNSetType(f3,FNEXP);
150: coeffs[0] = -tau;
151: FNSetParameters(f3,1,coeffs,0,NULL);
153: /*
154: Set the split operator. Note that A is passed first so that
155: SUBSET_NONZERO_PATTERN can be used
156: */
157: mats[0] = A; funs[0] = f2;
158: mats[1] = Id; funs[1] = f1;
159: mats[2] = B; funs[2] = f3;
160: NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN);
162: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
163: Customize nonlinear solver; set runtime options
164: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: NEPSetTolerances(nep,0,1e-9,0,0,0);
167: NEPSetDimensions(nep,1,0,0);
168: NEPSetLagPreconditioner(nep,0);
170: /*
171: Set solver parameters at runtime
172: */
173: NEPSetFromOptions(nep);
175: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: Solve the eigensystem
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: NEPSolve(nep);
180: NEPGetIterationNumber(nep,&its);
181: PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %D\n\n",its);
183: /*
184: Optional: Get some information from the solver and display it
185: */
186: NEPGetType(nep,&type);
187: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
188: NEPGetDimensions(nep,&nev,NULL,NULL);
189: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
191: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: Display solution and clean up
193: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
195: /*
196: Get number of converged approximate eigenpairs
197: */
198: NEPGetConverged(nep,&nconv);
199: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);
201: if (nconv>0) {
202: /*
203: Display eigenvalues and relative errors
204: */
205: PetscPrintf(PETSC_COMM_WORLD,
206: " k ||T(k)x||\n"
207: " ----------------- ------------------\n");
208: for (i=0;i<nconv;i++) {
209: NEPGetEigenpair(nep,i,&lambda,NULL);
210: NEPComputeRelativeError(nep,i,&norm);
211: #if defined(PETSC_USE_COMPLEX)
212: re = PetscRealPart(lambda);
213: im = PetscImaginaryPart(lambda);
214: #else
215: re = lambda;
216: im = 0.0;
217: #endif
218: if (im!=0.0) {
219: PetscPrintf(PETSC_COMM_WORLD," %9F%+9F j %12G\n",re,im,norm);
220: } else {
221: PetscPrintf(PETSC_COMM_WORLD," %12F %12G\n",re,norm);
222: }
223: }
224: PetscPrintf(PETSC_COMM_WORLD,"\n");
225: }
227: NEPDestroy(&nep);
228: MatDestroy(&Id);
229: MatDestroy(&A);
230: MatDestroy(&B);
231: FNDestroy(&f1);
232: FNDestroy(&f2);
233: FNDestroy(&f3);
234: SlepcFinalize();
235: return 0;
236: }