Actual source code: ex21.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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(*)(void))MatMult_Fun);
92: MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
93: MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
94: MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(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(*)(void))MatMult_Jac);
112: MatShellSetOperation(J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Jac);
113: MatShellSetOperation(J,MATOP_DESTROY,(void(*)(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: */
218: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
219: {
221: MatCtx *ctxJ;
224: MatShellGetContext(jac,(void**)&ctxJ);
225: ctxJ->lambda = lambda;
226: return(0);
227: }
229: /* ------------------------------------------------------------------- */
230: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
231: {
232: PetscErrorCode ierr;
233: MatCtx *ctx;
234: PetscInt i,n;
235: const PetscScalar *px;
236: PetscScalar *py,c,d,de,oe;
237: PetscReal h;
240: MatShellGetContext(A,(void**)&ctx);
241: VecGetArrayRead(x,&px);
242: VecGetArray(y,&py);
244: VecGetSize(x,&n);
245: h = ctx->h;
246: c = ctx->kappa/(ctx->lambda-ctx->kappa);
247: d = n;
248: de = 2.0*(d-ctx->lambda*h/3.0); /* diagonal entry */
249: oe = -d-ctx->lambda*h/6.0; /* offdiagonal entry */
250: py[0] = de*px[0] + oe*px[1];
251: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
252: de = d-ctx->lambda*h/3.0+ctx->lambda*c; /* diagonal entry of last row */
253: py[n-1] = oe*px[n-2] + de*px[n-1];
255: VecRestoreArrayRead(x,&px);
256: VecRestoreArray(y,&py);
257: return(0);
258: }
260: /* ------------------------------------------------------------------- */
261: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
262: {
264: MatCtx *ctx;
265: PetscInt n;
266: PetscScalar *pd,c,d;
267: PetscReal h;
270: MatShellGetContext(A,(void**)&ctx);
271: VecGetSize(diag,&n);
272: h = ctx->h;
273: c = ctx->kappa/(ctx->lambda-ctx->kappa);
274: d = n;
275: VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
276: VecGetArray(diag,&pd);
277: pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
278: VecRestoreArray(diag,&pd);
279: return(0);
280: }
282: /* ------------------------------------------------------------------- */
283: PetscErrorCode MatDestroy_Fun(Mat A)
284: {
285: MatCtx *ctx;
289: MatShellGetContext(A,(void**)&ctx);
290: PetscFree(ctx);
291: return(0);
292: }
294: /* ------------------------------------------------------------------- */
295: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
296: {
297: MatCtx *actx,*bctx;
298: PetscInt n;
299: MPI_Comm comm;
303: MatShellGetContext(A,(void**)&actx);
304: MatGetSize(A,&n,NULL);
306: PetscNew(&bctx);
307: bctx->h = actx->h;
308: bctx->kappa = actx->kappa;
309: bctx->lambda = actx->lambda;
311: PetscObjectGetComm((PetscObject)A,&comm);
312: MatCreateShell(comm,n,n,n,n,(void*)bctx,B);
313: MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_Fun);
314: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
315: MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
316: MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
317: return(0);
318: }
320: /* ------------------------------------------------------------------- */
321: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
322: {
323: PetscErrorCode ierr;
324: MatCtx *ctx;
325: PetscInt i,n;
326: const PetscScalar *px;
327: PetscScalar *py,c,de,oe;
328: PetscReal h;
331: MatShellGetContext(A,(void**)&ctx);
332: VecGetArrayRead(x,&px);
333: VecGetArray(y,&py);
335: VecGetSize(x,&n);
336: h = ctx->h;
337: c = ctx->kappa/(ctx->lambda-ctx->kappa);
338: de = -2.0*h/3.0; /* diagonal entry */
339: oe = -h/6.0; /* offdiagonal entry */
340: py[0] = de*px[0] + oe*px[1];
341: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
342: de = -h/3.0-c*c; /* diagonal entry of last row */
343: py[n-1] = oe*px[n-2] + de*px[n-1];
345: VecRestoreArrayRead(x,&px);
346: VecRestoreArray(y,&py);
347: return(0);
348: }
350: /* ------------------------------------------------------------------- */
351: PetscErrorCode MatGetDiagonal_Jac(Mat A,Vec diag)
352: {
354: MatCtx *ctx;
355: PetscInt n;
356: PetscScalar *pd,c;
357: PetscReal h;
360: MatShellGetContext(A,(void**)&ctx);
361: VecGetSize(diag,&n);
362: h = ctx->h;
363: c = ctx->kappa/(ctx->lambda-ctx->kappa);
364: VecSet(diag,-2.0*h/3.0);
365: VecGetArray(diag,&pd);
366: pd[n-1] = -h/3.0-c*c;
367: VecRestoreArray(diag,&pd);
368: return(0);
369: }
371: /* ------------------------------------------------------------------- */
372: PetscErrorCode MatDestroy_Jac(Mat A)
373: {
374: MatCtx *ctx;
378: MatShellGetContext(A,(void**)&ctx);
379: PetscFree(ctx);
380: return(0);
381: }
383: /*TEST
385: testset:
386: args: -terse
387: requires: !single
388: output_file: output/ex21_1.out
389: test:
390: suffix: 1_rii
391: args: -nep_type rii -nep_target 4
392: test:
393: suffix: 1_slp
394: args: -nep_type slp -nep_slp_pc_type jacobi -nep_slp_ksp_type bcgs -nep_target 10
396: TEST*/