Actual source code: ex10.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[] = "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 STCreate_User(SampleShellST**);
37: PetscErrorCode STSetUp_User(SampleShellST*,ST);
38: PetscErrorCode STApply_User(ST,Vec,Vec);
39: PetscErrorCode STBackTransform_User(ST,PetscInt,PetscScalar*,PetscScalar*);
40: PetscErrorCode STDestroy_User(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: EPSType type;
51: PetscScalar value[3];
52: PetscInt n=30,i,col[3],Istart,Iend,FirstBlock=0,LastBlock=0,nev;
53: PetscBool isShell;
56: SlepcInitialize(&argc,&argv,(char*)0,help);
58: PetscOptionsGetInt(NULL,"-n",&n,NULL);
59: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%D\n\n",n);
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Compute the operator matrix that defines the eigensystem, Ax=kx
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: MatCreate(PETSC_COMM_WORLD,&A);
66: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
67: MatSetFromOptions(A);
68: MatSetUp(A);
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,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: PetscObjectTypeCompare((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: STCreate_User(&shell);
120: /* (Required) Set the user-defined routine for applying the operator */
121: STShellSetApply(st,STApply_User);
122: STShellSetContext(st,shell);
124: /* (Optional) Set the user-defined routine for back-transformation */
125: STShellSetBackTransform(st,STBackTransform_User);
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: STSetUp_User(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,NULL,NULL);
146: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Display solution and clean up
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: EPSPrintSolution(eps,NULL);
153: if (isShell) {
154: STDestroy_User(shell);
155: }
156: EPSDestroy(&eps);
157: MatDestroy(&A);
158: SlepcFinalize();
159: return 0;
160: }
162: /***********************************************************************/
163: /* Routines for a user-defined shell spectral transformation */
164: /***********************************************************************/
168: /*
169: STCreate_User - This routine creates a user-defined
170: spectral transformation context.
172: Output Parameter:
173: . shell - user-defined spectral transformation context
174: */
175: PetscErrorCode STCreate_User(SampleShellST **shell)
176: {
177: SampleShellST *newctx;
181: PetscNew(SampleShellST,&newctx);
182: KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
183: KSPAppendOptionsPrefix(newctx->ksp,"st_");
184: *shell = newctx;
185: return(0);
186: }
187: /* ------------------------------------------------------------------- */
190: /*
191: STSetUp_User - This routine sets up a user-defined
192: spectral transformation context.
194: Input Parameters:
195: . shell - user-defined spectral transformation context
196: . st - spectral transformation context containing the operator matrices
198: Output Parameter:
199: . shell - fully set up user-defined transformation context
201: Notes:
202: In this example, the user-defined transformation is simply OP=A^-1.
203: Therefore, the eigenpairs converge in reversed order. The KSP object
204: used for the solution of linear systems with A is handled via the
205: user-defined context SampleShellST.
206: */
207: PetscErrorCode STSetUp_User(SampleShellST *shell,ST st)
208: {
209: Mat A;
213: STGetOperators(st,0,&A);
214: KSPSetOperators(shell->ksp,A,A,DIFFERENT_NONZERO_PATTERN);
215: KSPSetFromOptions(shell->ksp);
216: return(0);
217: }
218: /* ------------------------------------------------------------------- */
221: /*
222: STApply_User - This routine demonstrates the use of a
223: user-provided spectral transformation.
225: Input Parameters:
226: . ctx - optional user-defined context, as set by STShellSetContext()
227: . x - input vector
229: Output Parameter:
230: . y - output vector
232: Notes:
233: The transformation implemented in this code is just OP=A^-1 and
234: therefore it is of little use, merely as an example of working with
235: a STSHELL.
236: */
237: PetscErrorCode STApply_User(ST st,Vec x,Vec y)
238: {
239: SampleShellST *shell;
243: STShellGetContext(st,(void**)&shell);
244: KSPSolve(shell->ksp,x,y);
245: return(0);
246: }
247: /* ------------------------------------------------------------------- */
250: /*
251: STBackTransform_User - This routine demonstrates the use of a
252: user-provided spectral transformation.
254: Input Parameters:
255: . ctx - optional user-defined context, as set by STShellSetContext()
256: . eigr - pointer to real part of eigenvalues
257: . eigi - pointer to imaginary part of eigenvalues
259: Output Parameters:
260: . eigr - modified real part of eigenvalues
261: . eigi - modified imaginary part of eigenvalues
263: Notes:
264: This code implements the back transformation of eigenvalues in
265: order to retrieve the eigenvalues of the original problem. In this
266: example, simply set k_i = 1/k_i.
267: */
268: PetscErrorCode STBackTransform_User(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
269: {
270: PetscInt j;
273: for (j=0;j<n;j++) {
274: eigr[j] = 1.0 / eigr[j];
275: }
276: return(0);
277: }
278: /* ------------------------------------------------------------------- */
281: /*
282: STDestroy_User - This routine destroys a user-defined
283: spectral transformation context.
285: Input Parameter:
286: . shell - user-defined spectral transformation context
287: */
288: PetscErrorCode STDestroy_User(SampleShellST *shell)
289: {
293: KSPDestroy(&shell->ksp);
294: PetscFree(shell);
295: return(0);
296: }