Actual source code: test10.c

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

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: static char help[] = "Tests multiple calls to NEPSolve() with different matrix size.\n\n"
 12:   "The command line options are:\n"
 13:   "  -n <n>, where <n> = number of grid subdivisions.\n"
 14:   "  -tau <tau>, where <tau> is the delay parameter.\n"
 15:   "  -split <0/1>, to select the split form in the problem definition (enabled by default).\n";

 17: /* Based on ex22.c (delay) */

 19: #include <slepcnep.h>

 21: /*
 22:    User-defined application context
 23: */
 24: typedef struct {
 25:   PetscScalar tau;
 26:   PetscReal   a;
 27: } ApplicationCtx;

 29: /*
 30:    Create problem matrices in split form
 31: */
 32: PetscErrorCode BuildSplitMatrices(PetscInt n,PetscReal a,Mat *Id,Mat *A,Mat *B)
 33: {
 34:   PetscInt       i,Istart,Iend;
 35:   PetscReal      h,xi;
 36:   PetscScalar    b;

 39:   h = PETSC_PI/(PetscReal)(n+1);

 41:   /* Identity matrix */
 42:   MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,Id);
 43:   MatSetOption(*Id,MAT_HERMITIAN,PETSC_TRUE);

 45:   /* A = 1/h^2*tridiag(1,-2,1) + a*I */
 46:   MatCreate(PETSC_COMM_WORLD,A);
 47:   MatSetSizes(*A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 48:   MatSetFromOptions(*A);
 49:   MatSetUp(*A);
 50:   MatGetOwnershipRange(*A,&Istart,&Iend);
 51:   for (i=Istart;i<Iend;i++) {
 52:     if (i>0) MatSetValue(*A,i,i-1,1.0/(h*h),INSERT_VALUES);
 53:     if (i<n-1) MatSetValue(*A,i,i+1,1.0/(h*h),INSERT_VALUES);
 54:     MatSetValue(*A,i,i,-2.0/(h*h)+a,INSERT_VALUES);
 55:   }
 56:   MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
 57:   MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
 58:   MatSetOption(*A,MAT_HERMITIAN,PETSC_TRUE);

 60:   /* B = diag(b(xi)) */
 61:   MatCreate(PETSC_COMM_WORLD,B);
 62:   MatSetSizes(*B,PETSC_DECIDE,PETSC_DECIDE,n,n);
 63:   MatSetFromOptions(*B);
 64:   MatSetUp(*B);
 65:   MatGetOwnershipRange(*B,&Istart,&Iend);
 66:   for (i=Istart;i<Iend;i++) {
 67:     xi = (i+1)*h;
 68:     b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
 69:     MatSetValue(*B,i,i,b,INSERT_VALUES);
 70:   }
 71:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
 72:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
 73:   MatSetOption(*B,MAT_HERMITIAN,PETSC_TRUE);
 74:   PetscFunctionReturn(0);
 75: }

 77: /*
 78:    Compute Function matrix  T(lambda)
 79: */
 80: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
 81: {
 82:   ApplicationCtx *user = (ApplicationCtx*)ctx;
 83:   PetscInt       i,n,Istart,Iend;
 84:   PetscReal      h,xi;
 85:   PetscScalar    b;

 88:   MatGetSize(fun,&n,NULL);
 89:   h = PETSC_PI/(PetscReal)(n+1);
 90:   MatGetOwnershipRange(fun,&Istart,&Iend);
 91:   for (i=Istart;i<Iend;i++) {
 92:     if (i>0) MatSetValue(fun,i,i-1,1.0/(h*h),INSERT_VALUES);
 93:     if (i<n-1) MatSetValue(fun,i,i+1,1.0/(h*h),INSERT_VALUES);
 94:     xi = (i+1)*h;
 95:     b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
 96:     MatSetValue(fun,i,i,-lambda-2.0/(h*h)+user->a+PetscExpScalar(-user->tau*lambda)*b,INSERT_VALUES);
 97:   }
 98:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 99:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
100:   if (fun != B) {
101:     MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
102:     MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
103:   }
104:   PetscFunctionReturn(0);
105: }

107: /*
108:    Compute Jacobian matrix  T'(lambda)
109: */
110: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
111: {
112:   ApplicationCtx *user = (ApplicationCtx*)ctx;
113:   PetscInt       i,n,Istart,Iend;
114:   PetscReal      h,xi;
115:   PetscScalar    b;

118:   MatGetSize(jac,&n,NULL);
119:   h = PETSC_PI/(PetscReal)(n+1);
120:   MatGetOwnershipRange(jac,&Istart,&Iend);
121:   for (i=Istart;i<Iend;i++) {
122:     xi = (i+1)*h;
123:     b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
124:     MatSetValue(jac,i,i,-1.0-user->tau*PetscExpScalar(-user->tau*lambda)*b,INSERT_VALUES);
125:   }
126:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
127:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
128:   PetscFunctionReturn(0);
129: }

131: int main(int argc,char **argv)
132: {
133:   NEP            nep;             /* nonlinear eigensolver context */
134:   Mat            Id,A,B,J,F;      /* problem matrices */
135:   FN             f1,f2,f3;        /* functions to define the nonlinear operator */
136:   ApplicationCtx ctx;             /* user-defined context */
137:   Mat            mats[3];
138:   FN             funs[3];
139:   PetscScalar    coeffs[2];
140:   PetscInt       n=128;
141:   PetscReal      tau=0.001,a=20;
142:   PetscBool      split=PETSC_TRUE;

144:   SlepcInitialize(&argc,&argv,(char*)0,help);
145:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
146:   PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL);
147:   PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL);
148:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:               Create nonlinear eigensolver and set options
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

154:   NEPCreate(PETSC_COMM_WORLD,&nep);
155:   NEPSetTolerances(nep,1e-9,PETSC_DEFAULT);

157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158:                       First solve
159:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

161:   if (split) {
162:     BuildSplitMatrices(n,a,&Id,&A,&B);
163:     /* f1=-lambda */
164:     FNCreate(PETSC_COMM_WORLD,&f1);
165:     FNSetType(f1,FNRATIONAL);
166:     coeffs[0] = -1.0; coeffs[1] = 0.0;
167:     FNRationalSetNumerator(f1,2,coeffs);
168:     /* f2=1.0 */
169:     FNCreate(PETSC_COMM_WORLD,&f2);
170:     FNSetType(f2,FNRATIONAL);
171:     coeffs[0] = 1.0;
172:     FNRationalSetNumerator(f2,1,coeffs);
173:     /* f3=exp(-tau*lambda) */
174:     FNCreate(PETSC_COMM_WORLD,&f3);
175:     FNSetType(f3,FNEXP);
176:     FNSetScale(f3,-tau,1.0);
177:     mats[0] = A;  funs[0] = f2;
178:     mats[1] = Id; funs[1] = f1;
179:     mats[2] = B;  funs[2] = f3;
180:     NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN);
181:   } else {
182:     /* callback form  */
183:     ctx.tau = tau;
184:     ctx.a   = a;
185:     MatCreate(PETSC_COMM_WORLD,&F);
186:     MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
187:     MatSetFromOptions(F);
188:     MatSeqAIJSetPreallocation(F,3,NULL);
189:     MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
190:     MatSetUp(F);
191:     NEPSetFunction(nep,F,F,FormFunction,&ctx);
192:     MatCreate(PETSC_COMM_WORLD,&J);
193:     MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
194:     MatSetFromOptions(J);
195:     MatSeqAIJSetPreallocation(J,3,NULL);
196:     MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
197:     MatSetUp(J);
198:     NEPSetJacobian(nep,J,FormJacobian,&ctx);
199:   }

201:   /* Set solver parameters at runtime */
202:   NEPSetFromOptions(nep);

204:   /* Solve the eigensystem */
205:   NEPSolve(nep);
206:   NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);

208:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209:                Second solve, with problem matrices of size 2*n
210:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

212:   n *= 2;
213:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau);
214:   if (split) {
215:     MatDestroy(&Id);
216:     MatDestroy(&A);
217:     MatDestroy(&B);
218:     BuildSplitMatrices(n,a,&Id,&A,&B);
219:     mats[0] = A;
220:     mats[1] = Id;
221:     mats[2] = B;
222:     NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN);
223:   } else {
224:     /* callback form  */
225:     MatDestroy(&F);
226:     MatDestroy(&J);
227:     MatCreate(PETSC_COMM_WORLD,&F);
228:     MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
229:     MatSetFromOptions(F);
230:     MatSeqAIJSetPreallocation(F,3,NULL);
231:     MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
232:     MatSetUp(F);
233:     NEPSetFunction(nep,F,F,FormFunction,&ctx);
234:     MatCreate(PETSC_COMM_WORLD,&J);
235:     MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
236:     MatSetFromOptions(J);
237:     MatSeqAIJSetPreallocation(J,3,NULL);
238:     MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
239:     MatSetUp(J);
240:     NEPSetJacobian(nep,J,FormJacobian,&ctx);
241:   }

243:   /* Solve the eigensystem */
244:   NEPSolve(nep);
245:   NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);

247:   NEPDestroy(&nep);
248:   if (split) {
249:     MatDestroy(&Id);
250:     MatDestroy(&A);
251:     MatDestroy(&B);
252:     FNDestroy(&f1);
253:     FNDestroy(&f2);
254:     FNDestroy(&f3);
255:   } else {
256:     MatDestroy(&F);
257:     MatDestroy(&J);
258:   }
259:   SlepcFinalize();
260:   return 0;
261: }

263: /*TEST

265:    testset:
266:       nsize: 2
267:       requires: !single
268:       output_file: output/test10_1.out
269:       test:
270:          suffix: 1
271:          args: -nep_type narnoldi -nep_target 0.55
272:       test:
273:          suffix: 1_rii
274:          args: -nep_type rii -nep_target 0.55 -nep_rii_hermitian -split {{0 1}}
275:       test:
276:          suffix: 1_narnoldi
277:          args: -nep_type narnoldi -nep_target 0.55 -nep_narnoldi_lag_preconditioner 2
278:       test:
279:          suffix: 1_slp
280:          args: -nep_type slp -nep_slp_st_pc_type redundant -split {{0 1}}
281:       test:
282:          suffix: 1_interpol
283:          args: -nep_type interpol -rg_type interval -rg_interval_endpoints .5,1,-.1,.1 -nep_target .7 -nep_interpol_st_pc_type redundant
284:       test:
285:          suffix: 1_narnoldi_sync
286:          args: -nep_type narnoldi -ds_parallel synchronized

288:    testset:
289:       args: -nep_nev 2 -rg_type interval -rg_interval_endpoints .5,15,-.1,.1 -nep_target .7
290:       requires: !single
291:       output_file: output/test10_2.out
292:       filter: sed -e "s/[+-]0\.0*i//g"
293:       test:
294:          suffix: 2_interpol
295:          args: -nep_type interpol -nep_interpol_pep_type jd -nep_interpol_st_pc_type sor
296:       test:
297:          suffix: 2_nleigs
298:          args: -nep_type nleigs -split {{0 1}}
299:          requires: complex
300:       test:
301:          suffix: 2_nleigs_real
302:          args: -nep_type nleigs -rg_interval_endpoints .5,15 -split {{0 1}}
303:          requires: !complex

305:    test:
306:       suffix: 3
307:       requires: complex !single
308:       args: -nep_type ciss -rg_type ellipse -rg_ellipse_center 10 -rg_ellipse_radius 9.5 -rg_ellipse_vscale 0.1 -split {{0 1}}
309:       timeoutfactor: 2

311: TEST*/