Actual source code: ex10.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2010, 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(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      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(ST st,Vec x,Vec y)
287: {
288:   SampleShellST  *shell;

292:   STShellGetContext(st,(void**)&shell);
293:   KSPSolve(shell->ksp,x,y);
294:   return(0);
295: }
296: /* ------------------------------------------------------------------- */
299: /*
300:    SampleShellSTBackTransform - This routine demonstrates the use of a
301:    user-provided spectral transformation.

303:    Input Parameters:
304: .  ctx  - optional user-defined context, as set by STShellSetContext()
305: .  eigr - pointer to real part of eigenvalues
306: .  eigi - pointer to imaginary part of eigenvalues

308:    Output Parameters:
309: .  eigr - modified real part of eigenvalues
310: .  eigi - modified imaginary part of eigenvalues

312:    Notes:
313:    This code implements the back transformation of eigenvalues in
314:    order to retrieve the eigenvalues of the original problem. In this
315:    example, simply set k_i = 1/k_i.
316: */
317: PetscErrorCode SampleShellSTBackTransform(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
318: {
319:   PetscInt j;
321:   for (j=0;j<n;j++) {
322:     eigr[j] = 1.0 / eigr[j];
323:   }
324:   return(0);
325: }
326: /* ------------------------------------------------------------------- */
329: /*
330:    SampleShellSTDestroy - This routine destroys a user-defined
331:    spectral transformation context.

333:    Input Parameter:
334: .  shell - user-defined spectral transformation context
335: */
336: PetscErrorCode SampleShellSTDestroy(SampleShellST *shell)
337: {

341:   KSPDestroy(shell->ksp);
342:   PetscFree(shell);
343:   return(0);
344: }