Actual source code: ex20.c

slepc-3.17.1 2022-04-11
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[] = "Simple 1-D nonlinear eigenproblem.\n\n"
 12:   "The command line options are:\n"
 13:   "  -n <n>, where <n> = number of grid subdivisions.\n"
 14:   "  -draw_sol, to draw the computed solution.\n\n";

 16: /*
 17:    Solve 1-D PDE
 18:             -u'' = lambda*u
 19:    on [0,1] subject to
 20:             u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
 21: */

 23: #include <slepcnep.h>

 25: /*
 26:    User-defined routines
 27: */
 28: PetscErrorCode FormInitialGuess(Vec);
 29: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
 30: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
 31: PetscErrorCode CheckSolution(PetscScalar,Vec,PetscReal*,void*);
 32: PetscErrorCode FixSign(Vec);

 34: /*
 35:    User-defined application context
 36: */
 37: typedef struct {
 38:   PetscScalar kappa;   /* ratio between stiffness of spring and attached mass */
 39:   PetscReal   h;       /* mesh spacing */
 40: } ApplicationCtx;

 42: int main(int argc,char **argv)
 43: {
 44:   NEP            nep;             /* nonlinear eigensolver context */
 45:   Vec            x;               /* eigenvector */
 46:   PetscScalar    lambda;          /* eigenvalue */
 47:   Mat            F,J;             /* Function and Jacobian matrices */
 48:   ApplicationCtx ctx;             /* user-defined context */
 49:   NEPType        type;
 50:   PetscInt       n=128,nev,i,its,maxit,nconv;
 51:   PetscReal      re,im,tol,norm,error;
 52:   PetscBool      draw_sol=PETSC_FALSE;

 54:   SlepcInitialize(&argc,&argv,(char*)0,help);
 55:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 56:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n);
 57:   ctx.h = 1.0/(PetscReal)n;
 58:   ctx.kappa = 1.0;
 59:   PetscOptionsGetBool(NULL,NULL,"-draw_sol",&draw_sol,NULL);

 61:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 62:      Create nonlinear eigensolver context
 63:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 65:   NEPCreate(PETSC_COMM_WORLD,&nep);

 67:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68:      Create matrix data structure; set Function evaluation routine
 69:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 71:   MatCreate(PETSC_COMM_WORLD,&F);
 72:   MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
 73:   MatSetFromOptions(F);
 74:   MatSeqAIJSetPreallocation(F,3,NULL);
 75:   MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
 76:   MatSetUp(F);

 78:   /*
 79:      Set Function matrix data structure and default Function evaluation
 80:      routine
 81:   */
 82:   NEPSetFunction(nep,F,F,FormFunction,&ctx);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Create matrix data structure; set Jacobian evaluation routine
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 88:   MatCreate(PETSC_COMM_WORLD,&J);
 89:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
 90:   MatSetFromOptions(J);
 91:   MatSeqAIJSetPreallocation(J,3,NULL);
 92:   MatMPIAIJSetPreallocation(J,3,NULL,1,NULL);
 93:   MatSetUp(J);

 95:   /*
 96:      Set Jacobian matrix data structure and default Jacobian evaluation
 97:      routine
 98:   */
 99:   NEPSetJacobian(nep,J,FormJacobian,&ctx);

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Customize nonlinear solver; set runtime options
103:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

105:   NEPSetTolerances(nep,1e-9,PETSC_DEFAULT);
106:   NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT);

108:   /*
109:      Set solver parameters at runtime
110:   */
111:   NEPSetFromOptions(nep);

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:                       Initialize application
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

117:   /*
118:      Evaluate initial guess
119:   */
120:   MatCreateVecs(F,&x,NULL);
121:   FormInitialGuess(x);
122:   NEPSetInitialSpace(nep,1,&x);

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:                       Solve the eigensystem
126:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

128:   NEPSolve(nep);
129:   NEPGetIterationNumber(nep,&its);
130:   PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %" PetscInt_FMT "\n\n",its);

132:   /*
133:      Optional: Get some information from the solver and display it
134:   */
135:   NEPGetType(nep,&type);
136:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
137:   NEPGetDimensions(nep,&nev,NULL,NULL);
138:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
139:   NEPGetTolerances(nep,&tol,&maxit);
140:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit);

142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143:                     Display solution and clean up
144:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

146:   /*
147:      Get number of converged approximate eigenpairs
148:   */
149:   NEPGetConverged(nep,&nconv);
150:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %" PetscInt_FMT "\n\n",nconv);

152:   if (nconv>0) {
153:     /*
154:        Display eigenvalues and relative errors
155:     */
156:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,
157:          "           k              ||T(k)x||           error\n"
158:          "   ----------------- ------------------ ------------------\n"));
159:     for (i=0;i<nconv;i++) {
160:       /*
161:         Get converged eigenpairs (in this example they are always real)
162:       */
163:       NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL);
164:       FixSign(x);
165:       /*
166:          Compute residual norm and error
167:       */
168:       NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm);
169:       CheckSolution(lambda,x,&error,&ctx);

171: #if defined(PETSC_USE_COMPLEX)
172:       re = PetscRealPart(lambda);
173:       im = PetscImaginaryPart(lambda);
174: #else
175:       re = lambda;
176:       im = 0.0;
177: #endif
178:       if (im!=0.0) PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g     %12g\n",(double)re,(double)im,(double)norm,(double)error);
179:       else PetscPrintf(PETSC_COMM_WORLD,"   %12f         %12g     %12g\n",(double)re,(double)norm,(double)error);
180:       if (draw_sol) {
181:         PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
182:         VecView(x,PETSC_VIEWER_DRAW_WORLD);
183:       }
184:     }
185:     PetscPrintf(PETSC_COMM_WORLD,"\n");
186:   }

188:   NEPDestroy(&nep);
189:   MatDestroy(&F);
190:   MatDestroy(&J);
191:   VecDestroy(&x);
192:   SlepcFinalize();
193:   return 0;
194: }

196: /* ------------------------------------------------------------------- */
197: /*
198:    FormInitialGuess - Computes initial guess.

200:    Input/Output Parameter:
201: .  x - the solution vector
202: */
203: PetscErrorCode FormInitialGuess(Vec x)
204: {
206:   VecSet(x,1.0);
207:   PetscFunctionReturn(0);
208: }

210: /* ------------------------------------------------------------------- */
211: /*
212:    FormFunction - Computes Function matrix  T(lambda)

214:    Input Parameters:
215: .  nep    - the NEP context
216: .  lambda - the scalar argument
217: .  ctx    - optional user-defined context, as set by NEPSetFunction()

219:    Output Parameters:
220: .  fun - Function matrix
221: .  B   - optionally different preconditioning matrix
222: */
223: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
224: {
225:   ApplicationCtx *user = (ApplicationCtx*)ctx;
226:   PetscScalar    A[3],c,d;
227:   PetscReal      h;
228:   PetscInt       i,n,j[3],Istart,Iend;
229:   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;

232:   /*
233:      Compute Function entries and insert into matrix
234:   */
235:   MatGetSize(fun,&n,NULL);
236:   MatGetOwnershipRange(fun,&Istart,&Iend);
237:   if (Istart==0) FirstBlock=PETSC_TRUE;
238:   if (Iend==n) LastBlock=PETSC_TRUE;
239:   h = user->h;
240:   c = user->kappa/(lambda-user->kappa);
241:   d = n;

243:   /*
244:      Interior grid points
245:   */
246:   for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
247:     j[0] = i-1; j[1] = i; j[2] = i+1;
248:     A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
249:     MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
250:   }

252:   /*
253:      Boundary points
254:   */
255:   if (FirstBlock) {
256:     i = 0;
257:     j[0] = 0; j[1] = 1;
258:     A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
259:     MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
260:   }

262:   if (LastBlock) {
263:     i = n-1;
264:     j[0] = n-2; j[1] = n-1;
265:     A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
266:     MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
267:   }

269:   /*
270:      Assemble matrix
271:   */
272:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
273:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
274:   if (fun != B) {
275:     MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
276:     MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
277:   }
278:   PetscFunctionReturn(0);
279: }

281: /* ------------------------------------------------------------------- */
282: /*
283:    FormJacobian - Computes Jacobian matrix  T'(lambda)

285:    Input Parameters:
286: .  nep    - the NEP context
287: .  lambda - the scalar argument
288: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

290:    Output Parameters:
291: .  jac - Jacobian matrix
292: */
293: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
294: {
295:   ApplicationCtx *user = (ApplicationCtx*)ctx;
296:   PetscScalar    A[3],c;
297:   PetscReal      h;
298:   PetscInt       i,n,j[3],Istart,Iend;
299:   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;

302:   /*
303:      Compute Jacobian entries and insert into matrix
304:   */
305:   MatGetSize(jac,&n,NULL);
306:   MatGetOwnershipRange(jac,&Istart,&Iend);
307:   if (Istart==0) FirstBlock=PETSC_TRUE;
308:   if (Iend==n) LastBlock=PETSC_TRUE;
309:   h = user->h;
310:   c = user->kappa/(lambda-user->kappa);

312:   /*
313:      Interior grid points
314:   */
315:   for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
316:     j[0] = i-1; j[1] = i; j[2] = i+1;
317:     A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
318:     MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
319:   }

321:   /*
322:      Boundary points
323:   */
324:   if (FirstBlock) {
325:     i = 0;
326:     j[0] = 0; j[1] = 1;
327:     A[0] = -2.0*h/3.0; A[1] = -h/6.0;
328:     MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
329:   }

331:   if (LastBlock) {
332:     i = n-1;
333:     j[0] = n-2; j[1] = n-1;
334:     A[0] = -h/6.0; A[1] = -h/3.0-c*c;
335:     MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
336:   }

338:   /*
339:      Assemble matrix
340:   */
341:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
342:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
343:   PetscFunctionReturn(0);
344: }

346: /* ------------------------------------------------------------------- */
347: /*
348:    CheckSolution - Given a computed solution (lambda,x) check if it
349:    satisfies the analytic solution.

351:    Input Parameters:
352: +  lambda - the computed eigenvalue
353: -  y      - the computed eigenvector

355:    Output Parameter:
356: .  error - norm of difference between the computed and exact eigenvector
357: */
358: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
359: {
360:   PetscScalar    *uu;
361:   PetscInt       i,n,Istart,Iend;
362:   PetscReal      nu,x;
363:   Vec            u;
364:   ApplicationCtx *user = (ApplicationCtx*)ctx;

367:   nu = PetscSqrtReal(PetscRealPart(lambda));
368:   VecDuplicate(y,&u);
369:   VecGetSize(u,&n);
370:   VecGetOwnershipRange(y,&Istart,&Iend);
371:   VecGetArray(u,&uu);
372:   for (i=Istart;i<Iend;i++) {
373:     x = (i+1)*user->h;
374:     uu[i-Istart] = PetscSinReal(nu*x);
375:   }
376:   VecRestoreArray(u,&uu);
377:   VecNormalize(u,NULL);
378:   VecAXPY(u,-1.0,y);
379:   VecNorm(u,NORM_2,error);
380:   VecDestroy(&u);
381:   PetscFunctionReturn(0);
382: }

384: /* ------------------------------------------------------------------- */
385: /*
386:    FixSign - Force the eigenfunction to be real and positive, since
387:    some eigensolvers may return the eigenvector multiplied by a
388:    complex number of modulus one.

390:    Input/Output Parameter:
391: .  x - the computed vector
392: */
393: PetscErrorCode FixSign(Vec x)
394: {
395:   PetscMPIInt       rank;
396:   PetscScalar       sign=0.0;
397:   const PetscScalar *xx;

400:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
401:   if (!rank) {
402:     VecGetArrayRead(x,&xx);
403:     sign = *xx/PetscAbsScalar(*xx);
404:     VecRestoreArrayRead(x,&xx);
405:   }
406:   MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
407:   VecScale(x,1.0/sign);
408:   PetscFunctionReturn(0);
409: }

411: /*TEST

413:    build:
414:       requires: !single

416:    test:
417:       suffix: 1
418:       args: -nep_target 4
419:       filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
420:       requires: !single

422:    testset:
423:       args: -nep_type nleigs -nep_target 10 -nep_nev 4
424:       test:
425:          suffix: 2
426:          filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
427:          args: -rg_interval_endpoints 4,200
428:          requires: !single !complex
429:       test:
430:          suffix: 2_complex
431:          filter: sed -e "s/[+-]0.[0-9]*i//" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
432:          output_file: output/ex20_2.out
433:          args: -rg_interval_endpoints 4,200,-.1,.1 -nep_target_real
434:          requires: !single complex

436:    testset:
437:       args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_ncv 18 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
438:       test:
439:          suffix: 3
440:          filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
441:          output_file: output/ex20_2.out
442:          args: -rg_interval_endpoints 4,200
443:          requires: !single !complex
444:       test:
445:          suffix: 3_complex
446:          filter: sed -e "s/[+-]0.[0-9]*i//" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
447:          output_file: output/ex20_2.out
448:          args: -rg_interval_endpoints 4,200,-.1,.1
449:          requires: !single complex

451:    test:
452:       suffix: 4
453:       args: -n 256 -nep_nev 2 -nep_target 10
454:       filter: sed -e "s/[+-]0\.0*i//g" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
455:       requires: !single

457: TEST*/