Actual source code: ex10.c

slepc-3.7.2 2016-07-19
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2016, 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:   PetscInt       n=30,i,Istart,Iend,nev;
 52:   PetscBool      isShell,terse;

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

 57:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 58:   PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%D\n\n",n);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:      Compute the operator matrix that defines the eigensystem, Ax=kx
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   MatCreate(PETSC_COMM_WORLD,&A);
 65:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 66:   MatSetFromOptions(A);
 67:   MatSetUp(A);

 69:   MatGetOwnershipRange(A,&Istart,&Iend);
 70:   for (i=Istart;i<Iend;i++) {
 71:     if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
 72:     if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
 73:     MatSetValue(A,i,i,2.0,INSERT_VALUES);
 74:   }
 75:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 76:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:                 Create the eigensolver and set various options
 80:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 82:   /*
 83:      Create eigensolver context
 84:   */
 85:   EPSCreate(PETSC_COMM_WORLD,&eps);

 87:   /*
 88:      Set operators. In this case, it is a standard eigenvalue problem
 89:   */
 90:   EPSSetOperators(eps,A,NULL);
 91:   EPSSetProblemType(eps,EPS_HEP);

 93:   /*
 94:      Set solver parameters at runtime
 95:   */
 96:   EPSSetFromOptions(eps);

 98:   /*
 99:      Initialize shell spectral transformation if selected by user
100:   */
101:   EPSGetST(eps,&st);
102:   PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell);
103:   if (isShell) {
104:     /* Change sorting criterion since this ST example computes values
105:        closest to 0 */
106:     EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);

108:     /* (Optional) Create a context for the user-defined spectral transform;
109:        this context can be defined to contain any application-specific data. */
110:     STCreate_User(&shell);

112:     /* (Required) Set the user-defined routine for applying the operator */
113:     STShellSetApply(st,STApply_User);
114:     STShellSetContext(st,shell);

116:     /* (Optional) Set the user-defined routine for back-transformation */
117:     STShellSetBackTransform(st,STBackTransform_User);

119:     /* (Optional) Set a name for the transformation, used for STView() */
120:     PetscObjectSetName((PetscObject)st,"MyTransformation");

122:     /* (Optional) Do any setup required for the new transformation */
123:     STSetUp_User(shell,st);
124:   }

126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:                       Solve the eigensystem
128:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

130:   EPSSolve(eps);

132:   /*
133:      Optional: Get some information from the solver and display it
134:   */
135:   EPSGetType(eps,&type);
136:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
137:   EPSGetDimensions(eps,&nev,NULL,NULL);
138:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:                     Display solution and clean up
142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

144:   /* show detailed info unless -terse option is given by user */
145:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
146:   if (terse) {
147:     EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
148:   } else {
149:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
150:     EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
151:     EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
152:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
153:   }
154:   if (isShell) {
155:     STDestroy_User(shell);
156:   }
157:   EPSDestroy(&eps);
158:   MatDestroy(&A);
159:   SlepcFinalize();
160:   return ierr;
161: }

163: /***********************************************************************/
164: /*     Routines for a user-defined shell spectral transformation       */
165: /***********************************************************************/

169: /*
170:    STCreate_User - This routine creates a user-defined
171:    spectral transformation context.

173:    Output Parameter:
174: .  shell - user-defined spectral transformation context
175: */
176: PetscErrorCode STCreate_User(SampleShellST **shell)
177: {
178:   SampleShellST  *newctx;

182:   PetscNew(&newctx);
183:   KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
184:   KSPAppendOptionsPrefix(newctx->ksp,"st_");
185:   *shell = newctx;
186:   return(0);
187: }
188: /* ------------------------------------------------------------------- */
191: /*
192:    STSetUp_User - This routine sets up a user-defined
193:    spectral transformation context.

195:    Input Parameters:
196: .  shell - user-defined spectral transformation context
197: .  st    - spectral transformation context containing the operator matrices

199:    Output Parameter:
200: .  shell - fully set up user-defined transformation context

202:    Notes:
203:    In this example, the user-defined transformation is simply OP=A^-1.
204:    Therefore, the eigenpairs converge in reversed order. The KSP object
205:    used for the solution of linear systems with A is handled via the
206:    user-defined context SampleShellST.
207: */
208: PetscErrorCode STSetUp_User(SampleShellST *shell,ST st)
209: {
210:   Mat            A;

214:   STGetOperators(st,0,&A);
215:   KSPSetOperators(shell->ksp,A,A);
216:   KSPSetFromOptions(shell->ksp);
217:   return(0);
218: }
219: /* ------------------------------------------------------------------- */
222: /*
223:    STApply_User - This routine demonstrates the use of a
224:    user-provided spectral transformation.

226:    Input Parameters:
227: .  ctx - optional user-defined context, as set by STShellSetContext()
228: .  x - input vector

230:    Output Parameter:
231: .  y - output vector

233:    Notes:
234:    The transformation implemented in this code is just OP=A^-1 and
235:    therefore it is of little use, merely as an example of working with
236:    a STSHELL.
237: */
238: PetscErrorCode STApply_User(ST st,Vec x,Vec y)
239: {
240:   SampleShellST  *shell;

244:   STShellGetContext(st,(void**)&shell);
245:   KSPSolve(shell->ksp,x,y);
246:   return(0);
247: }
248: /* ------------------------------------------------------------------- */
251: /*
252:    STBackTransform_User - This routine demonstrates the use of a
253:    user-provided spectral transformation.

255:    Input Parameters:
256: +  ctx  - optional user-defined context, as set by STShellSetContext()
257: .  eigr - pointer to real part of eigenvalues
258: -  eigi - pointer to imaginary part of eigenvalues

260:    Output Parameters:
261: +  eigr - modified real part of eigenvalues
262: -  eigi - modified imaginary part of eigenvalues

264:    Notes:
265:    This code implements the back transformation of eigenvalues in
266:    order to retrieve the eigenvalues of the original problem. In this
267:    example, simply set k_i = 1/k_i.
268: */
269: PetscErrorCode STBackTransform_User(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
270: {
271:   PetscInt j;

274:   for (j=0;j<n;j++) {
275:     eigr[j] = 1.0 / eigr[j];
276:   }
277:   return(0);
278: }
279: /* ------------------------------------------------------------------- */
282: /*
283:    STDestroy_User - This routine destroys a user-defined
284:    spectral transformation context.

286:    Input Parameter:
287: .  shell - user-defined spectral transformation context
288: */
289: PetscErrorCode STDestroy_User(SampleShellST *shell)
290: {

294:   KSPDestroy(&shell->ksp);
295:   PetscFree(shell);
296:   return(0);
297: }