Actual source code: ex21.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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 (matrix-free version, sequential).\n\n"
 12:   "The command line options are:\n"
 13:   "  -n <n>, where <n> = number of grid subdivisions\n\n";

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

 22: #include <slepcnep.h>

 24: /*
 25:    User-defined routines
 26: */
 27: PetscErrorCode FormInitialGuess(Vec);
 28: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
 29: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);

 31: /*
 32:    Matrix operations and context
 33: */
 34: PetscErrorCode MatMult_Fun(Mat,Vec,Vec);
 35: PetscErrorCode MatGetDiagonal_Fun(Mat,Vec);
 36: PetscErrorCode MatDestroy_Fun(Mat);
 37: PetscErrorCode MatDuplicate_Fun(Mat,MatDuplicateOption,Mat*);
 38: PetscErrorCode MatMult_Jac(Mat,Vec,Vec);
 39: PetscErrorCode MatGetDiagonal_Jac(Mat,Vec);
 40: PetscErrorCode MatDestroy_Jac(Mat);

 42: typedef struct {
 43:   PetscScalar lambda,kappa;
 44:   PetscReal   h;
 45: } MatCtx;

 47: /*
 48:    User-defined application context
 49: */
 50: typedef struct {
 51:   PetscScalar kappa;   /* ratio between stiffness of spring and attached mass */
 52:   PetscReal   h;       /* mesh spacing */
 53: } ApplicationCtx;

 55: int main(int argc,char **argv)
 56: {
 57:   NEP            nep;             /* nonlinear eigensolver context */
 58:   Mat            F,J;             /* Function and Jacobian matrices */
 59:   ApplicationCtx ctx;             /* user-defined context */
 60:   MatCtx         *ctxF,*ctxJ;     /* contexts for shell matrices */
 61:   PetscInt       n=128,nev;
 62:   KSP            ksp;
 63:   PC             pc;
 64:   PetscMPIInt    size;
 65:   PetscBool      terse;

 68:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 69:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 70:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 71:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 72:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%D\n\n",n);
 73:   ctx.h = 1.0/(PetscReal)n;
 74:   ctx.kappa = 1.0;

 76:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 77:      Create nonlinear eigensolver context
 78:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 80:   NEPCreate(PETSC_COMM_WORLD,&nep);

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:      Create matrix data structure; set Function evaluation routine
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 86:   PetscNew(&ctxF);
 87:   ctxF->h = ctx.h;
 88:   ctxF->kappa = ctx.kappa;

 90:   MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxF,&F);
 91:   MatShellSetOperation(F,MATOP_MULT,(void(*)(void))MatMult_Fun);
 92:   MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
 93:   MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
 94:   MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);

 96:   /*
 97:      Set Function matrix data structure and default Function evaluation
 98:      routine
 99:   */
100:   NEPSetFunction(nep,F,F,FormFunction,NULL);

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:      Create matrix data structure; set Jacobian evaluation routine
104:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   PetscNew(&ctxJ);
107:   ctxJ->h = ctx.h;
108:   ctxJ->kappa = ctx.kappa;

110:   MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxJ,&J);
111:   MatShellSetOperation(J,MATOP_MULT,(void(*)(void))MatMult_Jac);
112:   MatShellSetOperation(J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Jac);
113:   MatShellSetOperation(J,MATOP_DESTROY,(void(*)(void))MatDestroy_Jac);

115:   /*
116:      Set Jacobian matrix data structure and default Jacobian evaluation
117:      routine
118:   */
119:   NEPSetJacobian(nep,J,FormJacobian,NULL);

121:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122:      Customize nonlinear solver; set runtime options
123:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

125:   NEPSetType(nep,NEPRII);
126:   NEPRIISetLagPreconditioner(nep,0);
127:   NEPRIIGetKSP(nep,&ksp);
128:   KSPSetType(ksp,KSPBCGS);
129:   KSPGetPC(ksp,&pc);
130:   PCSetType(pc,PCJACOBI);

132:   /*
133:      Set solver parameters at runtime
134:   */
135:   NEPSetFromOptions(nep);

137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:                       Solve the eigensystem
139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

141:   NEPSolve(nep);
142:   NEPGetDimensions(nep,&nev,NULL,NULL);
143:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:                     Display solution and clean up
147:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

149:   /* show detailed info unless -terse option is given by user */
150:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
151:   if (terse) {
152:     NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
153:   } else {
154:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
155:     NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
156:     NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
157:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
158:   }
159:   NEPDestroy(&nep);
160:   MatDestroy(&F);
161:   MatDestroy(&J);
162:   SlepcFinalize();
163:   return ierr;
164: }

166: /* ------------------------------------------------------------------- */
167: /*
168:    FormInitialGuess - Computes initial guess.

170:    Input/Output Parameter:
171: .  x - the solution vector
172: */
173: PetscErrorCode FormInitialGuess(Vec x)
174: {

178:   VecSet(x,1.0);
179:   return(0);
180: }

182: /* ------------------------------------------------------------------- */
183: /*
184:    FormFunction - Computes Function matrix  T(lambda)

186:    Input Parameters:
187: .  nep    - the NEP context
188: .  lambda - real part of the scalar argument
189: .  ctx    - optional user-defined context, as set by NEPSetFunction()

191:    Output Parameters:
192: .  fun - Function matrix
193: .  B   - optionally different preconditioning matrix
194: */
195: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
196: {
198:   MatCtx         *ctxF;

201:   MatShellGetContext(fun,(void**)&ctxF);
202:   ctxF->lambda = lambda;
203:   return(0);
204: }

206: /* ------------------------------------------------------------------- */
207: /*
208:    FormJacobian - Computes Jacobian matrix  T'(lambda)

210:    Input Parameters:
211: .  nep    - the NEP context
212: .  lambda - real part of the scalar argument
213: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

215:    Output Parameters:
216: .  jac - Jacobian matrix
217: */
218: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
219: {
221:   MatCtx         *ctxJ;

224:   MatShellGetContext(jac,(void**)&ctxJ);
225:   ctxJ->lambda = lambda;
226:   return(0);
227: }

229: /* ------------------------------------------------------------------- */
230: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
231: {
232:   PetscErrorCode    ierr;
233:   MatCtx            *ctx;
234:   PetscInt          i,n;
235:   const PetscScalar *px;
236:   PetscScalar       *py,c,d,de,oe;
237:   PetscReal         h;

240:   MatShellGetContext(A,(void**)&ctx);
241:   VecGetArrayRead(x,&px);
242:   VecGetArray(y,&py);

244:   VecGetSize(x,&n);
245:   h = ctx->h;
246:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
247:   d = n;
248:   de = 2.0*(d-ctx->lambda*h/3.0);   /* diagonal entry */
249:   oe = -d-ctx->lambda*h/6.0;        /* offdiagonal entry */
250:   py[0] = de*px[0] + oe*px[1];
251:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
252:   de = d-ctx->lambda*h/3.0+ctx->lambda*c;   /* diagonal entry of last row */
253:   py[n-1] = oe*px[n-2] + de*px[n-1];

255:   VecRestoreArrayRead(x,&px);
256:   VecRestoreArray(y,&py);
257:   return(0);
258: }

260: /* ------------------------------------------------------------------- */
261: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
262: {
264:   MatCtx         *ctx;
265:   PetscInt       n;
266:   PetscScalar    *pd,c,d;
267:   PetscReal      h;

270:   MatShellGetContext(A,(void**)&ctx);
271:   VecGetSize(diag,&n);
272:   h = ctx->h;
273:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
274:   d = n;
275:   VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
276:   VecGetArray(diag,&pd);
277:   pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
278:   VecRestoreArray(diag,&pd);
279:   return(0);
280: }

282: /* ------------------------------------------------------------------- */
283: PetscErrorCode MatDestroy_Fun(Mat A)
284: {
285:   MatCtx         *ctx;

289:   MatShellGetContext(A,(void**)&ctx);
290:   PetscFree(ctx);
291:   return(0);
292: }

294: /* ------------------------------------------------------------------- */
295: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
296: {
297:   MatCtx         *actx,*bctx;
298:   PetscInt       n;
299:   MPI_Comm       comm;

303:   MatShellGetContext(A,(void**)&actx);
304:   MatGetSize(A,&n,NULL);

306:   PetscNew(&bctx);
307:   bctx->h      = actx->h;
308:   bctx->kappa  = actx->kappa;
309:   bctx->lambda = actx->lambda;

311:   PetscObjectGetComm((PetscObject)A,&comm);
312:   MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
313:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_Fun);
314:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
315:   MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
316:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
317:   return(0);
318: }

320: /* ------------------------------------------------------------------- */
321: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
322: {
323:   PetscErrorCode    ierr;
324:   MatCtx            *ctx;
325:   PetscInt          i,n;
326:   const PetscScalar *px;
327:   PetscScalar       *py,c,de,oe;
328:   PetscReal         h;

331:   MatShellGetContext(A,(void**)&ctx);
332:   VecGetArrayRead(x,&px);
333:   VecGetArray(y,&py);

335:   VecGetSize(x,&n);
336:   h = ctx->h;
337:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
338:   de = -2.0*h/3.0;    /* diagonal entry */
339:   oe = -h/6.0;        /* offdiagonal entry */
340:   py[0] = de*px[0] + oe*px[1];
341:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
342:   de = -h/3.0-c*c;    /* diagonal entry of last row */
343:   py[n-1] = oe*px[n-2] + de*px[n-1];

345:   VecRestoreArrayRead(x,&px);
346:   VecRestoreArray(y,&py);
347:   return(0);
348: }

350: /* ------------------------------------------------------------------- */
351: PetscErrorCode MatGetDiagonal_Jac(Mat A,Vec diag)
352: {
354:   MatCtx         *ctx;
355:   PetscInt       n;
356:   PetscScalar    *pd,c;
357:   PetscReal      h;

360:   MatShellGetContext(A,(void**)&ctx);
361:   VecGetSize(diag,&n);
362:   h = ctx->h;
363:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
364:   VecSet(diag,-2.0*h/3.0);
365:   VecGetArray(diag,&pd);
366:   pd[n-1] = -h/3.0-c*c;
367:   VecRestoreArray(diag,&pd);
368:   return(0);
369: }

371: /* ------------------------------------------------------------------- */
372: PetscErrorCode MatDestroy_Jac(Mat A)
373: {
374:   MatCtx         *ctx;

378:   MatShellGetContext(A,(void**)&ctx);
379:   PetscFree(ctx);
380:   return(0);
381: }

383: /*TEST

385:    testset:
386:       args: -terse
387:       requires: !single
388:       output_file: output/ex21_1.out
389:       test:
390:          suffix: 1_rii
391:          args: -nep_type rii -nep_target 4
392:       test:
393:          suffix: 1_slp
394:          args: -nep_type slp -nep_slp_pc_type jacobi -nep_slp_ksp_type bcgs -nep_target 10

396: TEST*/