Actual source code: ex21.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, 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;
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,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: NEPSetType(nep,NEPRII);
138: NEPRIISetLagPreconditioner(nep,0);
139: NEPRIIGetKSP(nep,&ksp);
140: KSPSetType(ksp,KSPBCGS);
141: KSPGetPC(ksp,&pc);
142: PCSetType(pc,PCJACOBI);
144: /*
145: Set solver parameters at runtime
146: */
147: NEPSetFromOptions(nep);
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Solve the eigensystem
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: NEPSolve(nep);
154: NEPGetType(nep,&type);
155: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
156: NEPGetDimensions(nep,&nev,NULL,NULL);
157: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
159: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160: Display solution and clean up
161: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: /* show detailed info unless -terse option is given by user */
164: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
165: if (terse) {
166: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
167: } else {
168: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
169: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
170: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
171: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
172: }
173: NEPDestroy(&nep);
174: MatDestroy(&F);
175: MatDestroy(&J);
176: SlepcFinalize();
177: return ierr;
178: }
180: /* ------------------------------------------------------------------- */
183: /*
184: FormInitialGuess - Computes initial guess.
186: Input/Output Parameter:
187: . x - the solution vector
188: */
189: PetscErrorCode FormInitialGuess(Vec x)
190: {
194: VecSet(x,1.0);
195: return(0);
196: }
198: /* ------------------------------------------------------------------- */
201: /*
202: FormFunction - Computes Function matrix T(lambda)
204: Input Parameters:
205: . nep - the NEP context
206: . lambda - real part of the scalar argument
207: . ctx - optional user-defined context, as set by NEPSetFunction()
209: Output Parameters:
210: . fun - Function matrix
211: . B - optionally different preconditioning matrix
212: */
213: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
214: {
216: MatCtx *ctxF;
219: MatShellGetContext(fun,(void**)&ctxF);
220: ctxF->lambda = lambda;
221: return(0);
222: }
224: /* ------------------------------------------------------------------- */
227: /*
228: FormJacobian - Computes Jacobian matrix T'(lambda)
230: Input Parameters:
231: . nep - the NEP context
232: . lambda - real part of the scalar argument
233: . ctx - optional user-defined context, as set by NEPSetJacobian()
235: Output Parameters:
236: . jac - Jacobian matrix
237: . B - optionally different preconditioning matrix
238: */
239: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
240: {
242: MatCtx *ctxJ;
245: MatShellGetContext(jac,(void**)&ctxJ);
246: ctxJ->lambda = lambda;
247: return(0);
248: }
250: /* ------------------------------------------------------------------- */
253: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
254: {
255: PetscErrorCode ierr;
256: MatCtx *ctx;
257: PetscInt i,n;
258: const PetscScalar *px;
259: PetscScalar *py,c,d,de,oe;
260: PetscReal h;
263: MatShellGetContext(A,(void**)&ctx);
264: VecGetArrayRead(x,&px);
265: VecGetArray(y,&py);
267: VecGetSize(x,&n);
268: h = ctx->h;
269: c = ctx->kappa/(ctx->lambda-ctx->kappa);
270: d = n;
271: de = 2.0*(d-ctx->lambda*h/3.0); /* diagonal entry */
272: oe = -d-ctx->lambda*h/6.0; /* offdiagonal entry */
273: py[0] = de*px[0] + oe*px[1];
274: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
275: de = d-ctx->lambda*h/3.0+ctx->lambda*c; /* diagonal entry of last row */
276: py[n-1] = oe*px[n-2] + de*px[n-1];
278: VecRestoreArrayRead(x,&px);
279: VecRestoreArray(y,&py);
280: return(0);
281: }
283: /* ------------------------------------------------------------------- */
286: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
287: {
288: PetscErrorCode ierr;
289: MatCtx *ctx;
290: PetscInt n;
291: PetscScalar *pd,c,d;
292: PetscReal h;
295: MatShellGetContext(A,(void**)&ctx);
296: VecGetSize(diag,&n);
297: h = ctx->h;
298: c = ctx->kappa/(ctx->lambda-ctx->kappa);
299: d = n;
300: VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
301: VecGetArray(diag,&pd);
302: pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
303: VecRestoreArray(diag,&pd);
304: return(0);
305: }
307: /* ------------------------------------------------------------------- */
310: PetscErrorCode MatDestroy_Fun(Mat A)
311: {
312: MatCtx *ctx;
316: MatShellGetContext(A,(void**)&ctx);
317: PetscFree(ctx);
318: return(0);
319: }
321: /* ------------------------------------------------------------------- */
324: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
325: {
326: MatCtx *actx,*bctx;
327: PetscInt n;
328: MPI_Comm comm;
332: MatShellGetContext(A,(void**)&actx);
333: MatGetSize(A,&n,NULL);
335: PetscNew(&bctx);
336: bctx->h = actx->h;
337: bctx->kappa = actx->kappa;
338: bctx->lambda = actx->lambda;
340: PetscObjectGetComm((PetscObject)A,&comm);
341: MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
342: MatShellSetOperation(*B,MATOP_MULT,(void(*)())MatMult_Fun);
343: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
344: MatShellSetOperation(*B,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
345: MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);
346: return(0);
347: }
349: /* ------------------------------------------------------------------- */
352: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
353: {
354: PetscErrorCode ierr;
355: MatCtx *ctx;
356: PetscInt i,n;
357: const PetscScalar *px;
358: PetscScalar *py,c,de,oe;
359: PetscReal h;
362: MatShellGetContext(A,(void**)&ctx);
363: VecGetArrayRead(x,&px);
364: VecGetArray(y,&py);
366: VecGetSize(x,&n);
367: h = ctx->h;
368: c = ctx->kappa/(ctx->lambda-ctx->kappa);
369: de = -2.0*h/3.0; /* diagonal entry */
370: oe = -h/6.0; /* offdiagonal entry */
371: py[0] = de*px[0] + oe*px[1];
372: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
373: de = -h/3.0-c*c; /* diagonal entry of last row */
374: py[n-1] = oe*px[n-2] + de*px[n-1];
376: VecRestoreArrayRead(x,&px);
377: VecRestoreArray(y,&py);
378: return(0);
379: }
381: /* ------------------------------------------------------------------- */
384: PetscErrorCode MatDestroy_Jac(Mat A)
385: {
386: MatCtx *ctx;
390: MatShellGetContext(A,(void**)&ctx);
391: PetscFree(ctx);
392: return(0);
393: }