Actual source code: ex20.c

slepc-3.19.1 2023-05-15
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:   PetscFunctionBeginUser;
 55:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
 56:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
 57:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n));
 58:   ctx.h = 1.0/(PetscReal)n;
 59:   ctx.kappa = 1.0;
 60:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-draw_sol",&draw_sol,NULL));

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

 66:   PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

153:   if (nconv>0) {
154:     /*
155:        Display eigenvalues and relative errors
156:     */
157:     PetscCall(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:       PetscCall(NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL));
165:       PetscCall(FixSign(x));
166:       /*
167:          Compute residual norm and error
168:       */
169:       PetscCall(NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm));
170:       PetscCall(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) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g     %12g\n",(double)re,(double)im,(double)norm,(double)error));
180:       else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"   %12f         %12g     %12g\n",(double)re,(double)norm,(double)error));
181:       if (draw_sol) {
182:         PetscCall(PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1));
183:         PetscCall(VecView(x,PETSC_VIEWER_DRAW_WORLD));
184:       }
185:     }
186:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
187:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

412: /*TEST

414:    build:
415:       requires: !single

417:    test:
418:       suffix: 1
419:       args: -nep_target 4
420:       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 = /"
421:       requires: !single

423:    testset:
424:       args: -nep_type nleigs -nep_target 10 -nep_nev 4
425:       test:
426:          suffix: 2
427:          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 = /"
428:          args: -rg_interval_endpoints 4,200
429:          requires: !single !complex
430:       test:
431:          suffix: 2_complex
432:          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 = /"
433:          output_file: output/ex20_2.out
434:          args: -rg_interval_endpoints 4,200,-.1,.1 -nep_target_real
435:          requires: !single complex

437:    testset:
438:       args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_ncv 18 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
439:       test:
440:          suffix: 3
441:          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 = /"
442:          output_file: output/ex20_2.out
443:          args: -rg_interval_endpoints 4,200
444:          requires: !single !complex
445:       test:
446:          suffix: 3_complex
447:          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 = /"
448:          output_file: output/ex20_2.out
449:          args: -rg_interval_endpoints 4,200,-.1,.1
450:          requires: !single complex

452:    test:
453:       suffix: 4
454:       args: -n 256 -nep_nev 2 -nep_target 10
455:       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 = /"
456:       requires: !single

458: TEST*/