Actual source code: ex21.c
slepc-3.8.3 2018-04-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2017, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: static char help[] = "Simple 1-D nonlinear eigenproblem (matrix-free version, sequential).\n\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = number of grid subdivisions\n\n";
15: /*
16: Solve 1-D PDE
17: -u'' = lambda*u
18: on [0,1] subject to
19: u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
20: */
22: #include <slepcnep.h>
24: /*
25: User-defined routines
26: */
27: PetscErrorCode FormInitialGuess(Vec);
28: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
29: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
31: /*
32: Matrix operations and context
33: */
34: PetscErrorCode MatMult_Fun(Mat,Vec,Vec);
35: PetscErrorCode MatGetDiagonal_Fun(Mat,Vec);
36: PetscErrorCode MatDestroy_Fun(Mat);
37: PetscErrorCode MatDuplicate_Fun(Mat,MatDuplicateOption,Mat*);
38: PetscErrorCode MatMult_Jac(Mat,Vec,Vec);
39: PetscErrorCode MatGetDiagonal_Jac(Mat,Vec);
40: PetscErrorCode MatDestroy_Jac(Mat);
42: typedef struct {
43: PetscScalar lambda,kappa;
44: PetscReal h;
45: } MatCtx;
47: /*
48: User-defined application context
49: */
50: typedef struct {
51: PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
52: PetscReal h; /* mesh spacing */
53: } ApplicationCtx;
55: int main(int argc,char **argv)
56: {
57: NEP nep; /* nonlinear eigensolver context */
58: Mat F,J; /* Function and Jacobian matrices */
59: ApplicationCtx ctx; /* user-defined context */
60: MatCtx *ctxF,*ctxJ; /* contexts for shell matrices */
61: PetscInt n=128,nev;
62: KSP ksp;
63: PC pc;
64: PetscMPIInt size;
65: PetscBool terse;
68: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
69: MPI_Comm_size(PETSC_COMM_WORLD,&size);
70: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
71: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
72: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%D\n\n",n);
73: ctx.h = 1.0/(PetscReal)n;
74: ctx.kappa = 1.0;
76: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: Create nonlinear eigensolver context
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: NEPCreate(PETSC_COMM_WORLD,&nep);
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Create matrix data structure; set Function evaluation routine
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: PetscNew(&ctxF);
87: ctxF->h = ctx.h;
88: ctxF->kappa = ctx.kappa;
90: MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxF,&F);
91: MatShellSetOperation(F,MATOP_MULT,(void(*)())MatMult_Fun);
92: MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
93: MatShellSetOperation(F,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
94: MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);
96: /*
97: Set Function matrix data structure and default Function evaluation
98: routine
99: */
100: NEPSetFunction(nep,F,F,FormFunction,NULL);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create matrix data structure; set Jacobian evaluation routine
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: PetscNew(&ctxJ);
107: ctxJ->h = ctx.h;
108: ctxJ->kappa = ctx.kappa;
110: MatCreateShell(PETSC_COMM_WORLD,n,n,n,n,(void*)ctxJ,&J);
111: MatShellSetOperation(J,MATOP_MULT,(void(*)())MatMult_Jac);
112: MatShellSetOperation(J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Jac);
113: MatShellSetOperation(J,MATOP_DESTROY,(void(*)())MatDestroy_Jac);
115: /*
116: Set Jacobian matrix data structure and default Jacobian evaluation
117: routine
118: */
119: NEPSetJacobian(nep,J,FormJacobian,NULL);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Customize nonlinear solver; set runtime options
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125: NEPSetType(nep,NEPRII);
126: NEPRIISetLagPreconditioner(nep,0);
127: NEPRIIGetKSP(nep,&ksp);
128: KSPSetType(ksp,KSPBCGS);
129: KSPGetPC(ksp,&pc);
130: PCSetType(pc,PCJACOBI);
132: /*
133: Set solver parameters at runtime
134: */
135: NEPSetFromOptions(nep);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Solve the eigensystem
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: NEPSolve(nep);
142: NEPGetDimensions(nep,&nev,NULL,NULL);
143: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Display solution and clean up
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: /* show detailed info unless -terse option is given by user */
150: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
151: if (terse) {
152: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
153: } else {
154: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
155: NEPReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
156: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
157: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
158: }
159: NEPDestroy(&nep);
160: MatDestroy(&F);
161: MatDestroy(&J);
162: SlepcFinalize();
163: return ierr;
164: }
166: /* ------------------------------------------------------------------- */
167: /*
168: FormInitialGuess - Computes initial guess.
170: Input/Output Parameter:
171: . x - the solution vector
172: */
173: PetscErrorCode FormInitialGuess(Vec x)
174: {
178: VecSet(x,1.0);
179: return(0);
180: }
182: /* ------------------------------------------------------------------- */
183: /*
184: FormFunction - Computes Function matrix T(lambda)
186: Input Parameters:
187: . nep - the NEP context
188: . lambda - real part of the scalar argument
189: . ctx - optional user-defined context, as set by NEPSetFunction()
191: Output Parameters:
192: . fun - Function matrix
193: . B - optionally different preconditioning matrix
194: */
195: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
196: {
198: MatCtx *ctxF;
201: MatShellGetContext(fun,(void**)&ctxF);
202: ctxF->lambda = lambda;
203: return(0);
204: }
206: /* ------------------------------------------------------------------- */
207: /*
208: FormJacobian - Computes Jacobian matrix T'(lambda)
210: Input Parameters:
211: . nep - the NEP context
212: . lambda - real part of the scalar argument
213: . ctx - optional user-defined context, as set by NEPSetJacobian()
215: Output Parameters:
216: . jac - Jacobian matrix
217: . B - optionally different preconditioning matrix
218: */
219: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
220: {
222: MatCtx *ctxJ;
225: MatShellGetContext(jac,(void**)&ctxJ);
226: ctxJ->lambda = lambda;
227: return(0);
228: }
230: /* ------------------------------------------------------------------- */
231: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
232: {
233: PetscErrorCode ierr;
234: MatCtx *ctx;
235: PetscInt i,n;
236: const PetscScalar *px;
237: PetscScalar *py,c,d,de,oe;
238: PetscReal h;
241: MatShellGetContext(A,(void**)&ctx);
242: VecGetArrayRead(x,&px);
243: VecGetArray(y,&py);
245: VecGetSize(x,&n);
246: h = ctx->h;
247: c = ctx->kappa/(ctx->lambda-ctx->kappa);
248: d = n;
249: de = 2.0*(d-ctx->lambda*h/3.0); /* diagonal entry */
250: oe = -d-ctx->lambda*h/6.0; /* offdiagonal entry */
251: py[0] = de*px[0] + oe*px[1];
252: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
253: de = d-ctx->lambda*h/3.0+ctx->lambda*c; /* diagonal entry of last row */
254: py[n-1] = oe*px[n-2] + de*px[n-1];
256: VecRestoreArrayRead(x,&px);
257: VecRestoreArray(y,&py);
258: return(0);
259: }
261: /* ------------------------------------------------------------------- */
262: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
263: {
265: MatCtx *ctx;
266: PetscInt n;
267: PetscScalar *pd,c,d;
268: PetscReal h;
271: MatShellGetContext(A,(void**)&ctx);
272: VecGetSize(diag,&n);
273: h = ctx->h;
274: c = ctx->kappa/(ctx->lambda-ctx->kappa);
275: d = n;
276: VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
277: VecGetArray(diag,&pd);
278: pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
279: VecRestoreArray(diag,&pd);
280: return(0);
281: }
283: /* ------------------------------------------------------------------- */
284: PetscErrorCode MatDestroy_Fun(Mat A)
285: {
286: MatCtx *ctx;
290: MatShellGetContext(A,(void**)&ctx);
291: PetscFree(ctx);
292: return(0);
293: }
295: /* ------------------------------------------------------------------- */
296: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
297: {
298: MatCtx *actx,*bctx;
299: PetscInt n;
300: MPI_Comm comm;
304: MatShellGetContext(A,(void**)&actx);
305: MatGetSize(A,&n,NULL);
307: PetscNew(&bctx);
308: bctx->h = actx->h;
309: bctx->kappa = actx->kappa;
310: bctx->lambda = actx->lambda;
312: PetscObjectGetComm((PetscObject)A,&comm);
313: MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
314: MatShellSetOperation(*B,MATOP_MULT,(void(*)())MatMult_Fun);
315: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_Fun);
316: MatShellSetOperation(*B,MATOP_DESTROY,(void(*)())MatDestroy_Fun);
317: MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)())MatDuplicate_Fun);
318: return(0);
319: }
321: /* ------------------------------------------------------------------- */
322: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
323: {
324: PetscErrorCode ierr;
325: MatCtx *ctx;
326: PetscInt i,n;
327: const PetscScalar *px;
328: PetscScalar *py,c,de,oe;
329: PetscReal h;
332: MatShellGetContext(A,(void**)&ctx);
333: VecGetArrayRead(x,&px);
334: VecGetArray(y,&py);
336: VecGetSize(x,&n);
337: h = ctx->h;
338: c = ctx->kappa/(ctx->lambda-ctx->kappa);
339: de = -2.0*h/3.0; /* diagonal entry */
340: oe = -h/6.0; /* offdiagonal entry */
341: py[0] = de*px[0] + oe*px[1];
342: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
343: de = -h/3.0-c*c; /* diagonal entry of last row */
344: py[n-1] = oe*px[n-2] + de*px[n-1];
346: VecRestoreArrayRead(x,&px);
347: VecRestoreArray(y,&py);
348: return(0);
349: }
351: /* ------------------------------------------------------------------- */
352: PetscErrorCode MatGetDiagonal_Jac(Mat A,Vec diag)
353: {
355: MatCtx *ctx;
356: PetscInt n;
357: PetscScalar *pd,c;
358: PetscReal h;
361: MatShellGetContext(A,(void**)&ctx);
362: VecGetSize(diag,&n);
363: h = ctx->h;
364: c = ctx->kappa/(ctx->lambda-ctx->kappa);
365: VecSet(diag,-2.0*h/3.0);
366: VecGetArray(diag,&pd);
367: pd[n-1] = -h/3.0-c*c;
368: VecRestoreArray(diag,&pd);
369: return(0);
370: }
372: /* ------------------------------------------------------------------- */
373: PetscErrorCode MatDestroy_Jac(Mat A)
374: {
375: MatCtx *ctx;
379: MatShellGetContext(A,(void**)&ctx);
380: PetscFree(ctx);
381: return(0);
382: }