Actual source code: ex20.c
slepc-3.7.4 2017-05-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2016, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Simple 1-D nonlinear eigenproblem.\n\n"
23: "The command line options are:\n"
24: " -n <n>, where <n> = number of grid subdivisions.\n"
25: " -draw_sol, to draw the computed solution.\n\n";
27: /*
28: Solve 1-D PDE
29: -u'' = lambda*u
30: on [0,1] subject to
31: u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
32: */
34: #include <slepcnep.h>
36: /*
37: User-defined routines
38: */
39: PetscErrorCode FormInitialGuess(Vec);
40: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
41: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
42: PetscErrorCode CheckSolution(PetscScalar,Vec,PetscReal*,void*);
43: PetscErrorCode FixSign(Vec);
45: /*
46: User-defined application context
47: */
48: typedef struct {
49: PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
50: PetscReal h; /* mesh spacing */
51: } ApplicationCtx;
55: int main(int argc,char **argv)
56: {
57: NEP nep; /* nonlinear eigensolver context */
58: Vec x; /* eigenvector */
59: PetscScalar lambda; /* eigenvalue */
60: Mat F,J; /* Function and Jacobian matrices */
61: ApplicationCtx ctx; /* user-defined context */
62: NEPType type;
63: PetscInt n=128,nev,i,its,maxit,nconv;
64: PetscReal re,im,tol,norm,error;
65: PetscBool draw_sol;
68: SlepcInitialize(&argc,&argv,(char*)0,help);
69: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
70: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%D\n\n",n);
71: ctx.h = 1.0/(PetscReal)n;
72: ctx.kappa = 1.0;
73: PetscOptionsHasName(NULL,NULL,"-draw_sol",&draw_sol);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create nonlinear eigensolver context
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: NEPCreate(PETSC_COMM_WORLD,&nep);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create matrix data structure; set Function evaluation routine
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: MatCreate(PETSC_COMM_WORLD,&F);
86: MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n);
87: MatSetFromOptions(F);
88: MatSeqAIJSetPreallocation(F,3,NULL);
89: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
90: MatSetUp(F);
92: /*
93: Set Function matrix data structure and default Function evaluation
94: routine
95: */
96: NEPSetFunction(nep,F,F,FormFunction,&ctx);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Create matrix data structure; set Jacobian evaluation routine
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: MatCreate(PETSC_COMM_WORLD,&J);
103: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
104: MatSetFromOptions(J);
105: MatSeqAIJSetPreallocation(J,3,NULL);
106: MatMPIAIJSetPreallocation(F,3,NULL,1,NULL);
107: MatSetUp(J);
109: /*
110: Set Jacobian matrix data structure and default Jacobian evaluation
111: routine
112: */
113: NEPSetJacobian(nep,J,FormJacobian,&ctx);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Customize nonlinear solver; set runtime options
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: NEPSetTolerances(nep,1e-9,PETSC_DEFAULT);
120: NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT);
122: /*
123: Set solver parameters at runtime
124: */
125: NEPSetFromOptions(nep);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Initialize application
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: /*
132: Evaluate initial guess
133: */
134: MatCreateVecs(F,&x,NULL);
135: FormInitialGuess(x);
136: NEPSetInitialSpace(nep,1,&x);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Solve the eigensystem
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: NEPSolve(nep);
143: NEPGetIterationNumber(nep,&its);
144: PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %D\n\n",its);
146: /*
147: Optional: Get some information from the solver and display it
148: */
149: NEPGetType(nep,&type);
150: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
151: NEPGetDimensions(nep,&nev,NULL,NULL);
152: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
153: NEPGetTolerances(nep,&tol,&maxit);
154: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%g, maxit=%D\n",(double)tol,maxit);
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Display solution and clean up
158: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: /*
161: Get number of converged approximate eigenpairs
162: */
163: NEPGetConverged(nep,&nconv);
164: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);
166: if (nconv>0) {
167: /*
168: Display eigenvalues and relative errors
169: */
170: PetscPrintf(PETSC_COMM_WORLD,
171: " k ||T(k)x|| error\n"
172: " ----------------- ------------------ ------------------\n");
173: for (i=0;i<nconv;i++) {
174: /*
175: Get converged eigenpairs (in this example they are always real)
176: */
177: NEPGetEigenpair(nep,i,&lambda,NULL,x,NULL);
178: FixSign(x);
179: /*
180: Compute residual norm and error
181: */
182: NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&norm);
183: CheckSolution(lambda,x,&error,&ctx);
185: #if defined(PETSC_USE_COMPLEX)
186: re = PetscRealPart(lambda);
187: im = PetscImaginaryPart(lambda);
188: #else
189: re = lambda;
190: im = 0.0;
191: #endif
192: if (im!=0.0) {
193: PetscPrintf(PETSC_COMM_WORLD," %9f%+9fi %12g %12g\n",(double)re,(double)im,(double)norm,(double)error);
194: } else {
195: PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)error);
196: }
197: if (draw_sol) {
198: PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
199: VecView(x,PETSC_VIEWER_DRAW_WORLD);
200: }
201: }
202: PetscPrintf(PETSC_COMM_WORLD,"\n");
203: }
205: NEPDestroy(&nep);
206: MatDestroy(&F);
207: MatDestroy(&J);
208: VecDestroy(&x);
209: SlepcFinalize();
210: return ierr;
211: }
213: /* ------------------------------------------------------------------- */
216: /*
217: FormInitialGuess - Computes initial guess.
219: Input/Output Parameter:
220: . x - the solution vector
221: */
222: PetscErrorCode FormInitialGuess(Vec x)
223: {
227: VecSet(x,1.0);
228: return(0);
229: }
231: /* ------------------------------------------------------------------- */
234: /*
235: FormFunction - Computes Function matrix T(lambda)
237: Input Parameters:
238: . nep - the NEP context
239: . lambda - the scalar argument
240: . ctx - optional user-defined context, as set by NEPSetFunction()
242: Output Parameters:
243: . fun - Function matrix
244: . B - optionally different preconditioning matrix
245: */
246: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
247: {
249: ApplicationCtx *user = (ApplicationCtx*)ctx;
250: PetscScalar A[3],c,d;
251: PetscReal h;
252: PetscInt i,n,j[3],Istart,Iend;
253: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
256: /*
257: Compute Function entries and insert into matrix
258: */
259: MatGetSize(fun,&n,NULL);
260: MatGetOwnershipRange(fun,&Istart,&Iend);
261: if (Istart==0) FirstBlock=PETSC_TRUE;
262: if (Iend==n) LastBlock=PETSC_TRUE;
263: h = user->h;
264: c = user->kappa/(lambda-user->kappa);
265: d = n;
267: /*
268: Interior grid points
269: */
270: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
271: j[0] = i-1; j[1] = i; j[2] = i+1;
272: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
273: MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
274: }
276: /*
277: Boundary points
278: */
279: if (FirstBlock) {
280: i = 0;
281: j[0] = 0; j[1] = 1;
282: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
283: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
284: }
286: if (LastBlock) {
287: i = n-1;
288: j[0] = n-2; j[1] = n-1;
289: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
290: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
291: }
293: /*
294: Assemble matrix
295: */
296: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
297: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
298: if (fun != B) {
299: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
300: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
301: }
302: return(0);
303: }
305: /* ------------------------------------------------------------------- */
308: /*
309: FormJacobian - Computes Jacobian matrix T'(lambda)
311: Input Parameters:
312: . nep - the NEP context
313: . lambda - the scalar argument
314: . ctx - optional user-defined context, as set by NEPSetJacobian()
316: Output Parameters:
317: . jac - Jacobian matrix
318: . B - optionally different preconditioning matrix
319: */
320: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
321: {
323: ApplicationCtx *user = (ApplicationCtx*)ctx;
324: PetscScalar A[3],c;
325: PetscReal h;
326: PetscInt i,n,j[3],Istart,Iend;
327: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
330: /*
331: Compute Jacobian entries and insert into matrix
332: */
333: MatGetSize(jac,&n,NULL);
334: MatGetOwnershipRange(jac,&Istart,&Iend);
335: if (Istart==0) FirstBlock=PETSC_TRUE;
336: if (Iend==n) LastBlock=PETSC_TRUE;
337: h = user->h;
338: c = user->kappa/(lambda-user->kappa);
340: /*
341: Interior grid points
342: */
343: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
344: j[0] = i-1; j[1] = i; j[2] = i+1;
345: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
346: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
347: }
349: /*
350: Boundary points
351: */
352: if (FirstBlock) {
353: i = 0;
354: j[0] = 0; j[1] = 1;
355: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
356: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
357: }
359: if (LastBlock) {
360: i = n-1;
361: j[0] = n-2; j[1] = n-1;
362: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
363: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
364: }
366: /*
367: Assemble matrix
368: */
369: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
370: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
371: return(0);
372: }
374: /* ------------------------------------------------------------------- */
377: /*
378: CheckSolution - Given a computed solution (lambda,x) check if it
379: satisfies the analytic solution.
381: Input Parameters:
382: + lambda - the computed eigenvalue
383: - y - the computed eigenvector
385: Output Parameter:
386: . error - norm of difference between the computed and exact eigenvector
387: */
388: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
389: {
391: PetscScalar nu,*uu;
392: PetscInt i,n,Istart,Iend;
393: PetscReal x;
394: Vec u;
395: ApplicationCtx *user = (ApplicationCtx*)ctx;
398: nu = PetscSqrtScalar(lambda);
399: VecDuplicate(y,&u);
400: VecGetSize(u,&n);
401: VecGetOwnershipRange(y,&Istart,&Iend);
402: VecGetArray(u,&uu);
403: for (i=Istart;i<Iend;i++) {
404: x = (i+1)*user->h;
405: uu[i-Istart] = PetscSinReal(nu*x);
406: }
407: VecRestoreArray(u,&uu);
408: VecNormalize(u,NULL);
409: VecAXPY(u,-1.0,y);
410: VecNorm(u,NORM_2,error);
411: VecDestroy(&u);
412: return(0);
413: }
415: /* ------------------------------------------------------------------- */
418: /*
419: FixSign - Force the eigenfunction to be real and positive, since
420: some eigensolvers may return the eigenvector multiplied by a
421: complex number of modulus one.
423: Input/Output Parameter:
424: . x - the computed vector
425: */
426: PetscErrorCode FixSign(Vec x)
427: {
428: PetscErrorCode ierr;
429: PetscMPIInt rank;
430: PetscScalar sign;
431: const PetscScalar *xx;
434: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
435: if (!rank) {
436: VecGetArrayRead(x,&xx);
437: sign = *xx/PetscAbsScalar(*xx);
438: VecRestoreArrayRead(x,&xx);
439: }
440: MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
441: VecScale(x,1.0/sign);
442: return(0);
443: }