Actual source code: ex20.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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);if (ierr) return ierr;
56: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
57: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%D\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 = %D\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: %D\n",nev);
140: NEPGetTolerances(nep,&tol,&maxit);
141: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%g, maxit=%D\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: %D\n\n",nconv);
153: if (nconv>0) {
154: /*
155: Display eigenvalues and relative errors
156: */
157: 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) {
180: PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g %12g\n",(double)re,(double)im,(double)norm,(double)error);
181: } else {
182: PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)error);
183: }
184: if (draw_sol) {
185: PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
186: VecView(x,PETSC_VIEWER_DRAW_WORLD);
187: }
188: }
189: PetscPrintf(PETSC_COMM_WORLD,"\n");
190: }
192: NEPDestroy(&nep);
193: MatDestroy(&F);
194: MatDestroy(&J);
195: VecDestroy(&x);
196: SlepcFinalize();
197: return ierr;
198: }
200: /* ------------------------------------------------------------------- */
201: /*
202: FormInitialGuess - Computes initial guess.
204: Input/Output Parameter:
205: . x - the solution vector
206: */
207: PetscErrorCode FormInitialGuess(Vec x)
208: {
212: VecSet(x,1.0);
213: return(0);
214: }
216: /* ------------------------------------------------------------------- */
217: /*
218: FormFunction - Computes Function matrix T(lambda)
220: Input Parameters:
221: . nep - the NEP context
222: . lambda - the scalar argument
223: . ctx - optional user-defined context, as set by NEPSetFunction()
225: Output Parameters:
226: . fun - Function matrix
227: . B - optionally different preconditioning matrix
228: */
229: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
230: {
232: ApplicationCtx *user = (ApplicationCtx*)ctx;
233: PetscScalar A[3],c,d;
234: PetscReal h;
235: PetscInt i,n,j[3],Istart,Iend;
236: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
239: /*
240: Compute Function entries and insert into matrix
241: */
242: MatGetSize(fun,&n,NULL);
243: MatGetOwnershipRange(fun,&Istart,&Iend);
244: if (Istart==0) FirstBlock=PETSC_TRUE;
245: if (Iend==n) LastBlock=PETSC_TRUE;
246: h = user->h;
247: c = user->kappa/(lambda-user->kappa);
248: d = n;
250: /*
251: Interior grid points
252: */
253: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
254: j[0] = i-1; j[1] = i; j[2] = i+1;
255: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
256: MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
257: }
259: /*
260: Boundary points
261: */
262: if (FirstBlock) {
263: i = 0;
264: j[0] = 0; j[1] = 1;
265: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
266: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
267: }
269: if (LastBlock) {
270: i = n-1;
271: j[0] = n-2; j[1] = n-1;
272: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
273: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
274: }
276: /*
277: Assemble matrix
278: */
279: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
280: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
281: if (fun != B) {
282: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
283: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
284: }
285: return(0);
286: }
288: /* ------------------------------------------------------------------- */
289: /*
290: FormJacobian - Computes Jacobian matrix T'(lambda)
292: Input Parameters:
293: . nep - the NEP context
294: . lambda - the scalar argument
295: . ctx - optional user-defined context, as set by NEPSetJacobian()
297: Output Parameters:
298: . jac - Jacobian matrix
299: */
300: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
301: {
303: ApplicationCtx *user = (ApplicationCtx*)ctx;
304: PetscScalar A[3],c;
305: PetscReal h;
306: PetscInt i,n,j[3],Istart,Iend;
307: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
310: /*
311: Compute Jacobian entries and insert into matrix
312: */
313: MatGetSize(jac,&n,NULL);
314: MatGetOwnershipRange(jac,&Istart,&Iend);
315: if (Istart==0) FirstBlock=PETSC_TRUE;
316: if (Iend==n) LastBlock=PETSC_TRUE;
317: h = user->h;
318: c = user->kappa/(lambda-user->kappa);
320: /*
321: Interior grid points
322: */
323: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
324: j[0] = i-1; j[1] = i; j[2] = i+1;
325: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
326: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
327: }
329: /*
330: Boundary points
331: */
332: if (FirstBlock) {
333: i = 0;
334: j[0] = 0; j[1] = 1;
335: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
336: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
337: }
339: if (LastBlock) {
340: i = n-1;
341: j[0] = n-2; j[1] = n-1;
342: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
343: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
344: }
346: /*
347: Assemble matrix
348: */
349: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
350: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
351: return(0);
352: }
354: /* ------------------------------------------------------------------- */
355: /*
356: CheckSolution - Given a computed solution (lambda,x) check if it
357: satisfies the analytic solution.
359: Input Parameters:
360: + lambda - the computed eigenvalue
361: - y - the computed eigenvector
363: Output Parameter:
364: . error - norm of difference between the computed and exact eigenvector
365: */
366: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
367: {
369: PetscScalar nu,*uu;
370: PetscInt i,n,Istart,Iend;
371: PetscReal x;
372: Vec u;
373: ApplicationCtx *user = (ApplicationCtx*)ctx;
376: nu = PetscSqrtScalar(lambda);
377: VecDuplicate(y,&u);
378: VecGetSize(u,&n);
379: VecGetOwnershipRange(y,&Istart,&Iend);
380: VecGetArray(u,&uu);
381: for (i=Istart;i<Iend;i++) {
382: x = (i+1)*user->h;
383: uu[i-Istart] = PetscSinReal(nu*x);
384: }
385: VecRestoreArray(u,&uu);
386: VecNormalize(u,NULL);
387: VecAXPY(u,-1.0,y);
388: VecNorm(u,NORM_2,error);
389: VecDestroy(&u);
390: return(0);
391: }
393: /* ------------------------------------------------------------------- */
394: /*
395: FixSign - Force the eigenfunction to be real and positive, since
396: some eigensolvers may return the eigenvector multiplied by a
397: complex number of modulus one.
399: Input/Output Parameter:
400: . x - the computed vector
401: */
402: PetscErrorCode FixSign(Vec x)
403: {
404: PetscErrorCode ierr;
405: PetscMPIInt rank;
406: PetscScalar sign=0.0;
407: const PetscScalar *xx;
410: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
411: if (!rank) {
412: VecGetArrayRead(x,&xx);
413: sign = *xx/PetscAbsScalar(*xx);
414: VecRestoreArrayRead(x,&xx);
415: }
416: MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
417: VecScale(x,1.0/sign);
418: return(0);
419: }
421: /*TEST
423: build:
424: requires: !single
426: test:
427: suffix: 1
428: args: -nep_target 4
429: 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 = /"
430: requires: !single
432: testset:
433: args: -nep_type nleigs -nep_target 10 -nep_nev 4
434: test:
435: suffix: 2
436: 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 = /"
437: args: -rg_interval_endpoints 4,900
438: requires: !single !complex
439: test:
440: suffix: 2_complex
441: 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 = /"
442: output_file: output/ex20_2.out
443: args: -rg_interval_endpoints 4,900,-.1,.1
444: requires: !single complex
446: testset:
447: args: -nep_type nleigs -nep_target 10 -nep_nev 4 -nep_two_sided {{0 1}} -nep_nleigs_full_basis
448: test:
449: suffix: 3
450: 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 = /"
451: output_file: output/ex20_2.out
452: args: -rg_interval_endpoints 4,900
453: requires: !single !complex
454: test:
455: suffix: 3_complex
456: 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 = /"
457: output_file: output/ex20_2.out
458: args: -rg_interval_endpoints 4,900,-.1,.1
459: requires: !single complex
461: TEST*/