Actual source code: ex20.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2013, 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*,MatStructure*,void*);
41: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat*,MatStructure*,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,maxf,nconv;
64: PetscReal re,im,abstol,rtol,stol,norm,error;
65: PetscBool draw_sol;
68: SlepcInitialize(&argc,&argv,(char*)0,help);
69: PetscOptionsGetInt(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,"-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,0,1e-9,0,0,0);
120: NEPSetDimensions(nep,1,0,0);
121: NEPSetLagPreconditioner(nep,0);
123: /*
124: Set solver parameters at runtime
125: */
126: NEPSetFromOptions(nep);
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Initialize application
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: /*
133: Evaluate initial guess
134: */
135: MatGetVecs(F,&x,NULL);
136: FormInitialGuess(x);
137: NEPSetInitialSpace(nep,1,&x);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Solve the eigensystem
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: NEPSolve(nep);
144: NEPGetIterationNumber(nep,&its);
145: PetscPrintf(PETSC_COMM_WORLD," Number of NEP iterations = %D\n\n",its);
147: /*
148: Optional: Get some information from the solver and display it
149: */
150: NEPGetType(nep,&type);
151: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type);
152: NEPGetDimensions(nep,&nev,NULL,NULL);
153: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
154: NEPGetTolerances(nep,&abstol,&rtol,&stol,&maxit,&maxf);
155: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: atol=%G, rtol=%G, stol=%G, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Display solution and clean up
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: /*
162: Get number of converged approximate eigenpairs
163: */
164: NEPGetConverged(nep,&nconv);
165: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);
167: if (nconv>0) {
168: /*
169: Display eigenvalues and relative errors
170: */
171: PetscPrintf(PETSC_COMM_WORLD,
172: " k ||T(k)x|| error\n"
173: " ----------------- ------------------ ------------------\n");
174: for (i=0;i<nconv;i++) {
175: /*
176: Get converged eigenpairs (in this example they are always real)
177: */
178: NEPGetEigenpair(nep,i,&lambda,x);
179: FixSign(x);
180: /*
181: Compute residual norm and error
182: */
183: NEPComputeRelativeError(nep,i,&norm);
184: CheckSolution(lambda,x,&error,&ctx);
186: #if defined(PETSC_USE_COMPLEX)
187: re = PetscRealPart(lambda);
188: im = PetscImaginaryPart(lambda);
189: #else
190: re = lambda;
191: im = 0.0;
192: #endif
193: if (im!=0.0) {
194: PetscPrintf(PETSC_COMM_WORLD," %9F%+9F j %12G %12G\n",re,im,norm,error);
195: } else {
196: PetscPrintf(PETSC_COMM_WORLD," %12F %12G %12G\n",re,norm,error);
197: }
198: if (draw_sol) {
199: PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
200: VecView(x,PETSC_VIEWER_DRAW_WORLD);
201: }
202: }
203: PetscPrintf(PETSC_COMM_WORLD,"\n");
204: }
206: NEPDestroy(&nep);
207: MatDestroy(&F);
208: MatDestroy(&J);
209: VecDestroy(&x);
210: SlepcFinalize();
211: return 0;
212: }
214: /* ------------------------------------------------------------------- */
217: /*
218: FormInitialGuess - Computes initial guess.
220: Input/Output Parameter:
221: . x - the solution vector
222: */
223: PetscErrorCode FormInitialGuess(Vec x)
224: {
228: VecSet(x,1.0);
229: return(0);
230: }
232: /* ------------------------------------------------------------------- */
235: /*
236: FormFunction - Computes Function matrix T(lambda)
238: Input Parameters:
239: . nep - the NEP context
240: . lambda - the scalar argument
241: . ctx - optional user-defined context, as set by NEPSetJacobian()
243: Output Parameters:
244: . fun - Function matrix
245: . B - optionally different preconditioning matrix
246: . flg - flag indicating matrix structure
247: */
248: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat *fun,Mat *B,MatStructure *flg,void *ctx)
249: {
251: ApplicationCtx *user = (ApplicationCtx*)ctx;
252: PetscScalar A[3],c,d;
253: PetscReal h;
254: PetscInt i,n,j[3],Istart,Iend;
255: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
258: /*
259: Compute Function entries and insert into matrix
260: */
261: MatGetSize(*fun,&n,NULL);
262: MatGetOwnershipRange(*fun,&Istart,&Iend);
263: if (Istart==0) FirstBlock=PETSC_TRUE;
264: if (Iend==n) LastBlock=PETSC_TRUE;
265: h = user->h;
266: c = user->kappa/(lambda-user->kappa);
267: d = n;
269: /*
270: Interior grid points
271: */
272: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
273: j[0] = i-1; j[1] = i; j[2] = i+1;
274: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
275: MatSetValues(*fun,1,&i,3,j,A,INSERT_VALUES);
276: }
278: /*
279: Boundary points
280: */
281: if (FirstBlock) {
282: i = 0;
283: j[0] = 0; j[1] = 1;
284: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
285: MatSetValues(*fun,1,&i,2,j,A,INSERT_VALUES);
286: }
288: if (LastBlock) {
289: i = n-1;
290: j[0] = n-2; j[1] = n-1;
291: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
292: MatSetValues(*fun,1,&i,2,j,A,INSERT_VALUES);
293: }
295: /*
296: Assemble matrix
297: */
298: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
299: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
300: if (*fun != *B) {
301: MatAssemblyBegin(*fun,MAT_FINAL_ASSEMBLY);
302: MatAssemblyEnd(*fun,MAT_FINAL_ASSEMBLY);
303: }
304: return(0);
305: }
307: /* ------------------------------------------------------------------- */
310: /*
311: FormJacobian - Computes Jacobian matrix T'(lambda)
313: Input Parameters:
314: . nep - the NEP context
315: . lambda - the scalar argument
316: . ctx - optional user-defined context, as set by NEPSetJacobian()
318: Output Parameters:
319: . jac - Jacobian matrix
320: . B - optionally different preconditioning matrix
321: . flg - flag indicating matrix structure
322: */
323: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat *jac,MatStructure *flg,void *ctx)
324: {
326: ApplicationCtx *user = (ApplicationCtx*)ctx;
327: PetscScalar A[3],c;
328: PetscReal h;
329: PetscInt i,n,j[3],Istart,Iend;
330: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
333: /*
334: Compute Jacobian entries and insert into matrix
335: */
336: MatGetSize(*jac,&n,NULL);
337: MatGetOwnershipRange(*jac,&Istart,&Iend);
338: if (Istart==0) FirstBlock=PETSC_TRUE;
339: if (Iend==n) LastBlock=PETSC_TRUE;
340: h = user->h;
341: c = user->kappa/(lambda-user->kappa);
343: /*
344: Interior grid points
345: */
346: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
347: j[0] = i-1; j[1] = i; j[2] = i+1;
348: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
349: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
350: }
352: /*
353: Boundary points
354: */
355: if (FirstBlock) {
356: i = 0;
357: j[0] = 0; j[1] = 1;
358: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
359: MatSetValues(*jac,1,&i,2,j,A,INSERT_VALUES);
360: }
362: if (LastBlock) {
363: i = n-1;
364: j[0] = n-2; j[1] = n-1;
365: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
366: MatSetValues(*jac,1,&i,2,j,A,INSERT_VALUES);
367: }
369: /*
370: Assemble matrix
371: */
372: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
373: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
374: return(0);
375: }
377: /* ------------------------------------------------------------------- */
380: /*
381: CheckSolution - Given a computed solution (lambda,x) check if it
382: satisfies the analytic solution.
384: Input Parameters:
385: + lambda - the computed eigenvalue
386: - y - the computed eigenvector
388: Output Parameter:
389: . error - norm of difference between the computed and exact eigenvector
390: */
391: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
392: {
394: PetscScalar nu,*uu;
395: PetscInt i,n,Istart,Iend;
396: PetscReal x;
397: Vec u;
398: ApplicationCtx *user = (ApplicationCtx*)ctx;
401: nu = PetscSqrtScalar(lambda);
402: VecDuplicate(y,&u);
403: VecGetSize(u,&n);
404: VecGetOwnershipRange(y,&Istart,&Iend);
405: VecGetArray(u,&uu);
406: for (i=Istart;i<Iend;i++) {
407: x = (i+1)*user->h;
408: uu[i-Istart] = sin(nu*x);
409: }
410: VecRestoreArray(u,&uu);
411: VecNormalize(u,NULL);
412: VecAXPY(u,-1.0,y);
413: VecNorm(u,NORM_2,error);
414: VecDestroy(&u);
415: return(0);
416: }
418: /* ------------------------------------------------------------------- */
421: /*
422: FixSign - Force the eigenfunction to be real and positive, since
423: some eigensolvers may return the eigenvector multiplied by a
424: complex number of modulus one.
426: Input/Output Parameter:
427: . x - the computed vector
428: */
429: PetscErrorCode FixSign(Vec x)
430: {
432: PetscMPIInt rank;
433: PetscScalar sign,*xx;
436: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
437: if (!rank) {
438: VecGetArray(x,&xx);
439: sign = *xx/PetscAbsScalar(*xx);
440: VecRestoreArray(x,&xx);
441: }
442: MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
443: VecScale(x,1.0/sign);
444: return(0);
445: }