Actual source code: ex21.c

slepc-3.12.2 2020-01-13
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 FormFunction(NEP,PetscScalar,Mat,Mat,void*);
 28: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);

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

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

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

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

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

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

 79:   NEPCreate(PETSC_COMM_WORLD,&nep);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

165: /* ------------------------------------------------------------------- */
166: /*
167:    FormFunction - Computes Function matrix  T(lambda)

169:    Input Parameters:
170: .  nep    - the NEP context
171: .  lambda - real part of the scalar argument
172: .  ctx    - optional user-defined context, as set by NEPSetFunction()

174:    Output Parameters:
175: .  fun - Function matrix
176: .  B   - optionally different preconditioning matrix
177: */
178: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
179: {
181:   MatCtx         *ctxF;

184:   MatShellGetContext(fun,(void**)&ctxF);
185:   ctxF->lambda = lambda;
186:   return(0);
187: }

189: /* ------------------------------------------------------------------- */
190: /*
191:    FormJacobian - Computes Jacobian matrix  T'(lambda)

193:    Input Parameters:
194: .  nep    - the NEP context
195: .  lambda - real part of the scalar argument
196: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

198:    Output Parameters:
199: .  jac - Jacobian matrix
200: */
201: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
202: {
204:   MatCtx         *ctxJ;

207:   MatShellGetContext(jac,(void**)&ctxJ);
208:   ctxJ->lambda = lambda;
209:   return(0);
210: }

212: /* ------------------------------------------------------------------- */
213: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
214: {
215:   PetscErrorCode    ierr;
216:   MatCtx            *ctx;
217:   PetscInt          i,n;
218:   const PetscScalar *px;
219:   PetscScalar       *py,c,d,de,oe;
220:   PetscReal         h;

223:   MatShellGetContext(A,(void**)&ctx);
224:   VecGetArrayRead(x,&px);
225:   VecGetArray(y,&py);

227:   VecGetSize(x,&n);
228:   h = ctx->h;
229:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
230:   d = n;
231:   de = 2.0*(d-ctx->lambda*h/3.0);   /* diagonal entry */
232:   oe = -d-ctx->lambda*h/6.0;        /* offdiagonal entry */
233:   py[0] = de*px[0] + oe*px[1];
234:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
235:   de = d-ctx->lambda*h/3.0+ctx->lambda*c;   /* diagonal entry of last row */
236:   py[n-1] = oe*px[n-2] + de*px[n-1];

238:   VecRestoreArrayRead(x,&px);
239:   VecRestoreArray(y,&py);
240:   return(0);
241: }

243: /* ------------------------------------------------------------------- */
244: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
245: {
247:   MatCtx         *ctx;
248:   PetscInt       n;
249:   PetscScalar    *pd,c,d;
250:   PetscReal      h;

253:   MatShellGetContext(A,(void**)&ctx);
254:   VecGetSize(diag,&n);
255:   h = ctx->h;
256:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
257:   d = n;
258:   VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
259:   VecGetArray(diag,&pd);
260:   pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
261:   VecRestoreArray(diag,&pd);
262:   return(0);
263: }

265: /* ------------------------------------------------------------------- */
266: PetscErrorCode MatDestroy_Fun(Mat A)
267: {
268:   MatCtx         *ctx;

272:   MatShellGetContext(A,(void**)&ctx);
273:   PetscFree(ctx);
274:   return(0);
275: }

277: /* ------------------------------------------------------------------- */
278: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
279: {
280:   MatCtx         *actx,*bctx;
281:   PetscInt       n;
282:   MPI_Comm       comm;

286:   MatShellGetContext(A,(void**)&actx);
287:   MatGetSize(A,&n,NULL);

289:   PetscNew(&bctx);
290:   bctx->h      = actx->h;
291:   bctx->kappa  = actx->kappa;
292:   bctx->lambda = actx->lambda;

294:   PetscObjectGetComm((PetscObject)A,&comm);
295:   MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
296:   MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_Fun);
297:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
298:   MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
299:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
300:   return(0);
301: }

303: /* ------------------------------------------------------------------- */
304: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
305: {
306:   PetscErrorCode    ierr;
307:   MatCtx            *ctx;
308:   PetscInt          i,n;
309:   const PetscScalar *px;
310:   PetscScalar       *py,c,de,oe;
311:   PetscReal         h;

314:   MatShellGetContext(A,(void**)&ctx);
315:   VecGetArrayRead(x,&px);
316:   VecGetArray(y,&py);

318:   VecGetSize(x,&n);
319:   h = ctx->h;
320:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
321:   de = -2.0*h/3.0;    /* diagonal entry */
322:   oe = -h/6.0;        /* offdiagonal entry */
323:   py[0] = de*px[0] + oe*px[1];
324:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
325:   de = -h/3.0-c*c;    /* diagonal entry of last row */
326:   py[n-1] = oe*px[n-2] + de*px[n-1];

328:   VecRestoreArrayRead(x,&px);
329:   VecRestoreArray(y,&py);
330:   return(0);
331: }

333: /* ------------------------------------------------------------------- */
334: PetscErrorCode MatGetDiagonal_Jac(Mat A,Vec diag)
335: {
337:   MatCtx         *ctx;
338:   PetscInt       n;
339:   PetscScalar    *pd,c;
340:   PetscReal      h;

343:   MatShellGetContext(A,(void**)&ctx);
344:   VecGetSize(diag,&n);
345:   h = ctx->h;
346:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
347:   VecSet(diag,-2.0*h/3.0);
348:   VecGetArray(diag,&pd);
349:   pd[n-1] = -h/3.0-c*c;
350:   VecRestoreArray(diag,&pd);
351:   return(0);
352: }

354: /* ------------------------------------------------------------------- */
355: PetscErrorCode MatDestroy_Jac(Mat A)
356: {
357:   MatCtx         *ctx;

361:   MatShellGetContext(A,(void**)&ctx);
362:   PetscFree(ctx);
363:   return(0);
364: }

366: /*TEST

368:    testset:
369:       args: -terse
370:       requires: !single
371:       output_file: output/ex21_1.out
372:       test:
373:          suffix: 1_rii
374:          args: -nep_type rii -nep_target 4
375:       test:
376:          suffix: 1_slp
377:          args: -nep_type slp -nep_slp_pc_type jacobi -nep_slp_ksp_type bcgs -nep_target 10

379: TEST*/