Actual source code: ex21.c

slepc-3.6.3 2016-03-29
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.

  8:    SLEPc is free software: you can redistribute it and/or modify it under  the
  9:    terms of version 3 of the GNU Lesser General Public License as published by
 10:    the Free Software Foundation.

 12:    SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
 13:    WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
 14:    FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
 15:    more details.

 17:    You  should have received a copy of the GNU Lesser General  Public  License
 18:    along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
 19:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 20: */

 22: static char help[] = "Simple 1-D nonlinear eigenproblem (matrix-free version, sequential).\n\n"
 23:   "The command line options are:\n"
 24:   "  -n <n>, where <n> = number of grid subdivisions\n\n";

 26: /*
 27:    Solve 1-D PDE
 28:             -u'' = lambda*u
 29:    on [0,1] subject to
 30:             u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
 31: */

 33: #include <slepcnep.h>

 35: /*
 36:    User-defined routines
 37: */
 38: PetscErrorCode FormInitialGuess(Vec);
 39: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
 40: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);

 42: /*
 43:    Matrix operations and context
 44: */
 45: PetscErrorCode MatMult_Fun(Mat,Vec,Vec);
 46: PetscErrorCode MatGetDiagonal_Fun(Mat,Vec);
 47: PetscErrorCode MatDestroy_Fun(Mat);
 48: PetscErrorCode MatDuplicate_Fun(Mat,MatDuplicateOption,Mat*);
 49: PetscErrorCode MatMult_Jac(Mat,Vec,Vec);
 50: PetscErrorCode MatDestroy_Jac(Mat);

 52: typedef struct {
 53:   PetscScalar lambda,kappa;
 54:   PetscReal   h;
 55: } MatCtx;

 57: /*
 58:    User-defined application context
 59: */
 60: typedef struct {
 61:   PetscScalar kappa;   /* ratio between stiffness of spring and attached mass */
 62:   PetscReal   h;       /* mesh spacing */
 63: } ApplicationCtx;

 67: int main(int argc,char **argv)
 68: {
 69:   NEP            nep;             /* nonlinear eigensolver context */
 70:   Mat            F,J;             /* Function and Jacobian matrices */
 71:   ApplicationCtx ctx;             /* user-defined context */
 72:   MatCtx         *ctxF,*ctxJ;     /* contexts for shell matrices */
 73:   NEPType        type;
 74:   PetscInt       n=128,nev,its;
 75:   KSP            ksp;
 76:   PC             pc;
 77:   PetscMPIInt    size;
 78:   PetscBool      terse;

 81:   SlepcInitialize(&argc,&argv,(char*)0,help);
 82:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 83:   if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
 84:   PetscOptionsGetInt(NULL,"-n",&n,NULL);
 85:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%D\n\n",n);
 86:   ctx.h = 1.0/(PetscReal)n;
 87:   ctx.kappa = 1.0;

 89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90:      Create nonlinear eigensolver context
 91:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 93:   NEPCreate(PETSC_COMM_WORLD,&nep);

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:      Create matrix data structure; set Function evaluation routine
 97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 99:   PetscNew(&ctxF);
100:   ctxF->h = ctx.h;
101:   ctxF->kappa = ctx.kappa;

103:   MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxF,&F);
104:   MatShellSetOperation(F,MATOP_MULT,(void(*)())MatMult_Fun);
105:   MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
106:   MatShellSetOperation(F,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
107:   MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);

109:   /*
110:      Set Function matrix data structure and default Function evaluation
111:      routine
112:   */
113:   NEPSetFunction(nep,F,F,FormFunction,NULL);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Create matrix data structure; set Jacobian evaluation routine
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   PetscNew(&ctxJ);
120:   ctxJ->h = ctx.h;
121:   ctxJ->kappa = ctx.kappa;

123:   MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxJ,&J);
124:   MatShellSetOperation(J,MATOP_MULT,(void(*)())MatMult_Jac);
125:   MatShellSetOperation(J,MATOP_DESTROY,(void(*)())MatDestroy_Jac);

127:   /*
128:      Set Jacobian matrix data structure and default Jacobian evaluation
129:      routine
130:   */
131:   NEPSetJacobian(nep,J,FormJacobian,NULL);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:      Customize nonlinear solver; set runtime options
135:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

137:   NEPGetKSP(nep,&ksp);
138:   KSPSetType(ksp,KSPBCGS);
139:   KSPGetPC(ksp,&pc);
140:   PCSetType(pc,PCJACOBI);

142:   /*
143:      Set solver parameters at runtime
144:   */
145:   NEPSetFromOptions(nep);

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:                       Solve the eigensystem
149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

151:   NEPSolve(nep);
152:   NEPGetIterationNumber(nep,&its);
153:   PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %D\n\n",its);

155:   /*
156:      Optional: Get some information from the solver and display it
157:   */
158:   NEPGetType(nep,&type);
159:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
160:   NEPGetDimensions(nep,&nev,NULL,NULL);
161:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:                     Display solution and clean up
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

167:   /* show detailed info unless -terse option is given by user */
168:   PetscOptionsHasName(NULL,"-terse",&terse);
169:   if (terse) {
170:     NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
171:   } else {
172:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
173:     NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
174:     NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
175:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
176:   }
177:   NEPDestroy(&nep);
178:   MatDestroy(&F);
179:   MatDestroy(&J);
180:   SlepcFinalize();
181:   return 0;
182: }

184: /* ------------------------------------------------------------------- */
187: /*
188:    FormInitialGuess - Computes initial guess.

190:    Input/Output Parameter:
191: .  x - the solution vector
192: */
193: PetscErrorCode FormInitialGuess(Vec x)
194: {

198:   VecSet(x,1.0);
199:   return(0);
200: }

202: /* ------------------------------------------------------------------- */
205: /*
206:    FormFunction - Computes Function matrix  T(lambda)

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

213:    Output Parameters:
214: .  fun - Function matrix
215: .  B   - optionally different preconditioning matrix
216: */
217: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
218: {
220:   MatCtx         *ctxF;

223:   MatShellGetContext(fun,(void**)&ctxF);
224:   ctxF->lambda = lambda;
225:   return(0);
226: }

228: /* ------------------------------------------------------------------- */
231: /*
232:    FormJacobian - Computes Jacobian matrix  T'(lambda)

234:    Input Parameters:
235: .  nep    - the NEP context
236: .  lambda - real part of the scalar argument
237: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

239:    Output Parameters:
240: .  jac - Jacobian matrix
241: .  B   - optionally different preconditioning matrix
242: */
243: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
244: {
246:   MatCtx         *ctxJ;

249:   MatShellGetContext(jac,(void**)&ctxJ);
250:   ctxJ->lambda = lambda;
251:   return(0);
252: }

254: /* ------------------------------------------------------------------- */
257: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
258: {
259:   PetscErrorCode    ierr;
260:   MatCtx            *ctx;
261:   PetscInt          i,n;
262:   const PetscScalar *px;
263:   PetscScalar       *py,c,d,de,oe;
264:   PetscReal         h;

267:   MatShellGetContext(A,(void**)&ctx);
268:   VecGetArrayRead(x,&px);
269:   VecGetArray(y,&py);

271:   VecGetSize(x,&n);
272:   h = ctx->h;
273:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
274:   d = n;
275:   de = 2.0*(d-ctx->lambda*h/3.0);   /* diagonal entry */
276:   oe = -d-ctx->lambda*h/6.0;        /* offdiagonal entry */
277:   py[0] = de*px[0] + oe*px[1];
278:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
279:   de = d-ctx->lambda*h/3.0+ctx->lambda*c;   /* diagonal entry of last row */
280:   py[n-1] = oe*px[n-2] + de*px[n-1];

282:   VecRestoreArrayRead(x,&px);
283:   VecRestoreArray(y,&py);
284:   return(0);
285: }

287: /* ------------------------------------------------------------------- */
290: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
291: {
292:   PetscErrorCode    ierr;
293:   MatCtx            *ctx;
294:   PetscInt          n;
295:   PetscScalar       *pd,c,d;
296:   PetscReal         h;

299:   MatShellGetContext(A,(void**)&ctx);
300:   VecGetSize(diag,&n);
301:   h = ctx->h;
302:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
303:   d = n;
304:   VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
305:   VecGetArray(diag,&pd);
306:   pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
307:   VecRestoreArray(diag,&pd);
308:   return(0);
309: }

311: /* ------------------------------------------------------------------- */
314: PetscErrorCode MatDestroy_Fun(Mat A)
315: {
316:   MatCtx         *ctx;

320:   MatShellGetContext(A,(void**)&ctx);
321:   PetscFree(ctx);
322:   return(0);
323: }

325: /* ------------------------------------------------------------------- */
328: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
329: {
330:   MatCtx         *actx,*bctx;
331:   PetscInt       n;
332:   MPI_Comm       comm;

336:   MatShellGetContext(A,(void**)&actx);
337:   MatGetSize(A,&n,NULL);

339:   PetscNew(&bctx);
340:   bctx->h      = actx->h;
341:   bctx->kappa  = actx->kappa;
342:   bctx->lambda = actx->lambda;

344:   PetscObjectGetComm((PetscObject)A,&comm);
345:   MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
346:   MatShellSetOperation(*B,MATOP_MULT,(void(*)())MatMult_Fun);
347:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
348:   MatShellSetOperation(*B,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
349:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);
350:   return(0);
351: }

353: /* ------------------------------------------------------------------- */
356: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
357: {
358:   PetscErrorCode    ierr;
359:   MatCtx            *ctx;
360:   PetscInt          i,n;
361:   const PetscScalar *px;
362:   PetscScalar       *py,c,de,oe;
363:   PetscReal         h;

366:   MatShellGetContext(A,(void**)&ctx);
367:   VecGetArrayRead(x,&px);
368:   VecGetArray(y,&py);

370:   VecGetSize(x,&n);
371:   h = ctx->h;
372:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
373:   de = -2.0*h/3.0;    /* diagonal entry */
374:   oe = -h/6.0;        /* offdiagonal entry */
375:   py[0] = de*px[0] + oe*px[1];
376:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
377:   de = -h/3.0-c*c;    /* diagonal entry of last row */
378:   py[n-1] = oe*px[n-2] + de*px[n-1];

380:   VecRestoreArrayRead(x,&px);
381:   VecRestoreArray(y,&py);
382:   return(0);
383: }

385: /* ------------------------------------------------------------------- */
388: PetscErrorCode MatDestroy_Jac(Mat A)
389: {
390:   MatCtx         *ctx;

394:   MatShellGetContext(A,(void**)&ctx);
395:   PetscFree(ctx);
396:   return(0);
397: }