Actual source code: test10.c

slepc-3.18.3 2023-03-24
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:   return 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:   return 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:   return 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;

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

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

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

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

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

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

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

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

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

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

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

264: /*TEST

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

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

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

312: TEST*/