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: }