Actual source code: ex21.c
slepc-3.6.3 2016-03-29
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: }