Actual source code: ex21.c
slepc-3.18.2 2023-01-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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).\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 FormFunction(NEP,PetscScalar,Mat,Mat,void*);
28: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
30: /*
31: Matrix operations and context
32: */
33: PetscErrorCode MatMult_Fun(Mat,Vec,Vec);
34: PetscErrorCode MatGetDiagonal_Fun(Mat,Vec);
35: PetscErrorCode MatDestroy_Fun(Mat);
36: PetscErrorCode MatDuplicate_Fun(Mat,MatDuplicateOption,Mat*);
37: PetscErrorCode MatMult_Jac(Mat,Vec,Vec);
38: PetscErrorCode MatGetDiagonal_Jac(Mat,Vec);
39: PetscErrorCode MatDestroy_Jac(Mat);
41: typedef struct {
42: PetscScalar lambda,kappa;
43: PetscReal h;
44: PetscMPIInt next,prev;
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 rank,size;
65: PetscBool terse;
68: SlepcInitialize(&argc,&argv,(char*)0,help);
69: MPI_Comm_size(PETSC_COMM_WORLD,&size);
70: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
71: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
72: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\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;
89: ctxF->next = rank==size-1? MPI_PROC_NULL: rank+1;
90: ctxF->prev = rank==0? MPI_PROC_NULL: rank-1;
92: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctxF,&F);
93: MatShellSetOperation(F,MATOP_MULT,(void(*)(void))MatMult_Fun);
94: MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
95: MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
96: MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
98: /*
99: Set Function matrix data structure and default Function evaluation
100: routine
101: */
102: NEPSetFunction(nep,F,F,FormFunction,NULL);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Create matrix data structure; set Jacobian evaluation routine
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: PetscNew(&ctxJ);
109: ctxJ->h = ctx.h;
110: ctxJ->kappa = ctx.kappa;
111: ctxJ->next = rank==size-1? MPI_PROC_NULL: rank+1;
112: ctxJ->prev = rank==0? MPI_PROC_NULL: rank-1;
114: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctxJ,&J);
115: MatShellSetOperation(J,MATOP_MULT,(void(*)(void))MatMult_Jac);
116: MatShellSetOperation(J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Jac);
117: MatShellSetOperation(J,MATOP_DESTROY,(void(*)(void))MatDestroy_Jac);
119: /*
120: Set Jacobian matrix data structure and default Jacobian evaluation
121: routine
122: */
123: NEPSetJacobian(nep,J,FormJacobian,NULL);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Customize nonlinear solver; set runtime options
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: NEPSetType(nep,NEPRII);
130: NEPRIISetLagPreconditioner(nep,0);
131: NEPRIIGetKSP(nep,&ksp);
132: KSPSetType(ksp,KSPBCGS);
133: KSPGetPC(ksp,&pc);
134: PCSetType(pc,PCJACOBI);
136: /*
137: Set solver parameters at runtime
138: */
139: NEPSetFromOptions(nep);
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Solve the eigensystem
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: NEPSolve(nep);
146: NEPGetDimensions(nep,&nev,NULL,NULL);
147: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Display solution and clean up
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: /* show detailed info unless -terse option is given by user */
154: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
155: if (terse) NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
156: else {
157: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
158: NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
159: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
160: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
161: }
162: NEPDestroy(&nep);
163: MatDestroy(&F);
164: MatDestroy(&J);
165: SlepcFinalize();
166: return 0;
167: }
169: /* ------------------------------------------------------------------- */
170: /*
171: FormFunction - Computes Function matrix T(lambda)
173: Input Parameters:
174: . nep - the NEP context
175: . lambda - real part of the scalar argument
176: . ctx - optional user-defined context, as set by NEPSetFunction()
178: Output Parameters:
179: . fun - Function matrix
180: . B - optionally different preconditioning matrix
181: */
182: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
183: {
184: MatCtx *ctxF;
187: MatShellGetContext(fun,&ctxF);
188: ctxF->lambda = lambda;
189: return 0;
190: }
192: /* ------------------------------------------------------------------- */
193: /*
194: FormJacobian - Computes Jacobian matrix T'(lambda)
196: Input Parameters:
197: . nep - the NEP context
198: . lambda - real part of the scalar argument
199: . ctx - optional user-defined context, as set by NEPSetJacobian()
201: Output Parameters:
202: . jac - Jacobian matrix
203: */
204: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
205: {
206: MatCtx *ctxJ;
209: MatShellGetContext(jac,&ctxJ);
210: ctxJ->lambda = lambda;
211: return 0;
212: }
214: /* ------------------------------------------------------------------- */
215: PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
216: {
217: MatCtx *ctx;
218: PetscInt i,n,N;
219: const PetscScalar *px;
220: PetscScalar *py,c,d,de,oe,upper=0.0,lower=0.0;
221: PetscReal h;
222: MPI_Comm comm;
225: MatShellGetContext(A,&ctx);
226: VecGetArrayRead(x,&px);
227: VecGetArray(y,&py);
228: VecGetSize(x,&N);
229: VecGetLocalSize(x,&n);
231: PetscObjectGetComm((PetscObject)A,&comm);
232: MPI_Sendrecv(px,1,MPIU_SCALAR,ctx->prev,0,&lower,1,MPIU_SCALAR,ctx->next,0,comm,MPI_STATUS_IGNORE);
233: MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,ctx->next,0,&upper,1,MPIU_SCALAR,ctx->prev,0,comm,MPI_STATUS_IGNORE);
235: h = ctx->h;
236: c = ctx->kappa/(ctx->lambda-ctx->kappa);
237: d = N;
238: de = 2.0*(d-ctx->lambda*h/3.0); /* diagonal entry */
239: oe = -d-ctx->lambda*h/6.0; /* offdiagonal entry */
240: py[0] = oe*upper + de*px[0] + oe*px[1];
241: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
242: if (ctx->next==MPI_PROC_NULL) de = d-ctx->lambda*h/3.0+ctx->lambda*c; /* diagonal entry of last row */
243: py[n-1] = oe*px[n-2] + de*px[n-1] + oe*lower;
245: VecRestoreArrayRead(x,&px);
246: VecRestoreArray(y,&py);
247: return 0;
248: }
250: /* ------------------------------------------------------------------- */
251: PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
252: {
253: MatCtx *ctx;
254: PetscInt n,N;
255: PetscScalar *pd,c,d;
256: PetscReal h;
259: MatShellGetContext(A,&ctx);
260: VecGetSize(diag,&N);
261: VecGetLocalSize(diag,&n);
262: h = ctx->h;
263: c = ctx->kappa/(ctx->lambda-ctx->kappa);
264: d = N;
265: VecSet(diag,2.0*(d-ctx->lambda*h/3.0));
266: VecGetArray(diag,&pd);
267: pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
268: VecRestoreArray(diag,&pd);
269: return 0;
270: }
272: /* ------------------------------------------------------------------- */
273: PetscErrorCode MatDestroy_Fun(Mat A)
274: {
275: MatCtx *ctx;
277: MatShellGetContext(A,&ctx);
278: PetscFree(ctx);
279: return 0;
280: }
282: /* ------------------------------------------------------------------- */
283: PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
284: {
285: MatCtx *actx,*bctx;
286: PetscInt m,n,M,N;
287: MPI_Comm comm;
289: MatShellGetContext(A,&actx);
290: MatGetSize(A,&M,&N);
291: MatGetLocalSize(A,&m,&n);
293: PetscNew(&bctx);
294: bctx->h = actx->h;
295: bctx->kappa = actx->kappa;
296: bctx->lambda = actx->lambda;
297: bctx->next = actx->next;
298: bctx->prev = actx->prev;
300: PetscObjectGetComm((PetscObject)A,&comm);
301: MatCreateShell(comm,m,n,M,N,(void*)bctx,B);
302: MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_Fun);
303: MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun);
304: MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun);
305: MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun);
306: return 0;
307: }
309: /* ------------------------------------------------------------------- */
310: PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
311: {
312: MatCtx *ctx;
313: PetscInt i,n;
314: const PetscScalar *px;
315: PetscScalar *py,c,de,oe,upper=0.0,lower=0.0;
316: PetscReal h;
317: MPI_Comm comm;
320: MatShellGetContext(A,&ctx);
321: VecGetArrayRead(x,&px);
322: VecGetArray(y,&py);
323: VecGetLocalSize(x,&n);
325: PetscObjectGetComm((PetscObject)A,&comm);
326: MPI_Sendrecv(px,1,MPIU_SCALAR,ctx->prev,0,&lower,1,MPIU_SCALAR,ctx->next,0,comm,MPI_STATUS_IGNORE);
327: MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,ctx->next,0,&upper,1,MPIU_SCALAR,ctx->prev,0,comm,MPI_STATUS_IGNORE);
329: h = ctx->h;
330: c = ctx->kappa/(ctx->lambda-ctx->kappa);
331: de = -2.0*h/3.0; /* diagonal entry */
332: oe = -h/6.0; /* offdiagonal entry */
333: py[0] = oe*upper + de*px[0] + oe*px[1];
334: for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
335: if (ctx->next==MPI_PROC_NULL) de = -h/3.0-c*c; /* diagonal entry of last row */
336: py[n-1] = oe*px[n-2] + de*px[n-1] + oe*lower;
338: VecRestoreArrayRead(x,&px);
339: VecRestoreArray(y,&py);
340: return 0;
341: }
343: /* ------------------------------------------------------------------- */
344: PetscErrorCode MatGetDiagonal_Jac(Mat A,Vec diag)
345: {
346: MatCtx *ctx;
347: PetscInt n;
348: PetscScalar *pd,c;
349: PetscReal h;
352: MatShellGetContext(A,&ctx);
353: VecGetLocalSize(diag,&n);
354: h = ctx->h;
355: c = ctx->kappa/(ctx->lambda-ctx->kappa);
356: VecSet(diag,-2.0*h/3.0);
357: VecGetArray(diag,&pd);
358: pd[n-1] = -h/3.0-c*c;
359: VecRestoreArray(diag,&pd);
360: return 0;
361: }
363: /* ------------------------------------------------------------------- */
364: PetscErrorCode MatDestroy_Jac(Mat A)
365: {
366: MatCtx *ctx;
368: MatShellGetContext(A,&ctx);
369: PetscFree(ctx);
370: return 0;
371: }
373: /*TEST
375: testset:
376: nsize: {{1 2}}
377: args: -terse
378: requires: !single
379: output_file: output/ex21_1.out
380: filter: sed -e "s/[+-]0\.0*i//g" -e "s/+0i//g"
381: test:
382: suffix: 1_rii
383: args: -nep_type rii -nep_target 4
384: test:
385: suffix: 1_slp
386: args: -nep_type slp -nep_slp_pc_type jacobi -nep_slp_ksp_type bcgs -nep_target 10
388: TEST*/