Actual source code: ex10.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2011, Universitat 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(ST,Vec,Vec);
39: PetscErrorCode SampleShellSTBackTransform(ST,PetscInt,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 tol;
52: PetscScalar value[3];
53: PetscInt n=30,i,col[3],Istart,Iend,FirstBlock=0,LastBlock=0,nev,maxit;
54: PetscBool isShell;
57: SlepcInitialize(&argc,&argv,(char*)0,help);
59: PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
60: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%D\n\n",n);
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Compute the operator matrix that defines the eigensystem, Ax=kx
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: MatCreate(PETSC_COMM_WORLD,&A);
67: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
68: MatSetFromOptions(A);
69:
70: MatGetOwnershipRange(A,&Istart,&Iend);
71: if (Istart==0) FirstBlock=PETSC_TRUE;
72: if (Iend==n) LastBlock=PETSC_TRUE;
73: value[0]=-1.0; value[1]=2.0; value[2]=-1.0;
74: for (i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++) {
75: col[0]=i-1; col[1]=i; col[2]=i+1;
76: MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
77: }
78: if (LastBlock) {
79: i=n-1; col[0]=n-2; col[1]=n-1;
80: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
81: }
82: if (FirstBlock) {
83: i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
84: MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
85: }
87: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Create the eigensolver and set various options
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94: /*
95: Create eigensolver context
96: */
97: EPSCreate(PETSC_COMM_WORLD,&eps);
99: /*
100: Set operators. In this case, it is a standard eigenvalue problem
101: */
102: EPSSetOperators(eps,A,PETSC_NULL);
103: EPSSetProblemType(eps,EPS_HEP);
105: /*
106: Set solver parameters at runtime
107: */
108: EPSSetFromOptions(eps);
110: /*
111: Initialize shell spectral transformation if selected by user
112: */
113: EPSGetST(eps,&st);
114: PetscTypeCompare((PetscObject)st,STSHELL,&isShell);
115: if (isShell) {
116: /* (Optional) Create a context for the user-defined spectral tranform;
117: this context can be defined to contain any application-specific data. */
118: SampleShellSTCreate(&shell);
120: /* (Required) Set the user-defined routine for applying the operator */
121: STShellSetApply(st,SampleShellSTApply);
122: STShellSetContext(st,shell);
124: /* (Optional) Set the user-defined routine for back-transformation */
125: STShellSetBackTransform(st,SampleShellSTBackTransform);
127: /* (Optional) Set a name for the transformation, used for STView() */
128: PetscObjectSetName((PetscObject)st,"MyTransformation");
130: /* (Optional) Do any setup required for the new transformation */
131: SampleShellSTSetUp(shell,st);
132: }
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Solve the eigensystem
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: EPSSolve(eps);
140: /*
141: Optional: Get some information from the solver and display it
142: */
143: EPSGetType(eps,&type);
144: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
145: EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
146: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
147: EPSGetTolerances(eps,&tol,&maxit);
148: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4G, maxit=%D\n",tol,maxit);
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Display solution and clean up
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: EPSPrintSolution(eps,PETSC_NULL);
155: if (isShell) {
156: SampleShellSTDestroy(shell);
157: }
158: EPSDestroy(&eps);
159: MatDestroy(&A);
160: SlepcFinalize();
161: return 0;
162: }
164: /***********************************************************************/
165: /* Routines for a user-defined shell spectral transformation */
166: /***********************************************************************/
170: /*
171: SampleShellSTCreate - This routine creates a user-defined
172: spectral transformation context.
174: Output Parameter:
175: . shell - user-defined spectral transformation context
176: */
177: PetscErrorCode SampleShellSTCreate(SampleShellST **shell)
178: {
179: SampleShellST *newctx;
183: PetscNew(SampleShellST,&newctx);
184: KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
185: KSPAppendOptionsPrefix(newctx->ksp,"st_");
186: *shell = newctx;
187: return(0);
188: }
189: /* ------------------------------------------------------------------- */
192: /*
193: SampleShellSTSetUp - This routine sets up a user-defined
194: spectral transformation context.
196: Input Parameters:
197: . shell - user-defined spectral transformation context
198: . st - spectral transformation context containing the operator matrices
200: Output Parameter:
201: . shell - fully set up user-defined transformation context
203: Notes:
204: In this example, the user-defined transformation is simply OP=A^-1.
205: Therefore, the eigenpairs converge in reversed order. The KSP object
206: used for the solution of linear systems with A is handled via the
207: user-defined context SampleShellST.
208: */
209: PetscErrorCode SampleShellSTSetUp(SampleShellST *shell,ST st)
210: {
211: Mat A,B;
215: STGetOperators(st,&A,&B);
216: if (B) { PetscInfo(B,"This transformation is not intended for generalized problems, ignoring matrix B"); }
217: KSPSetOperators(shell->ksp,A,A,DIFFERENT_NONZERO_PATTERN);
218: KSPSetFromOptions(shell->ksp);
219: return(0);
220: }
221: /* ------------------------------------------------------------------- */
224: /*
225: SampleShellSTApply - This routine demonstrates the use of a
226: user-provided spectral transformation.
228: Input Parameters:
229: . ctx - optional user-defined context, as set by STShellSetContext()
230: . x - input vector
232: Output Parameter:
233: . y - output vector
235: Notes:
236: The transformation implemented in this code is just OP=A^-1 and
237: therefore it is of little use, merely as an example of working with
238: a STSHELL.
239: */
240: PetscErrorCode SampleShellSTApply(ST st,Vec x,Vec y)
241: {
242: SampleShellST *shell;
246: STShellGetContext(st,(void**)&shell);
247: KSPSolve(shell->ksp,x,y);
248: return(0);
249: }
250: /* ------------------------------------------------------------------- */
253: /*
254: SampleShellSTBackTransform - This routine demonstrates the use of a
255: user-provided spectral transformation.
257: Input Parameters:
258: . ctx - optional user-defined context, as set by STShellSetContext()
259: . eigr - pointer to real part of eigenvalues
260: . eigi - pointer to imaginary part of eigenvalues
262: Output Parameters:
263: . eigr - modified real part of eigenvalues
264: . eigi - modified imaginary part of eigenvalues
266: Notes:
267: This code implements the back transformation of eigenvalues in
268: order to retrieve the eigenvalues of the original problem. In this
269: example, simply set k_i = 1/k_i.
270: */
271: PetscErrorCode SampleShellSTBackTransform(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
272: {
273: PetscInt j;
276: for (j=0;j<n;j++) {
277: eigr[j] = 1.0 / eigr[j];
278: }
279: return(0);
280: }
281: /* ------------------------------------------------------------------- */
284: /*
285: SampleShellSTDestroy - This routine destroys a user-defined
286: spectral transformation context.
288: Input Parameter:
289: . shell - user-defined spectral transformation context
290: */
291: PetscErrorCode SampleShellSTDestroy(SampleShellST *shell)
292: {
296: KSPDestroy(&shell->ksp);
297: PetscFree(shell);
298: return(0);
299: }