Actual source code: ex20.c
slepc-3.17.1 2022-04-11
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.\n\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = number of grid subdivisions.\n"
14: " -draw_sol, to draw the computed solution.\n\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 FormInitialGuess(Vec);
29: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
30: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
31: PetscErrorCode CheckSolution(PetscScalar,Vec,PetscReal*,void*);
32: PetscErrorCode FixSign(Vec);
34: /*
35: User-defined application context
36: */
37: typedef struct {
38: PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
39: PetscReal h; /* mesh spacing */
40: } ApplicationCtx;
42: int main(int argc,char **argv)
43: {
44: NEP nep; /* nonlinear eigensolver context */
45: Vec x; /* eigenvector */
46: PetscScalar lambda; /* eigenvalue */
47: Mat F,J; /* Function and Jacobian matrices */
48: ApplicationCtx ctx; /* user-defined context */
49: NEPType type;
50: PetscInt n=128,nev,i,its,maxit,nconv;
51: PetscReal re,im,tol,norm,error;
52: PetscBool draw_sol=PETSC_FALSE;
54: SlepcInitialize(&argc,&argv,(char*)0,help);
55: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
56: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n);
57: ctx.h = 1.0/(PetscReal)n;
58: ctx.kappa = 1.0;
59: PetscOptionsGetBool(NULL,NULL,"-draw_sol",&draw_sol,NULL);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create nonlinear eigensolver context
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: NEPCreate(PETSC_COMM_WORLD,&nep);
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create matrix data structure; set Function evaluation routine
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: MatCreate(PETSC_COMM_WORLD,&F);
72: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
73: MatSetFromOptions(F);
74: MatSeqAIJSetPreallocation(F,3,NULL);
75: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
76: MatSetUp(F);
78: /*
79: Set Function matrix data structure and default Function evaluation
80: routine
81: */
82: NEPSetFunction(nep,F,F,FormFunction,&ctx);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create matrix data structure; set Jacobian evaluation routine
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: MatCreate(PETSC_COMM_WORLD,&J);
89: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
90: MatSetFromOptions(J);
91: MatSeqAIJSetPreallocation(J,3,NULL);
92: MatMPIAIJSetPreallocation(J,3,NULL,1,NULL);
93: MatSetUp(J);
95: /*
96: Set Jacobian matrix data structure and default Jacobian evaluation
97: routine
98: */
99: NEPSetJacobian(nep,J,FormJacobian,&ctx);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Customize nonlinear solver; set runtime options
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: NEPSetTolerances(nep,1e-9,PETSC_DEFAULT);
106: NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT);
108: /*
109: Set solver parameters at runtime
110: */
111: NEPSetFromOptions(nep);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Initialize application
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: /*
118: Evaluate initial guess
119: */
120: MatCreateVecs(F,&x,NULL);
121: FormInitialGuess(x);
122: NEPSetInitialSpace(nep,1,&x);
124: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Solve the eigensystem
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: NEPSolve(nep);
129: NEPGetIterationNumber(nep,&its);
130: PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %" PetscInt_FMT "\n\n",its);
132: /*
133: Optional: Get some information from the solver and display it
134: */
135: NEPGetType(nep,&type);
136: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
137: NEPGetDimensions(nep,&nev,NULL,NULL);
138: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
139: NEPGetTolerances(nep,&tol,&maxit);
140: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit);
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Display solution and clean up
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146: /*
147: Get number of converged approximate eigenpairs
148: */
149: NEPGetConverged(nep,&nconv);
150: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %" PetscInt_FMT "\n\n",nconv);
152: if (nconv>0) {
153: /*
154: Display eigenvalues and relative errors
155: */
156: PetscCall(PetscPrintf(PETSC_COMM_WORLD,
157: " k ||T(k)x|| error\n"
158: " ----------------- ------------------ ------------------\n"));
159: for (i=0;i<nconv;i++) {
160: /*
161: Get converged eigenpairs (in this example they are always real)
162: */
163: NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL);
164: FixSign(x);
165: /*
166: Compute residual norm and error
167: */
168: NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm);
169: CheckSolution(lambda,x,&error,&ctx);
171: #if defined(PETSC_USE_COMPLEX)
172: re = PetscRealPart(lambda);
173: im = PetscImaginaryPart(lambda);
174: #else
175: re = lambda;
176: im = 0.0;
177: #endif
178: if (im!=0.0) PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g %12g\n",(double)re,(double)im,(double)norm,(double)error);
179: else PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)error);
180: if (draw_sol) {
181: PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
182: VecView(x,PETSC_VIEWER_DRAW_WORLD);
183: }
184: }
185: PetscPrintf(PETSC_COMM_WORLD,"\n");
186: }
188: NEPDestroy(&nep);
189: MatDestroy(&F);
190: MatDestroy(&J);
191: VecDestroy(&x);
192: SlepcFinalize();
193: return 0;
194: }
196: /* ------------------------------------------------------------------- */
197: /*
198: FormInitialGuess - Computes initial guess.
200: Input/Output Parameter:
201: . x - the solution vector
202: */
203: PetscErrorCode FormInitialGuess(Vec x)
204: {
206: VecSet(x,1.0);
207: PetscFunctionReturn(0);
208: }
210: /* ------------------------------------------------------------------- */
211: /*
212: FormFunction - Computes Function matrix T(lambda)
214: Input Parameters:
215: . nep - the NEP context
216: . lambda - the scalar argument
217: . ctx - optional user-defined context, as set by NEPSetFunction()
219: Output Parameters:
220: . fun - Function matrix
221: . B - optionally different preconditioning matrix
222: */
223: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
224: {
225: ApplicationCtx *user = (ApplicationCtx*)ctx;
226: PetscScalar A[3],c,d;
227: PetscReal h;
228: PetscInt i,n,j[3],Istart,Iend;
229: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
232: /*
233: Compute Function entries and insert into matrix
234: */
235: MatGetSize(fun,&n,NULL);
236: MatGetOwnershipRange(fun,&Istart,&Iend);
237: if (Istart==0) FirstBlock=PETSC_TRUE;
238: if (Iend==n) LastBlock=PETSC_TRUE;
239: h = user->h;
240: c = user->kappa/(lambda-user->kappa);
241: d = n;
243: /*
244: Interior grid points
245: */
246: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
247: j[0] = i-1; j[1] = i; j[2] = i+1;
248: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
249: MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
250: }
252: /*
253: Boundary points
254: */
255: if (FirstBlock) {
256: i = 0;
257: j[0] = 0; j[1] = 1;
258: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
259: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
260: }
262: if (LastBlock) {
263: i = n-1;
264: j[0] = n-2; j[1] = n-1;
265: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
266: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
267: }
269: /*
270: Assemble matrix
271: */
272: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
273: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
274: if (fun != B) {
275: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
276: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
277: }
278: PetscFunctionReturn(0);
279: }
281: /* ------------------------------------------------------------------- */
282: /*
283: FormJacobian - Computes Jacobian matrix T'(lambda)
285: Input Parameters:
286: . nep - the NEP context
287: . lambda - the scalar argument
288: . ctx - optional user-defined context, as set by NEPSetJacobian()
290: Output Parameters:
291: . jac - Jacobian matrix
292: */
293: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
294: {
295: ApplicationCtx *user = (ApplicationCtx*)ctx;
296: PetscScalar A[3],c;
297: PetscReal h;
298: PetscInt i,n,j[3],Istart,Iend;
299: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
302: /*
303: Compute Jacobian entries and insert into matrix
304: */
305: MatGetSize(jac,&n,NULL);
306: MatGetOwnershipRange(jac,&Istart,&Iend);
307: if (Istart==0) FirstBlock=PETSC_TRUE;
308: if (Iend==n) LastBlock=PETSC_TRUE;
309: h = user->h;
310: c = user->kappa/(lambda-user->kappa);
312: /*
313: Interior grid points
314: */
315: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
316: j[0] = i-1; j[1] = i; j[2] = i+1;
317: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
318: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
319: }
321: /*
322: Boundary points
323: */
324: if (FirstBlock) {
325: i = 0;
326: j[0] = 0; j[1] = 1;
327: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
328: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
329: }
331: if (LastBlock) {
332: i = n-1;
333: j[0] = n-2; j[1] = n-1;
334: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
335: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
336: }
338: /*
339: Assemble matrix
340: */
341: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
342: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
343: PetscFunctionReturn(0);
344: }
346: /* ------------------------------------------------------------------- */
347: /*
348: CheckSolution - Given a computed solution (lambda,x) check if it
349: satisfies the analytic solution.
351: Input Parameters:
352: + lambda - the computed eigenvalue
353: - y - the computed eigenvector
355: Output Parameter:
356: . error - norm of difference between the computed and exact eigenvector
357: */
358: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
359: {
360: PetscScalar *uu;
361: PetscInt i,n,Istart,Iend;
362: PetscReal nu,x;
363: Vec u;
364: ApplicationCtx *user = (ApplicationCtx*)ctx;
367: nu = PetscSqrtReal(PetscRealPart(lambda));
368: VecDuplicate(y,&u);
369: VecGetSize(u,&n);
370: VecGetOwnershipRange(y,&Istart,&Iend);
371: VecGetArray(u,&uu);
372: for (i=Istart;i<Iend;i++) {
373: x = (i+1)*user->h;
374: uu[i-Istart] = PetscSinReal(nu*x);
375: }
376: VecRestoreArray(u,&uu);
377: VecNormalize(u,NULL);
378: VecAXPY(u,-1.0,y);
379: VecNorm(u,NORM_2,error);
380: VecDestroy(&u);
381: PetscFunctionReturn(0);
382: }
384: /* ------------------------------------------------------------------- */
385: /*
386: FixSign - Force the eigenfunction to be real and positive, since
387: some eigensolvers may return the eigenvector multiplied by a
388: complex number of modulus one.
390: Input/Output Parameter:
391: . x - the computed vector
392: */
393: PetscErrorCode FixSign(Vec x)
394: {
395: PetscMPIInt rank;
396: PetscScalar sign=0.0;
397: const PetscScalar *xx;
400: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
401: if (!rank) {
402: VecGetArrayRead(x,&xx);
403: sign = *xx/PetscAbsScalar(*xx);
404: VecRestoreArrayRead(x,&xx);
405: }
406: MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
407: VecScale(x,1.0/sign);
408: PetscFunctionReturn(0);
409: }
411: /*TEST
413: build:
414: requires: !single
416: test:
417: suffix: 1
418: args: -nep_target 4
419: filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
420: requires: !single
422: testset:
423: args: -nep_type nleigs -nep_target 10 -nep_nev 4
424: test:
425: suffix: 2
426: filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
427: args: -rg_interval_endpoints 4,200
428: requires: !single !complex
429: test:
430: suffix: 2_complex
431: filter: sed -e "s/[+-]0.[0-9]*i//" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
432: output_file: output/ex20_2.out
433: args: -rg_interval_endpoints 4,200,-.1,.1 -nep_target_real
434: requires: !single complex
436: testset:
437: args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_ncv 18 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
438: test:
439: suffix: 3
440: filter: sed -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
441: output_file: output/ex20_2.out
442: args: -rg_interval_endpoints 4,200
443: requires: !single !complex
444: test:
445: suffix: 3_complex
446: filter: sed -e "s/[+-]0.[0-9]*i//" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/0\.\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
447: output_file: output/ex20_2.out
448: args: -rg_interval_endpoints 4,200,-.1,.1
449: requires: !single complex
451: test:
452: suffix: 4
453: args: -n 256 -nep_nev 2 -nep_target 10
454: filter: sed -e "s/[+-]0\.0*i//g" -e "s/[0-9]\.[0-9]*e-\([0-9]*\)/removed/g" -e "s/ Number of NEP iterations = \([0-9]*\)/ Number of NEP iterations = /"
455: requires: !single
457: TEST*/