Actual source code: ex20.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.\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;
55: SlepcInitialize(&argc,&argv,(char*)0,help);
56: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
57: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n);
58: ctx.h = 1.0/(PetscReal)n;
59: ctx.kappa = 1.0;
60: PetscOptionsGetBool(NULL,NULL,"-draw_sol",&draw_sol,NULL);
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create nonlinear eigensolver context
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: NEPCreate(PETSC_COMM_WORLD,&nep);
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create matrix data structure; set Function evaluation routine
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: MatCreate(PETSC_COMM_WORLD,&F);
73: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
74: MatSetFromOptions(F);
75: MatSeqAIJSetPreallocation(F,3,NULL);
76: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
77: MatSetUp(F);
79: /*
80: Set Function matrix data structure and default Function evaluation
81: routine
82: */
83: NEPSetFunction(nep,F,F,FormFunction,&ctx);
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Create matrix data structure; set Jacobian evaluation routine
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: MatCreate(PETSC_COMM_WORLD,&J);
90: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
91: MatSetFromOptions(J);
92: MatSeqAIJSetPreallocation(J,3,NULL);
93: MatMPIAIJSetPreallocation(J,3,NULL,1,NULL);
94: MatSetUp(J);
96: /*
97: Set Jacobian matrix data structure and default Jacobian evaluation
98: routine
99: */
100: NEPSetJacobian(nep,J,FormJacobian,&ctx);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Customize nonlinear solver; set runtime options
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: NEPSetTolerances(nep,1e-9,PETSC_DEFAULT);
107: NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT);
109: /*
110: Set solver parameters at runtime
111: */
112: NEPSetFromOptions(nep);
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Initialize application
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: /*
119: Evaluate initial guess
120: */
121: MatCreateVecs(F,&x,NULL);
122: FormInitialGuess(x);
123: NEPSetInitialSpace(nep,1,&x);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Solve the eigensystem
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: NEPSolve(nep);
130: NEPGetIterationNumber(nep,&its);
131: PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %" PetscInt_FMT "\n\n",its);
133: /*
134: Optional: Get some information from the solver and display it
135: */
136: NEPGetType(nep,&type);
137: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
138: NEPGetDimensions(nep,&nev,NULL,NULL);
139: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
140: NEPGetTolerances(nep,&tol,&maxit);
141: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit);
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Display solution and clean up
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: /*
148: Get number of converged approximate eigenpairs
149: */
150: NEPGetConverged(nep,&nconv);
151: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %" PetscInt_FMT "\n\n",nconv);
153: if (nconv>0) {
154: /*
155: Display eigenvalues and relative errors
156: */
157: PetscCall(PetscPrintf(PETSC_COMM_WORLD,
158: " k ||T(k)x|| error\n"
159: " ----------------- ------------------ ------------------\n"));
160: for (i=0;i<nconv;i++) {
161: /*
162: Get converged eigenpairs (in this example they are always real)
163: */
164: NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL);
165: FixSign(x);
166: /*
167: Compute residual norm and error
168: */
169: NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm);
170: CheckSolution(lambda,x,&error,&ctx);
172: #if defined(PETSC_USE_COMPLEX)
173: re = PetscRealPart(lambda);
174: im = PetscImaginaryPart(lambda);
175: #else
176: re = lambda;
177: im = 0.0;
178: #endif
179: if (im!=0.0) PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g %12g\n",(double)re,(double)im,(double)norm,(double)error);
180: else PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)error);
181: if (draw_sol) {
182: PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
183: VecView(x,PETSC_VIEWER_DRAW_WORLD);
184: }
185: }
186: PetscPrintf(PETSC_COMM_WORLD,"\n");
187: }
189: NEPDestroy(&nep);
190: MatDestroy(&F);
191: MatDestroy(&J);
192: VecDestroy(&x);
193: SlepcFinalize();
194: return 0;
195: }
197: /* ------------------------------------------------------------------- */
198: /*
199: FormInitialGuess - Computes initial guess.
201: Input/Output Parameter:
202: . x - the solution vector
203: */
204: PetscErrorCode FormInitialGuess(Vec x)
205: {
207: VecSet(x,1.0);
208: return 0;
209: }
211: /* ------------------------------------------------------------------- */
212: /*
213: FormFunction - Computes Function matrix T(lambda)
215: Input Parameters:
216: . nep - the NEP context
217: . lambda - the scalar argument
218: . ctx - optional user-defined context, as set by NEPSetFunction()
220: Output Parameters:
221: . fun - Function matrix
222: . B - optionally different preconditioning matrix
223: */
224: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
225: {
226: ApplicationCtx *user = (ApplicationCtx*)ctx;
227: PetscScalar A[3],c,d;
228: PetscReal h;
229: PetscInt i,n,j[3],Istart,Iend;
230: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
233: /*
234: Compute Function entries and insert into matrix
235: */
236: MatGetSize(fun,&n,NULL);
237: MatGetOwnershipRange(fun,&Istart,&Iend);
238: if (Istart==0) FirstBlock=PETSC_TRUE;
239: if (Iend==n) LastBlock=PETSC_TRUE;
240: h = user->h;
241: c = user->kappa/(lambda-user->kappa);
242: d = n;
244: /*
245: Interior grid points
246: */
247: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
248: j[0] = i-1; j[1] = i; j[2] = i+1;
249: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
250: MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
251: }
253: /*
254: Boundary points
255: */
256: if (FirstBlock) {
257: i = 0;
258: j[0] = 0; j[1] = 1;
259: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
260: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
261: }
263: if (LastBlock) {
264: i = n-1;
265: j[0] = n-2; j[1] = n-1;
266: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
267: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
268: }
270: /*
271: Assemble matrix
272: */
273: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
274: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
275: if (fun != B) {
276: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
277: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
278: }
279: return 0;
280: }
282: /* ------------------------------------------------------------------- */
283: /*
284: FormJacobian - Computes Jacobian matrix T'(lambda)
286: Input Parameters:
287: . nep - the NEP context
288: . lambda - the scalar argument
289: . ctx - optional user-defined context, as set by NEPSetJacobian()
291: Output Parameters:
292: . jac - Jacobian matrix
293: */
294: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
295: {
296: ApplicationCtx *user = (ApplicationCtx*)ctx;
297: PetscScalar A[3],c;
298: PetscReal h;
299: PetscInt i,n,j[3],Istart,Iend;
300: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
303: /*
304: Compute Jacobian entries and insert into matrix
305: */
306: MatGetSize(jac,&n,NULL);
307: MatGetOwnershipRange(jac,&Istart,&Iend);
308: if (Istart==0) FirstBlock=PETSC_TRUE;
309: if (Iend==n) LastBlock=PETSC_TRUE;
310: h = user->h;
311: c = user->kappa/(lambda-user->kappa);
313: /*
314: Interior grid points
315: */
316: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
317: j[0] = i-1; j[1] = i; j[2] = i+1;
318: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
319: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
320: }
322: /*
323: Boundary points
324: */
325: if (FirstBlock) {
326: i = 0;
327: j[0] = 0; j[1] = 1;
328: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
329: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
330: }
332: if (LastBlock) {
333: i = n-1;
334: j[0] = n-2; j[1] = n-1;
335: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
336: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
337: }
339: /*
340: Assemble matrix
341: */
342: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
343: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
344: return 0;
345: }
347: /* ------------------------------------------------------------------- */
348: /*
349: CheckSolution - Given a computed solution (lambda,x) check if it
350: satisfies the analytic solution.
352: Input Parameters:
353: + lambda - the computed eigenvalue
354: - y - the computed eigenvector
356: Output Parameter:
357: . error - norm of difference between the computed and exact eigenvector
358: */
359: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
360: {
361: PetscScalar *uu;
362: PetscInt i,n,Istart,Iend;
363: PetscReal nu,x;
364: Vec u;
365: ApplicationCtx *user = (ApplicationCtx*)ctx;
368: nu = PetscSqrtReal(PetscRealPart(lambda));
369: VecDuplicate(y,&u);
370: VecGetSize(u,&n);
371: VecGetOwnershipRange(y,&Istart,&Iend);
372: VecGetArray(u,&uu);
373: for (i=Istart;i<Iend;i++) {
374: x = (i+1)*user->h;
375: uu[i-Istart] = PetscSinReal(nu*x);
376: }
377: VecRestoreArray(u,&uu);
378: VecNormalize(u,NULL);
379: VecAXPY(u,-1.0,y);
380: VecNorm(u,NORM_2,error);
381: VecDestroy(&u);
382: return 0;
383: }
385: /* ------------------------------------------------------------------- */
386: /*
387: FixSign - Force the eigenfunction to be real and positive, since
388: some eigensolvers may return the eigenvector multiplied by a
389: complex number of modulus one.
391: Input/Output Parameter:
392: . x - the computed vector
393: */
394: PetscErrorCode FixSign(Vec x)
395: {
396: PetscMPIInt rank;
397: PetscScalar sign=0.0;
398: const PetscScalar *xx;
401: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
402: if (!rank) {
403: VecGetArrayRead(x,&xx);
404: sign = *xx/PetscAbsScalar(*xx);
405: VecRestoreArrayRead(x,&xx);
406: }
407: MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
408: VecScale(x,1.0/sign);
409: return 0;
410: }
412: /*TEST
414: build:
415: requires: !single
417: test:
418: suffix: 1
419: args: -nep_target 4
420: 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 = /"
421: requires: !single
423: testset:
424: args: -nep_type nleigs -nep_target 10 -nep_nev 4
425: test:
426: suffix: 2
427: 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 = /"
428: args: -rg_interval_endpoints 4,200
429: requires: !single !complex
430: test:
431: suffix: 2_complex
432: 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 = /"
433: output_file: output/ex20_2.out
434: args: -rg_interval_endpoints 4,200,-.1,.1 -nep_target_real
435: requires: !single complex
437: testset:
438: args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_ncv 18 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
439: test:
440: suffix: 3
441: 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 = /"
442: output_file: output/ex20_2.out
443: args: -rg_interval_endpoints 4,200
444: requires: !single !complex
445: test:
446: suffix: 3_complex
447: 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 = /"
448: output_file: output/ex20_2.out
449: args: -rg_interval_endpoints 4,200,-.1,.1
450: requires: !single complex
452: test:
453: suffix: 4
454: args: -n 256 -nep_nev 2 -nep_target 10
455: 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 = /"
456: requires: !single
458: TEST*/