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