Actual source code: test4.c
slepc-3.16.0 2021-09-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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[] = "Test the RII solver with a user-provided KSP.\n\n"
12: "This is a simplified version of ex20.\n"
13: "The command line options are:\n"
14: " -n <n>, where <n> = number of grid subdivisions.\n";
16: /*
17: Solve 1-D PDE
18: -u'' = lambda*u
19: on [0,1] subject to
20: u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
21: */
23: #include <slepcnep.h>
25: /*
26: User-defined routines
27: */
28: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
29: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
31: /*
32: User-defined application context
33: */
34: typedef struct {
35: PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
36: PetscReal h; /* mesh spacing */
37: } ApplicationCtx;
39: int main(int argc,char **argv)
40: {
41: NEP nep;
42: KSP ksp;
43: PC pc;
44: Mat F,J;
45: ApplicationCtx ctx;
46: PetscInt n=128,lag,its;
47: PetscBool terse,flg,cct,herm;
48: PetscReal thres;
51: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
52: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
53: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%D\n\n",n);
54: ctx.h = 1.0/(PetscReal)n;
55: ctx.kappa = 1.0;
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create a standalone KSP with appropriate settings
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: KSPCreate(PETSC_COMM_WORLD,&ksp);
62: KSPSetType(ksp,KSPBCGS);
63: KSPGetPC(ksp,&pc);
64: PCSetType(pc,PCBJACOBI);
65: KSPSetFromOptions(ksp);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Prepare nonlinear eigensolver context
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: NEPCreate(PETSC_COMM_WORLD,&nep);
73: /* Create Function and Jacobian matrices; set evaluation routines */
74: MatCreate(PETSC_COMM_WORLD,&F);
75: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
76: MatSetFromOptions(F);
77: MatSeqAIJSetPreallocation(F,3,NULL);
78: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
79: MatSetUp(F);
80: NEPSetFunction(nep,F,F,FormFunction,&ctx);
82: MatCreate(PETSC_COMM_WORLD,&J);
83: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
84: MatSetFromOptions(J);
85: MatSeqAIJSetPreallocation(J,3,NULL);
86: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
87: MatSetUp(J);
88: NEPSetJacobian(nep,J,FormJacobian,&ctx);
90: NEPSetType(nep,NEPRII);
91: NEPRIISetKSP(nep,ksp);
92: NEPRIISetMaximumIterations(nep,6);
93: NEPRIISetConstCorrectionTol(nep,PETSC_TRUE);
94: NEPRIISetHermitian(nep,PETSC_TRUE);
95: NEPSetFromOptions(nep);
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Solve the eigensystem
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: NEPSolve(nep);
102: PetscObjectTypeCompare((PetscObject)nep,NEPRII,&flg);
103: if (flg) {
104: NEPRIIGetMaximumIterations(nep,&its);
105: NEPRIIGetLagPreconditioner(nep,&lag);
106: NEPRIIGetDeflationThreshold(nep,&thres);
107: NEPRIIGetConstCorrectionTol(nep,&cct);
108: NEPRIIGetHermitian(nep,&herm);
109: PetscPrintf(PETSC_COMM_WORLD," Maximum inner iterations of RII is %D\n",its);
110: PetscPrintf(PETSC_COMM_WORLD," Preconditioner rebuilt every %D iterations\n",lag);
111: if (thres>0.0) { PetscPrintf(PETSC_COMM_WORLD," Using deflation threshold=%g\n",(double)thres); }
112: if (cct) { PetscPrintf(PETSC_COMM_WORLD," Using a constant correction tolerance\n"); }
113: if (herm) { PetscPrintf(PETSC_COMM_WORLD," Hermitian version of scalar equation\n"); }
114: PetscPrintf(PETSC_COMM_WORLD,"\n");
115: }
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Display solution and clean up
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: /* show detailed info unless -terse option is given by user */
122: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
123: if (terse) {
124: NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL);
125: } else {
126: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
127: NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD);
128: NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
129: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
130: }
132: NEPDestroy(&nep);
133: KSPDestroy(&ksp);
134: MatDestroy(&F);
135: MatDestroy(&J);
136: SlepcFinalize();
137: return ierr;
138: }
140: /* ------------------------------------------------------------------- */
141: /*
142: FormFunction - Computes Function matrix T(lambda)
144: Input Parameters:
145: . nep - the NEP context
146: . lambda - the scalar argument
147: . ctx - optional user-defined context, as set by NEPSetFunction()
149: Output Parameters:
150: . fun - Function matrix
151: . B - optionally different preconditioning matrix
152: */
153: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
154: {
156: ApplicationCtx *user = (ApplicationCtx*)ctx;
157: PetscScalar A[3],c,d;
158: PetscReal h;
159: PetscInt i,n,j[3],Istart,Iend;
160: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
163: /*
164: Compute Function entries and insert into matrix
165: */
166: MatGetSize(fun,&n,NULL);
167: MatGetOwnershipRange(fun,&Istart,&Iend);
168: if (Istart==0) FirstBlock=PETSC_TRUE;
169: if (Iend==n) LastBlock=PETSC_TRUE;
170: h = user->h;
171: c = user->kappa/(lambda-user->kappa);
172: d = n;
174: /*
175: Interior grid points
176: */
177: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
178: j[0] = i-1; j[1] = i; j[2] = i+1;
179: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
180: MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
181: }
183: /*
184: Boundary points
185: */
186: if (FirstBlock) {
187: i = 0;
188: j[0] = 0; j[1] = 1;
189: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
190: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
191: }
193: if (LastBlock) {
194: i = n-1;
195: j[0] = n-2; j[1] = n-1;
196: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
197: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
198: }
200: /*
201: Assemble matrix
202: */
203: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
204: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
205: if (fun != B) {
206: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
207: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
208: }
209: return(0);
210: }
212: /* ------------------------------------------------------------------- */
213: /*
214: FormJacobian - Computes Jacobian matrix T'(lambda)
216: Input Parameters:
217: . nep - the NEP context
218: . lambda - the scalar argument
219: . ctx - optional user-defined context, as set by NEPSetJacobian()
221: Output Parameters:
222: . jac - Jacobian matrix
223: . B - optionally different preconditioning matrix
224: */
225: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
226: {
228: ApplicationCtx *user = (ApplicationCtx*)ctx;
229: PetscScalar A[3],c;
230: PetscReal h;
231: PetscInt i,n,j[3],Istart,Iend;
232: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
235: /*
236: Compute Jacobian entries and insert into matrix
237: */
238: MatGetSize(jac,&n,NULL);
239: MatGetOwnershipRange(jac,&Istart,&Iend);
240: if (Istart==0) FirstBlock=PETSC_TRUE;
241: if (Iend==n) LastBlock=PETSC_TRUE;
242: h = user->h;
243: c = user->kappa/(lambda-user->kappa);
245: /*
246: Interior grid points
247: */
248: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
249: j[0] = i-1; j[1] = i; j[2] = i+1;
250: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
251: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
252: }
254: /*
255: Boundary points
256: */
257: if (FirstBlock) {
258: i = 0;
259: j[0] = 0; j[1] = 1;
260: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
261: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
262: }
264: if (LastBlock) {
265: i = n-1;
266: j[0] = n-2; j[1] = n-1;
267: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
268: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
269: }
271: /*
272: Assemble matrix
273: */
274: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
275: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
276: return(0);
277: }
279: /*TEST
281: test:
282: suffix: 1
283: args: -nep_target 21 -nep_rii_lag_preconditioner 2 -terse
284: requires: !single
286: TEST*/