Actual source code: ex10.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7:
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[] = "Illustrates the use of shell spectral transformations. "
23: "The problem to be solved is the same as ex1.c and"
24: "corresponds to the Laplacian operator in 1 dimension.\n\n"
25: "The command line options are:\n"
26: " -n <n>, where <n> = number of grid subdivisions = matrix dimension.\n\n";
28: #include slepceps.h
30: /* Define context for user-provided spectral transformation */
31: typedef struct {
32: KSP ksp;
33: } SampleShellST;
35: /* Declare routines for user-provided spectral transformation */
36: PetscErrorCode SampleShellSTCreate(SampleShellST**);
37: PetscErrorCode SampleShellSTSetUp(SampleShellST*,ST);
38: PetscErrorCode SampleShellSTApply(void*,Vec,Vec);
39: PetscErrorCode SampleShellSTBackTransform(void*,PetscScalar*,PetscScalar*);
40: PetscErrorCode SampleShellSTDestroy(SampleShellST*);
44: int main( int argc, char **argv )
45: {
46: Mat A; /* operator matrix */
47: EPS eps; /* eigenproblem solver context */
48: ST st; /* spectral transformation context */
49: SampleShellST *shell; /* user-defined spectral transform context */
50: const EPSType type;
51: PetscReal error, tol, re, im;
52: PetscScalar kr, ki;
54: PetscInt n=30, i, col[3], Istart, Iend, FirstBlock=0, LastBlock=0, nev, maxit, its, nconv;
55: PetscScalar value[3];
56: PetscTruth isShell;
58: SlepcInitialize(&argc,&argv,(char*)0,help);
60: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
61: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%d\n\n",n);
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Compute the operator matrix that defines the eigensystem, Ax=kx
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: MatCreate(PETSC_COMM_WORLD,&A);
68: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
69: MatSetFromOptions(A);
70:
71: MatGetOwnershipRange(A,&Istart,&Iend);
72: if (Istart==0) FirstBlock=PETSC_TRUE;
73: if (Iend==n) LastBlock=PETSC_TRUE;
74: value[0]=-1.0; value[1]=2.0; value[2]=-1.0;
75: for( i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++ ) {
76: col[0]=i-1; col[1]=i; col[2]=i+1;
77: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
78: }
79: if (LastBlock) {
80: i=n-1; col[0]=n-2; col[1]=n-1;
81: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
82: }
83: if (FirstBlock) {
84: i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
85: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
86: }
88: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Create the eigensolver and set various options
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: /*
96: Create eigensolver context
97: */
98: EPSCreate(PETSC_COMM_WORLD,&eps);
100: /*
101: Set operators. In this case, it is a standard eigenvalue problem
102: */
103: EPSSetOperators(eps,A,PETSC_NULL);
104: EPSSetProblemType(eps,EPS_HEP);
106: /*
107: Set solver parameters at runtime
108: */
109: EPSSetFromOptions(eps);
111: /*
112: Initialize shell spectral transformation if selected by user
113: */
114: EPSGetST(eps,&st);
115: PetscTypeCompare((PetscObject)st,STSHELL,&isShell);
116: if (isShell) {
117: /* (Optional) Create a context for the user-defined spectral tranform;
118: this context can be defined to contain any application-specific data. */
119: SampleShellSTCreate(&shell);
121: /* (Required) Set the user-defined routine for applying the operator */
122: STShellSetApply(st,SampleShellSTApply);
123: STShellSetContext(st,shell);
125: /* (Optional) Set the user-defined routine for back-transformation */
126: STShellSetBackTransform(st,SampleShellSTBackTransform);
128: /* (Optional) Set a name for the transformation, used for STView() */
129: STShellSetName(st,"MyTransformation");
131: /* (Optional) Do any setup required for the new transformation */
132: SampleShellSTSetUp(shell,st);
133: }
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Solve the eigensystem
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: EPSSolve(eps);
140: EPSGetIterationNumber(eps, &its);
141: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);
143: /*
144: Optional: Get some information from the solver and display it
145: */
146: EPSGetType(eps,&type);
147: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
148: EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
149: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
150: EPSGetTolerances(eps,&tol,&maxit);
151: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);
153: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: Display solution and clean up
155: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157: /*
158: Get number of converged approximate eigenpairs
159: */
160: EPSGetConverged(eps,&nconv);
161: PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);
163: if (nconv>0) {
164: /*
165: Display eigenvalues and relative errors
166: */
167: PetscPrintf(PETSC_COMM_WORLD,
168: " k ||Ax-kx||/||kx||\n"
169: " ----------------- ------------------\n" );
171: for( i=0; i<nconv; i++ ) {
172: /*
173: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
174: ki (imaginary part)
175: */
176: EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
177: /*
178: Compute the relative error associated to each eigenpair
179: */
180: EPSComputeRelativeError(eps,i,&error);
182: #ifdef PETSC_USE_COMPLEX
183: re = PetscRealPart(kr);
184: im = PetscImaginaryPart(kr);
185: #else
186: re = kr;
187: im = ki;
188: #endif
189: if (im!=0.0) {
190: PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
191: } else {
192: PetscPrintf(PETSC_COMM_WORLD," %12f %12g\n",re,error);
193: }
194: }
195: PetscPrintf(PETSC_COMM_WORLD,"\n" );
196: }
197:
198: /*
199: Free work space
200: */
201: if (isShell) {
202: SampleShellSTDestroy(shell);
203: }
204: EPSDestroy(eps);
205: MatDestroy(A);
206: SlepcFinalize();
207: return 0;
208: }
210: /***********************************************************************/
211: /* Routines for a user-defined shell spectral transformation */
212: /***********************************************************************/
216: /*
217: SampleShellSTCreate - This routine creates a user-defined
218: spectral transformation context.
220: Output Parameter:
221: . shell - user-defined spectral transformation context
222: */
223: PetscErrorCode SampleShellSTCreate(SampleShellST **shell)
224: {
225: SampleShellST *newctx;
229: PetscNew(SampleShellST,&newctx);
230: KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
231: KSPAppendOptionsPrefix(newctx->ksp,"st_");
232: *shell = newctx;
233: return(0);
234: }
235: /* ------------------------------------------------------------------- */
238: /*
239: SampleShellSTSetUp - This routine sets up a user-defined
240: spectral transformation context.
242: Input Parameters:
243: . shell - user-defined spectral transformation context
244: . st - spectral transformation context containing the operator matrices
246: Output Parameter:
247: . shell - fully set up user-defined transformation context
249: Notes:
250: In this example, the user-defined transformation is simply OP=A^-1.
251: Therefore, the eigenpairs converge in reversed order. The KSP object
252: used for the solution of linear systems with A is handled via the
253: user-defined context SampleShellST.
254: */
255: PetscErrorCode SampleShellSTSetUp(SampleShellST *shell,ST st)
256: {
257: Mat A,B;
261: STGetOperators(st,&A,&B);
262: if (B) { SETERRQ(0,"Warning: This transformation is not intended for generalized problems"); }
263: KSPSetOperators(shell->ksp,A,A,DIFFERENT_NONZERO_PATTERN);
264: KSPSetFromOptions(shell->ksp);
265: return(0);
266: }
267: /* ------------------------------------------------------------------- */
270: /*
271: SampleShellSTApply - This routine demonstrates the use of a
272: user-provided spectral transformation.
274: Input Parameters:
275: . ctx - optional user-defined context, as set by STShellSetContext()
276: . x - input vector
278: Output Parameter:
279: . y - output vector
281: Notes:
282: The transformation implemented in this code is just OP=A^-1 and
283: therefore it is of little use, merely as an example of working with
284: a STSHELL.
285: */
286: PetscErrorCode SampleShellSTApply(void *ctx,Vec x,Vec y)
287: {
288: SampleShellST *shell = (SampleShellST*)ctx;
292: KSPSolve(shell->ksp,x,y);
293: return(0);
294: }
295: /* ------------------------------------------------------------------- */
298: /*
299: SampleShellSTBackTransform - This routine demonstrates the use of a
300: user-provided spectral transformation.
302: Input Parameters:
303: . ctx - optional user-defined context, as set by STShellSetContext()
304: . eigr - pointer to real part of eigenvalues
305: . eigi - pointer to imaginary part of eigenvalues
307: Output Parameters:
308: . eigr - modified real part of eigenvalues
309: . eigi - modified imaginary part of eigenvalues
311: Notes:
312: This code implements the back transformation of eigenvalues in
313: order to retrieve the eigenvalues of the original problem. In this
314: example, simply set k_i = 1/k_i.
315: */
316: PetscErrorCode SampleShellSTBackTransform(void *ctx,PetscScalar *eigr,PetscScalar *eigi)
317: {
319: *eigr = 1.0 / *eigr;
320: return(0);
321: }
322: /* ------------------------------------------------------------------- */
325: /*
326: SampleShellSTDestroy - This routine destroys a user-defined
327: spectral transformation context.
329: Input Parameter:
330: . shell - user-defined spectral transformation context
331: */
332: PetscErrorCode SampleShellSTDestroy(SampleShellST *shell)
333: {
337: KSPDestroy(shell->ksp);
338: PetscFree(shell);
339: return(0);
340: }