Actual source code: ex20.c
slepc-3.6.3 2016-03-29
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, 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,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,PETSC_DEFAULT,1e-9,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
120: NEPSetDimensions(nep,1,PETSC_DEFAULT,PETSC_DEFAULT);
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: MatCreateVecs(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",(double)abstol,(double)rtol,(double)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,NULL,x,NULL);
179: FixSign(x);
180: /*
181: Compute residual norm and error
182: */
183: NEPComputeError(nep,i,NEP_ERROR_RELATIVE,&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",(double)re,(double)im,(double)norm,(double)error);
195: } else {
196: PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)norm,(double)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: */
247: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
248: {
250: ApplicationCtx *user = (ApplicationCtx*)ctx;
251: PetscScalar A[3],c,d;
252: PetscReal h;
253: PetscInt i,n,j[3],Istart,Iend;
254: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
257: /*
258: Compute Function entries and insert into matrix
259: */
260: MatGetSize(fun,&n,NULL);
261: MatGetOwnershipRange(fun,&Istart,&Iend);
262: if (Istart==0) FirstBlock=PETSC_TRUE;
263: if (Iend==n) LastBlock=PETSC_TRUE;
264: h = user->h;
265: c = user->kappa/(lambda-user->kappa);
266: d = n;
268: /*
269: Interior grid points
270: */
271: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
272: j[0] = i-1; j[1] = i; j[2] = i+1;
273: A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
274: MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES);
275: }
277: /*
278: Boundary points
279: */
280: if (FirstBlock) {
281: i = 0;
282: j[0] = 0; j[1] = 1;
283: A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
284: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
285: }
287: if (LastBlock) {
288: i = n-1;
289: j[0] = n-2; j[1] = n-1;
290: A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
291: MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES);
292: }
294: /*
295: Assemble matrix
296: */
297: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
298: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
299: if (fun != B) {
300: MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY);
301: MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY);
302: }
303: return(0);
304: }
306: /* ------------------------------------------------------------------- */
309: /*
310: FormJacobian - Computes Jacobian matrix T'(lambda)
312: Input Parameters:
313: . nep - the NEP context
314: . lambda - the scalar argument
315: . ctx - optional user-defined context, as set by NEPSetJacobian()
317: Output Parameters:
318: . jac - Jacobian matrix
319: . B - optionally different preconditioning matrix
320: */
321: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
322: {
324: ApplicationCtx *user = (ApplicationCtx*)ctx;
325: PetscScalar A[3],c;
326: PetscReal h;
327: PetscInt i,n,j[3],Istart,Iend;
328: PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
331: /*
332: Compute Jacobian entries and insert into matrix
333: */
334: MatGetSize(jac,&n,NULL);
335: MatGetOwnershipRange(jac,&Istart,&Iend);
336: if (Istart==0) FirstBlock=PETSC_TRUE;
337: if (Iend==n) LastBlock=PETSC_TRUE;
338: h = user->h;
339: c = user->kappa/(lambda-user->kappa);
341: /*
342: Interior grid points
343: */
344: for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
345: j[0] = i-1; j[1] = i; j[2] = i+1;
346: A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
347: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
348: }
350: /*
351: Boundary points
352: */
353: if (FirstBlock) {
354: i = 0;
355: j[0] = 0; j[1] = 1;
356: A[0] = -2.0*h/3.0; A[1] = -h/6.0;
357: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
358: }
360: if (LastBlock) {
361: i = n-1;
362: j[0] = n-2; j[1] = n-1;
363: A[0] = -h/6.0; A[1] = -h/3.0-c*c;
364: MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES);
365: }
367: /*
368: Assemble matrix
369: */
370: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
371: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
372: return(0);
373: }
375: /* ------------------------------------------------------------------- */
378: /*
379: CheckSolution - Given a computed solution (lambda,x) check if it
380: satisfies the analytic solution.
382: Input Parameters:
383: + lambda - the computed eigenvalue
384: - y - the computed eigenvector
386: Output Parameter:
387: . error - norm of difference between the computed and exact eigenvector
388: */
389: PetscErrorCode CheckSolution(PetscScalar lambda,Vec y,PetscReal *error,void *ctx)
390: {
392: PetscScalar nu,*uu;
393: PetscInt i,n,Istart,Iend;
394: PetscReal x;
395: Vec u;
396: ApplicationCtx *user = (ApplicationCtx*)ctx;
399: nu = PetscSqrtScalar(lambda);
400: VecDuplicate(y,&u);
401: VecGetSize(u,&n);
402: VecGetOwnershipRange(y,&Istart,&Iend);
403: VecGetArray(u,&uu);
404: for (i=Istart;i<Iend;i++) {
405: x = (i+1)*user->h;
406: uu[i-Istart] = PetscSinReal(nu*x);
407: }
408: VecRestoreArray(u,&uu);
409: VecNormalize(u,NULL);
410: VecAXPY(u,-1.0,y);
411: VecNorm(u,NORM_2,error);
412: VecDestroy(&u);
413: return(0);
414: }
416: /* ------------------------------------------------------------------- */
419: /*
420: FixSign - Force the eigenfunction to be real and positive, since
421: some eigensolvers may return the eigenvector multiplied by a
422: complex number of modulus one.
424: Input/Output Parameter:
425: . x - the computed vector
426: */
427: PetscErrorCode FixSign(Vec x)
428: {
429: PetscErrorCode ierr;
430: PetscMPIInt rank;
431: PetscScalar sign;
432: const PetscScalar *xx;
435: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
436: if (!rank) {
437: VecGetArrayRead(x,&xx);
438: sign = *xx/PetscAbsScalar(*xx);
439: VecRestoreArrayRead(x,&xx);
440: }
441: MPI_Bcast(&sign,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
442: VecScale(x,1.0/sign);
443: return(0);
444: }