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: }