Actual source code: ex20.c

slepc-main 2024-11-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[] = "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,NULL,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));

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

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

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

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

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

104:   PetscCall(NEPSetTolerances(nep,1e-9,PETSC_CURRENT));
105:   PetscCall(NEPSetDimensions(nep,1,PETSC_DETERMINE,PETSC_DETERMINE));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

410: /*TEST

412:    build:
413:       requires: !single

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

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

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

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

456: TEST*/