Actual source code: ex21.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2013, 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*,MatStructure*,void*);
 40: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat*,MatStructure*,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:   PetscScalar    lambda;          /* eigenvalue */
 71:   Mat            F,J;             /* Function and Jacobian matrices */
 72:   ApplicationCtx ctx;             /* user-defined context */
 73:   MatCtx         *ctxF,*ctxJ;     /* contexts for shell matrices */
 74:   NEPType        type;
 75:   PetscInt       n=128,nev,i,its,nconv;
 76:   KSP            ksp;
 77:   PC             pc;
 78:   PetscMPIInt    size;
 79:   PetscReal      re,im,norm;

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

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

 94:   NEPCreate(PETSC_COMM_WORLD,&nep);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

168:   /*
169:      Get number of converged approximate eigenpairs
170:   */
171:   NEPGetConverged(nep,&nconv);
172:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);

174:   if (nconv>0) {
175:     /*
176:        Display eigenvalues and relative errors
177:     */
178:     PetscPrintf(PETSC_COMM_WORLD,
179:          "           k              ||T(k)x||\n"
180:          "   ----------------- ------------------\n");
181:     for (i=0;i<nconv;i++) {
182:       /*
183:         Get converged eigenpairs (in this example they are always real)
184:       */
185:       NEPGetEigenpair(nep,i,&lambda,NULL);
186:       /*
187:          Compute residual norm
188:       */
189:       NEPComputeRelativeError(nep,i,&norm);

191: #if defined(PETSC_USE_COMPLEX)
192:       re = PetscRealPart(lambda);
193:       im = PetscImaginaryPart(lambda);
194: #else
195:       re = lambda;
196:       im = 0.0;
197: #endif
198:       if (im!=0.0) {
199:         PetscPrintf(PETSC_COMM_WORLD," %9F%+9F j %12G\n",re,im,norm);
200:       } else {
201:         PetscPrintf(PETSC_COMM_WORLD,"   %12F         %12G\n",re,norm);
202:       }
203:     }
204:     PetscPrintf(PETSC_COMM_WORLD,"\n");
205:   }

207:   NEPDestroy(&nep);
208:   MatDestroy(&F);
209:   MatDestroy(&J);
210:   SlepcFinalize();
211:   return 0;
212: }

214: /* ------------------------------------------------------------------- */
217: /*
218:    FormInitialGuess - Computes initial guess.

220:    Input/Output Parameter:
221: .  x - the solution vector
222: */
223: PetscErrorCode FormInitialGuess(Vec x)
224: {

228:   VecSet(x,1.0);
229:   return(0);
230: }

232: /* ------------------------------------------------------------------- */
235: /*
236:    FormFunction - Computes Function matrix  T(lambda)

238:    Input Parameters:
239: .  nep    - the NEP context
240: .  lambda - real part of the scalar argument
241: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

243:    Output Parameters:
244: .  fun - Function matrix
245: .  B   - optionally different preconditioning matrix
246: .  flg - flag indicating matrix structure
247: */
248: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat *fun,Mat *B,MatStructure *flg,void *ctx)
249: {
251:   MatCtx         *ctxF;

254:   MatShellGetContext(*fun,(void**)&ctxF);
255:   ctxF->lambda = lambda;
256:   return(0);
257: }

259: /* ------------------------------------------------------------------- */
262: /*
263:    FormJacobian - Computes Jacobian matrix  T'(lambda)

265:    Input Parameters:
266: .  nep    - the NEP context
267: .  lambda - real part of the scalar argument
268: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

270:    Output Parameters:
271: .  jac - Jacobian matrix
272: .  B   - optionally different preconditioning matrix
273: .  flg - flag indicating matrix structure
274: */
275: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat *jac,MatStructure *flg,void *ctx)
276: {
278:   MatCtx         *ctxJ;

281:   MatShellGetContext(*jac,(void**)&ctxJ);
282:   ctxJ->lambda = lambda;
283:   return(0);
284: }

286: /* ------------------------------------------------------------------- */
289: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
290: {
291:   PetscErrorCode    ierr;
292:   MatCtx            *ctx;
293:   PetscInt          i,n;
294:   const PetscScalar *px;
295:   PetscScalar       *py,c,d,de,oe;
296:   PetscReal         h;

299:   MatShellGetContext(A,(void**)&ctx);
300:   VecGetArrayRead(x,&px);
301:   VecGetArray(y,&py);

303:   VecGetSize(x,&n);
304:   h = ctx->h;
305:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
306:   d = n;
307:   de = 2.0*(d-ctx->lambda*h/3.0);   /* diagonal entry */
308:   oe = -d-ctx->lambda*h/6.0;        /* offdiagonal entry */
309:   py[0] = de*px[0] + oe*px[1];
310:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
311:   de = d-ctx->lambda*h/3.0+ctx->lambda*c;   /* diagonal entry of last row */
312:   py[n-1] = oe*px[n-2] + de*px[n-1];

314:   VecRestoreArrayRead(x,&px);
315:   VecRestoreArray(y,&py);
316:   return(0);
317: }

319: /* ------------------------------------------------------------------- */
322: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
323: {
324:   PetscErrorCode    ierr;
325:   MatCtx            *ctx;
326:   PetscInt          n;
327:   PetscScalar       *pd,c,d;
328:   PetscReal         h;

331:   MatShellGetContext(A,(void**)&ctx);
332:   h = ctx->h;
333:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
334:   d = n;
335:   VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
336:   VecGetSize(diag,&n);
337:   VecGetArray(diag,&pd);
338:   pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
339:   VecRestoreArray(diag,&pd);
340:   return(0);
341: }

343: /* ------------------------------------------------------------------- */
346: PetscErrorCode MatDestroy_Fun(Mat A)
347: {
348:   MatCtx         *ctx;

352:   MatShellGetContext(A,(void**)&ctx);
353:   PetscFree(ctx);
354:   return(0);
355: }

357: /* ------------------------------------------------------------------- */
360: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
361: {
362:   MatCtx         *actx,*bctx;
363:   PetscInt       n;
364:   MPI_Comm       comm;

368:   MatShellGetContext(A,(void**)&actx);
369:   MatGetSize(A,&n,NULL);

371:   PetscNew(MatCtx,&bctx);
372:   bctx->h      = actx->h;
373:   bctx->kappa  = actx->kappa;
374:   bctx->lambda = actx->lambda;

376:   PetscObjectGetComm((PetscObject)A,&comm);
377:   MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
378:   MatShellSetOperation(*B,MATOP_MULT,(void(*)())MatMult_Fun);
379:   MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
380:   MatShellSetOperation(*B,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
381:   MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);
382:   return(0);
383: }

385: /* ------------------------------------------------------------------- */
388: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
389: {
390:   PetscErrorCode    ierr;
391:   MatCtx            *ctx;
392:   PetscInt          i,n;
393:   const PetscScalar *px;
394:   PetscScalar       *py,c,de,oe;
395:   PetscReal         h;

398:   MatShellGetContext(A,(void**)&ctx);
399:   VecGetArrayRead(x,&px);
400:   VecGetArray(y,&py);

402:   VecGetSize(x,&n);
403:   h = ctx->h;
404:   c = ctx->kappa/(ctx->lambda-ctx->kappa);
405:   de = -2.0*h/3.0;    /* diagonal entry */
406:   oe = -h/6.0;        /* offdiagonal entry */
407:   py[0] = de*px[0] + oe*px[1];
408:   for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
409:   de = -h/3.0-c*c;    /* diagonal entry of last row */
410:   py[n-1] = oe*px[n-2] + de*px[n-1];

412:   VecRestoreArrayRead(x,&px);
413:   VecRestoreArray(y,&py);
414:   return(0);
415: }

417: /* ------------------------------------------------------------------- */
420: PetscErrorCode MatDestroy_Jac(Mat A)
421: {
422:   MatCtx         *ctx;

426:   MatShellGetContext(A,(void**)&ctx);
427:   PetscFree(ctx);
428:   return(0);
429: }