Actual source code: ex20.c
slepc-3.19.1 2023-05-15
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: PetscFunctionBeginUser;
55: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
56: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
57: PetscCall(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: PetscCall(PetscOptionsGetBool(NULL,NULL,"-draw_sol",&draw_sol,NULL));
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create nonlinear eigensolver context
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Create matrix data structure; set Function evaluation routine
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
72: PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
73: PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
74: PetscCall(MatSetFromOptions(F));
75: PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
76: PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
77: PetscCall(MatSetUp(F));
79: /*
80: Set Function matrix data structure and default Function evaluation
81: routine
82: */
83: PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Create matrix data structure; set Jacobian evaluation routine
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
90: PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
91: PetscCall(MatSetFromOptions(J));
92: PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
93: PetscCall(MatMPIAIJSetPreallocation(J,3,NULL,1,NULL));
94: PetscCall(MatSetUp(J));
96: /*
97: Set Jacobian matrix data structure and default Jacobian evaluation
98: routine
99: */
100: PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Customize nonlinear solver; set runtime options
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: PetscCall(NEPSetTolerances(nep,1e-9,PETSC_DEFAULT));
107: PetscCall(NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT));
109: /*
110: Set solver parameters at runtime
111: */
112: PetscCall(NEPSetFromOptions(nep));
114: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115: Initialize application
116: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: /*
119: Evaluate initial guess
120: */
121: PetscCall(MatCreateVecs(F,&x,NULL));
122: PetscCall(FormInitialGuess(x));
123: PetscCall(NEPSetInitialSpace(nep,1,&x));
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Solve the eigensystem
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129: PetscCall(NEPSolve(nep));
130: PetscCall(NEPGetIterationNumber(nep,&its));
131: PetscCall(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: PetscCall(NEPGetType(nep,&type));
137: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
138: PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
139: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
140: PetscCall(NEPGetTolerances(nep,&tol,&maxit));
141: PetscCall(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: PetscCall(NEPGetConverged(nep,&nconv));
151: PetscCall(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: PetscCall(NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL));
165: PetscCall(FixSign(x));
166: /*
167: Compute residual norm and error
168: */
169: PetscCall(NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm));
170: PetscCall(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) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g %12g\n",(double)re,(double)im,(double)norm,(double)error));
180: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)error));
181: if (draw_sol) {
182: PetscCall(PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1));
183: PetscCall(VecView(x,PETSC_VIEWER_DRAW_WORLD));
184: }
185: }
186: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
187: }
189: PetscCall(NEPDestroy(&nep));
190: PetscCall(MatDestroy(&F));
191: PetscCall(MatDestroy(&J));
192: PetscCall(VecDestroy(&x));
193: PetscCall(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: {
206: PetscFunctionBeginUser;
207: PetscCall(VecSet(x,1.0));
208: PetscFunctionReturn(PETSC_SUCCESS);
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;
232: PetscFunctionBeginUser;
233: /*
234: Compute Function entries and insert into matrix
235: */
236: PetscCall(MatGetSize(fun,&n,NULL));
237: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
268: }
270: /*
271: Assemble matrix
272: */
273: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
274: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
275: if (fun != B) {
276: PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
277: PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
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;
302: PetscFunctionBeginUser;
303: /*
304: Compute Jacobian entries and insert into matrix
305: */
306: PetscCall(MatGetSize(jac,&n,NULL));
307: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
337: }
339: /*
340: Assemble matrix
341: */
342: PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
343: PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
344: PetscFunctionReturn(PETSC_SUCCESS);
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;
367: PetscFunctionBeginUser;
368: nu = PetscSqrtReal(PetscRealPart(lambda));
369: PetscCall(VecDuplicate(y,&u));
370: PetscCall(VecGetSize(u,&n));
371: PetscCall(VecGetOwnershipRange(y,&Istart,&Iend));
372: PetscCall(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: PetscCall(VecRestoreArray(u,&uu));
378: PetscCall(VecNormalize(u,NULL));
379: PetscCall(VecAXPY(u,-1.0,y));
380: PetscCall(VecNorm(u,NORM_2,error));
381: PetscCall(VecDestroy(&u));
382: PetscFunctionReturn(PETSC_SUCCESS);
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;
400: PetscFunctionBeginUser;
401: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
402: if (!rank) {
403: PetscCall(VecGetArrayRead(x,&xx));
404: sign = *xx/PetscAbsScalar(*xx);
405: PetscCall(VecRestoreArrayRead(x,&xx));
406: }
407: PetscCallMPI(MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD));
408: PetscCall(VecScale(x,1.0/sign));
409: PetscFunctionReturn(PETSC_SUCCESS);
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*/