Actual source code: ex10.c

  1: static char help[] = "Illustrates the use of shell spectral transformations. "
  2:   "The problem to be solved is the same as ex1.c and"
  3:   "corresponds to the Laplacian operator in 1 dimension.\n\n"
  4:   "The command line options are:\n"
  5:   "  -n <n>, where <n> = number of grid subdivisions = matrix dimension.\n\n";

 7:  #include slepceps.h

  9: /* Define context for user-provided spectral transformation */
 10: typedef struct {
 11:   KSP        ksp;
 12: } SampleShellST;

 14: /* Declare routines for user-provided spectral transformation */
 15: PetscErrorCode SampleShellSTCreate(SampleShellST**);
 16: PetscErrorCode SampleShellSTSetUp(SampleShellST*,ST);
 17: PetscErrorCode SampleShellSTApply(void*,Vec,Vec);
 18: PetscErrorCode SampleShellSTBackTransform(void*,PetscScalar*,PetscScalar*);
 19: PetscErrorCode SampleShellSTDestroy(SampleShellST*);

 23: int main( int argc, char **argv )
 24: {
 25:   Mat            A;                  /* operator matrix */
 26:   EPS            eps;                  /* eigenproblem solver context */
 27:   ST             st;                  /* spectral transformation context */
 28:   SampleShellST  *shell;          /* user-defined spectral transform context */
 29:   EPSType        type;
 30:   PetscReal      error, tol, re, im;
 31:   PetscScalar    kr, ki;
 33:   int            nev, maxit, its, nconv;
 34:   PetscInt       n=30, i, col[3], Istart, Iend, FirstBlock=0, LastBlock=0;
 35:   PetscScalar    value[3];
 36:   PetscTruth     isShell;

 38:   SlepcInitialize(&argc,&argv,(char*)0,help);

 40:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 41:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%d\n\n",n);

 43:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 44:      Compute the operator matrix that defines the eigensystem, Ax=kx
 45:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 47:   MatCreate(PETSC_COMM_WORLD,&A);
 48:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 49:   MatSetFromOptions(A);
 50: 
 51:   MatGetOwnershipRange(A,&Istart,&Iend);
 52:   if (Istart==0) FirstBlock=PETSC_TRUE;
 53:   if (Iend==n) LastBlock=PETSC_TRUE;
 54:   value[0]=-1.0; value[1]=2.0; value[2]=-1.0;
 55:   for( i=(FirstBlock? Istart+1: Istart); i<(LastBlock? Iend-1: Iend); i++ ) {
 56:     col[0]=i-1; col[1]=i; col[2]=i+1;
 57:     MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);
 58:   }
 59:   if (LastBlock) {
 60:     i=n-1; col[0]=n-2; col[1]=n-1;
 61:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
 62:   }
 63:   if (FirstBlock) {
 64:     i=0; col[0]=0; col[1]=1; value[0]=2.0; value[1]=-1.0;
 65:     MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);
 66:   }

 68:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 69:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 71:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 72:                 Create the eigensolver and set various options
 73:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 75:   /* 
 76:      Create eigensolver context
 77:   */
 78:   EPSCreate(PETSC_COMM_WORLD,&eps);

 80:   /* 
 81:      Set operators. In this case, it is a standard eigenvalue problem
 82:   */
 83:   EPSSetOperators(eps,A,PETSC_NULL);
 84:   EPSSetProblemType(eps,EPS_HEP);

 86:   /*
 87:      Set solver parameters at runtime
 88:   */
 89:   EPSSetFromOptions(eps);

 91:   /*
 92:      Initialize shell spectral transformation if selected by user
 93:   */
 94:   EPSGetST(eps,&st);
 95:   PetscTypeCompare((PetscObject)st,STSHELL,&isShell);
 96:   if (isShell) {
 97:     /* (Optional) Create a context for the user-defined spectral tranform;
 98:        this context can be defined to contain any application-specific data. */
 99:     SampleShellSTCreate(&shell);

101:     /* (Required) Set the user-defined routine for applying the operator */
102:     STShellSetApply(st,SampleShellSTApply);
103:     STShellSetContext(st,shell);

105:     /* (Optional) Set the user-defined routine for back-transformation */
106:     STShellSetBackTransform(st,SampleShellSTBackTransform);

108:     /* (Optional) Set a name for the transformation, used for STView() */
109:     STShellSetName(st,"MyTransformation");

111:     /* (Optional) Do any setup required for the new transformation */
112:     SampleShellSTSetUp(shell,st);
113:   }

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
116:                       Solve the eigensystem
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

119:   EPSSolve(eps);
120:   EPSGetIterationNumber(eps, &its);
121:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

123:   /*
124:      Optional: Get some information from the solver and display it
125:   */
126:   EPSGetType(eps,&type);
127:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
128:   EPSGetDimensions(eps,&nev,PETSC_NULL);
129:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
130:   EPSGetTolerances(eps,&tol,&maxit);
131:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
134:                     Display solution and clean up
135:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

137:   /* 
138:      Get number of converged approximate eigenpairs
139:   */
140:   EPSGetConverged(eps,&nconv);
141:   PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);

143:   if (nconv>0) {
144:     /*
145:        Display eigenvalues and relative errors
146:     */
147:     PetscPrintf(PETSC_COMM_WORLD,
148:          "           k          ||Ax-kx||/||kx||\n"
149:          "   ----------------- ------------------\n" );

151:     for( i=0; i<nconv; i++ ) {
152:       /* 
153:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
154:         ki (imaginary part)
155:       */
156:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
157:       /*
158:          Compute the relative error associated to each eigenpair
159:       */
160:       EPSComputeRelativeError(eps,i,&error);

162: #ifdef PETSC_USE_COMPLEX
163:       re = PetscRealPart(kr);
164:       im = PetscImaginaryPart(kr);
165: #else
166:       re = kr;
167:       im = ki;
168: #endif 
169:       if (im!=0.0) {
170:         PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
171:       } else {
172:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",re,error);
173:       }
174:     }
175:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
176:   }
177: 
178:   /* 
179:      Free work space
180:   */
181:   if (isShell) {
182:     SampleShellSTDestroy(shell);
183:   }
184:   EPSDestroy(eps);
185:   MatDestroy(A);
186:   SlepcFinalize();
187:   return 0;
188: }

190: /***********************************************************************/
191: /*     Routines for a user-defined shell spectral transformation       */
192: /***********************************************************************/

196: /*
197:    SampleShellSTCreate - This routine creates a user-defined
198:    spectral transformation context.

200:    Output Parameter:
201: .  shell - user-defined spectral transformation context
202: */
203: PetscErrorCode SampleShellSTCreate(SampleShellST **shell)
204: {
205:   SampleShellST  *newctx;

209:   PetscNew(SampleShellST,&newctx);
210:   KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
211:   KSPAppendOptionsPrefix(newctx->ksp,"st_");
212:   *shell = newctx;
213:   return(0);
214: }
215: /* ------------------------------------------------------------------- */
218: /*
219:    SampleShellSTSetUp - This routine sets up a user-defined
220:    spectral transformation context.  

222:    Input Parameters:
223: .  shell - user-defined spectral transformation context
224: .  st    - spectral transformation context containing the operator matrices

226:    Output Parameter:
227: .  shell - fully set up user-defined transformation context

229:    Notes:
230:    In this example, the user-defined transformation is simply OP=A^-1.
231:    Therefore, the eigenpairs converge in reversed order. The KSP object
232:    used for the solution of linear systems with A is handled via the
233:    user-defined context SampleShellST.
234: */
235: PetscErrorCode SampleShellSTSetUp(SampleShellST *shell,ST st)
236: {
237:   Mat            A,B;

241:   STGetOperators(st,&A,&B);
242:   if (B) { SETERRQ(0,"Warning: This transformation is not intended for generalized problems"); }
243:   KSPSetOperators(shell->ksp,A,A,DIFFERENT_NONZERO_PATTERN);
244:   KSPSetFromOptions(shell->ksp);
245:   return(0);
246: }
247: /* ------------------------------------------------------------------- */
250: /*
251:    SampleShellSTApply - 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: .  x - input vector

258:    Output Parameter:
259: .  y - output vector

261:    Notes:
262:    The transformation implemented in this code is just OP=A^-1 and
263:    therefore it is of little use, merely as an example of working with
264:    a STSHELL.
265: */
266: PetscErrorCode SampleShellSTApply(void *ctx,Vec x,Vec y)
267: {
268:   SampleShellST  *shell = (SampleShellST*)ctx;

272:   KSPSolve(shell->ksp,x,y);
273:   return(0);
274: }
275: /* ------------------------------------------------------------------- */
278: /*
279:    SampleShellSTBackTransform - This routine demonstrates the use of a
280:    user-provided spectral transformation.

282:    Input Parameters:
283: .  ctx  - optional user-defined context, as set by STShellSetContext()
284: .  eigr - pointer to real part of eigenvalues
285: .  eigi - pointer to imaginary part of eigenvalues

287:    Output Parameters:
288: .  eigr - modified real part of eigenvalues
289: .  eigi - modified imaginary part of eigenvalues

291:    Notes:
292:    This code implements the back transformation of eigenvalues in
293:    order to retrieve the eigenvalues of the original problem. In this
294:    example, simply set k_i = 1/k_i.
295: */
296: PetscErrorCode SampleShellSTBackTransform(void *ctx,PetscScalar *eigr,PetscScalar *eigi)
297: {
299:   *eigr = 1.0 / *eigr;
300:   return(0);
301: }
302: /* ------------------------------------------------------------------- */
305: /*
306:    SampleShellSTDestroy - This routine destroys a user-defined
307:    spectral transformation context.

309:    Input Parameter:
310: .  shell - user-defined spectral transformation context
311: */
312: PetscErrorCode SampleShellSTDestroy(SampleShellST *shell)
313: {

317:   KSPDestroy(shell->ksp);
318:   PetscFree(shell);
319:   return(0);
320: }