Actual source code: ex20.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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;

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

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

 66:   NEPCreate(PETSC_COMM_WORLD,&nep);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

192:   NEPDestroy(&nep);
193:   MatDestroy(&F);
194:   MatDestroy(&J);
195:   VecDestroy(&x);
196:   SlepcFinalize();
197:   return ierr;
198: }

200: /* ------------------------------------------------------------------- */
201: /*
202:    FormInitialGuess - Computes initial guess.

204:    Input/Output Parameter:
205: .  x - the solution vector
206: */
207: PetscErrorCode FormInitialGuess(Vec x)
208: {

212:   VecSet(x,1.0);
213:   return(0);
214: }

216: /* ------------------------------------------------------------------- */
217: /*
218:    FormFunction - Computes Function matrix  T(lambda)

220:    Input Parameters:
221: .  nep    - the NEP context
222: .  lambda - the scalar argument
223: .  ctx    - optional user-defined context, as set by NEPSetFunction()

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

239:   /*
240:      Compute Function entries and insert into matrix
241:   */
242:   MatGetSize(fun,&n,NULL);
243:   MatGetOwnershipRange(fun,&Istart,&Iend);
244:   if (Istart==0) FirstBlock=PETSC_TRUE;
245:   if (Iend==n) LastBlock=PETSC_TRUE;
246:   h = user->h;
247:   c = user->kappa/(lambda-user->kappa);
248:   d = n;

250:   /*
251:      Interior grid points
252:   */
253:   for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
254:     j[0] = i-1; j[1] = i; j[2] = i+1;
255:     A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
256:     MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
257:   }

259:   /*
260:      Boundary points
261:   */
262:   if (FirstBlock) {
263:     i = 0;
264:     j[0] = 0; j[1] = 1;
265:     A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
266:     MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
267:   }

269:   if (LastBlock) {
270:     i = n-1;
271:     j[0] = n-2; j[1] = n-1;
272:     A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
273:     MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
274:   }

276:   /*
277:      Assemble matrix
278:   */
279:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
280:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
281:   if (fun != B) {
282:     MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
283:     MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
284:   }
285:   return(0);
286: }

288: /* ------------------------------------------------------------------- */
289: /*
290:    FormJacobian - Computes Jacobian matrix  T'(lambda)

292:    Input Parameters:
293: .  nep    - the NEP context
294: .  lambda - the scalar argument
295: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

297:    Output Parameters:
298: .  jac - Jacobian matrix
299: */
300: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
301: {
303:   ApplicationCtx *user = (ApplicationCtx*)ctx;
304:   PetscScalar    A[3],c;
305:   PetscReal      h;
306:   PetscInt       i,n,j[3],Istart,Iend;
307:   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;

310:   /*
311:      Compute Jacobian entries and insert into matrix
312:   */
313:   MatGetSize(jac,&n,NULL);
314:   MatGetOwnershipRange(jac,&Istart,&Iend);
315:   if (Istart==0) FirstBlock=PETSC_TRUE;
316:   if (Iend==n) LastBlock=PETSC_TRUE;
317:   h = user->h;
318:   c = user->kappa/(lambda-user->kappa);

320:   /*
321:      Interior grid points
322:   */
323:   for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
324:     j[0] = i-1; j[1] = i; j[2] = i+1;
325:     A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
326:     MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
327:   }

329:   /*
330:      Boundary points
331:   */
332:   if (FirstBlock) {
333:     i = 0;
334:     j[0] = 0; j[1] = 1;
335:     A[0] = -2.0*h/3.0; A[1] = -h/6.0;
336:     MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
337:   }

339:   if (LastBlock) {
340:     i = n-1;
341:     j[0] = n-2; j[1] = n-1;
342:     A[0] = -h/6.0; A[1] = -h/3.0-c*c;
343:     MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
344:   }

346:   /*
347:      Assemble matrix
348:   */
349:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
350:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
351:   return(0);
352: }

354: /* ------------------------------------------------------------------- */
355: /*
356:    CheckSolution - Given a computed solution (lambda,x) check if it
357:    satisfies the analytic solution.

359:    Input Parameters:
360: +  lambda - the computed eigenvalue
361: -  y      - the computed eigenvector

363:    Output Parameter:
364: .  error - norm of difference between the computed and exact eigenvector
365: */
366: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
367: {
369:   PetscScalar    nu,*uu;
370:   PetscInt       i,n,Istart,Iend;
371:   PetscReal      x;
372:   Vec            u;
373:   ApplicationCtx *user = (ApplicationCtx*)ctx;

376:   nu = PetscSqrtScalar(lambda);
377:   VecDuplicate(y,&u);
378:   VecGetSize(u,&n);
379:   VecGetOwnershipRange(y,&Istart,&Iend);
380:   VecGetArray(u,&uu);
381:   for (i=Istart;i<Iend;i++) {
382:     x = (i+1)*user->h;
383:     uu[i-Istart] = PetscSinReal(nu*x);
384:   }
385:   VecRestoreArray(u,&uu);
386:   VecNormalize(u,NULL);
387:   VecAXPY(u,-1.0,y);
388:   VecNorm(u,NORM_2,error);
389:   VecDestroy(&u);
390:   return(0);
391: }

393: /* ------------------------------------------------------------------- */
394: /*
395:    FixSign - Force the eigenfunction to be real and positive, since
396:    some eigensolvers may return the eigenvector multiplied by a
397:    complex number of modulus one.

399:    Input/Output Parameter:
400: .  x - the computed vector
401: */
402: PetscErrorCode FixSign(Vec x)
403: {
404:   PetscErrorCode    ierr;
405:   PetscMPIInt       rank;
406:   PetscScalar       sign=0.0;
407:   const PetscScalar *xx;

410:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
411:   if (!rank) {
412:     VecGetArrayRead(x,&xx);
413:     sign = *xx/PetscAbsScalar(*xx);
414:     VecRestoreArrayRead(x,&xx);
415:   }
416:   MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
417:   VecScale(x,1.0/sign);
418:   return(0);
419: }

421: /*TEST

423:    build:
424:       requires: !single

426:    test:
427:       suffix: 1
428:       args: -nep_target 4
429:       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 = /"
430:       requires: !single

432:    testset:
433:       args: -nep_type nleigs -nep_target 10 -nep_nev 4
434:       test:
435:          suffix: 2
436:          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 = /"
437:          args: -rg_interval_endpoints 4,200
438:          requires: !single !complex
439:       test:
440:          suffix: 2_complex
441:          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 = /"
442:          output_file: output/ex20_2.out
443:          args: -rg_interval_endpoints 4,200,-.1,.1 -nep_target_real
444:          requires: !single complex

446:    testset:
447:       args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_ncv 18 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
448:       test:
449:          suffix: 3
450:          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 = /"
451:          output_file: output/ex20_2.out
452:          args: -rg_interval_endpoints 4,200
453:          requires: !single !complex
454:       test:
455:          suffix: 3_complex
456:          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 = /"
457:          output_file: output/ex20_2.out
458:          args: -rg_interval_endpoints 4,200,-.1,.1
459:          requires: !single complex

461:    test:
462:       suffix: 4
463:       args: -n 256 -nep_nev 2 -nep_target 10
464:       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 = /"
465:       requires: !single

467: TEST*/