Actual source code: ex10.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

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

 17: #include <slepceps.h>

 19: /* Define context for user-provided spectral transformation */
 20: typedef struct {
 21:   KSP        ksp;
 22: } SampleShellST;

 24: /* Declare routines for user-provided spectral transformation */
 25: PetscErrorCode STCreate_User(SampleShellST**);
 26: PetscErrorCode STSetUp_User(SampleShellST*,ST);
 27: PetscErrorCode STApply_User(ST,Vec,Vec);
 28: PetscErrorCode STApplyTranspose_User(ST,Vec,Vec);
 29: PetscErrorCode STBackTransform_User(ST,PetscInt,PetscScalar*,PetscScalar*);
 30: PetscErrorCode STDestroy_User(SampleShellST*);

 32: int main (int argc,char **argv)
 33: {
 34:   Mat            A;               /* operator matrix */
 35:   EPS            eps;             /* eigenproblem solver context */
 36:   ST             st;              /* spectral transformation context */
 37:   SampleShellST  *shell;          /* user-defined spectral transform context */
 38:   EPSType        type;
 39:   PetscInt       n=30,i,Istart,Iend,nev;
 40:   PetscBool      isShell,terse;

 43:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

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

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:      Compute the operator matrix that defines the eigensystem, Ax=kx
 50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 52:   MatCreate(PETSC_COMM_WORLD,&A);
 53:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
 54:   MatSetFromOptions(A);
 55:   MatSetUp(A);

 57:   MatGetOwnershipRange(A,&Istart,&Iend);
 58:   for (i=Istart;i<Iend;i++) {
 59:     if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
 60:     if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
 61:     MatSetValue(A,i,i,2.0,INSERT_VALUES);
 62:   }
 63:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 64:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 67:                 Create the eigensolver and set various options
 68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 70:   /*
 71:      Create eigensolver context
 72:   */
 73:   EPSCreate(PETSC_COMM_WORLD,&eps);

 75:   /*
 76:      Set operators. In this case, it is a standard eigenvalue problem
 77:   */
 78:   EPSSetOperators(eps,A,NULL);
 79:   EPSSetProblemType(eps,EPS_HEP);

 81:   /*
 82:      Set solver parameters at runtime
 83:   */
 84:   EPSSetFromOptions(eps);

 86:   /*
 87:      Initialize shell spectral transformation if selected by user
 88:   */
 89:   EPSGetST(eps,&st);
 90:   PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell);
 91:   if (isShell) {
 92:     /* Change sorting criterion since this ST example computes values
 93:        closest to 0 */
 94:     EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);

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

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

104:     /* (Optional) Set the user-defined routine for applying the transposed operator */
105:     STShellSetApplyTranspose(st,STApplyTranspose_User);

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

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

113:     /* (Optional) Do any setup required for the new transformation */
114:     STSetUp_User(shell,st);
115:   }

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:                       Solve the eigensystem
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   EPSSolve(eps);

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,NULL,NULL);
129:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132:                     Display solution and clean up
133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

135:   /* show detailed info unless -terse option is given by user */
136:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
137:   if (terse) {
138:     EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
139:   } else {
140:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
141:     EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
142:     EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
143:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
144:   }
145:   if (isShell) {
146:     STDestroy_User(shell);
147:   }
148:   EPSDestroy(&eps);
149:   MatDestroy(&A);
150:   SlepcFinalize();
151:   return ierr;
152: }

154: /***********************************************************************/
155: /*     Routines for a user-defined shell spectral transformation       */
156: /***********************************************************************/

158: /*
159:    STCreate_User - This routine creates a user-defined
160:    spectral transformation context.

162:    Output Parameter:
163: .  shell - user-defined spectral transformation context
164: */
165: PetscErrorCode STCreate_User(SampleShellST **shell)
166: {
167:   SampleShellST  *newctx;

171:   PetscNew(&newctx);
172:   KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
173:   KSPAppendOptionsPrefix(newctx->ksp,"st_");
174:   *shell = newctx;
175:   return(0);
176: }
177: /* ------------------------------------------------------------------- */
178: /*
179:    STSetUp_User - This routine sets up a user-defined
180:    spectral transformation context.

182:    Input Parameters:
183: +  shell - user-defined spectral transformation context
184: -  st    - spectral transformation context containing the operator matrices

186:    Output Parameter:
187: .  shell - fully set up user-defined transformation context

189:    Notes:
190:    In this example, the user-defined transformation is simply OP=A^-1.
191:    Therefore, the eigenpairs converge in reversed order. The KSP object
192:    used for the solution of linear systems with A is handled via the
193:    user-defined context SampleShellST.
194: */
195: PetscErrorCode STSetUp_User(SampleShellST *shell,ST st)
196: {
197:   Mat            A;

201:   STGetMatrix(st,0,&A);
202:   KSPSetOperators(shell->ksp,A,A);
203:   KSPSetFromOptions(shell->ksp);
204:   return(0);
205: }
206: /* ------------------------------------------------------------------- */
207: /*
208:    STApply_User - This routine demonstrates the use of a
209:    user-provided spectral transformation.

211:    Input Parameters:
212: +  st - spectral transformation context
213: -  x - input vector

215:    Output Parameter:
216: .  y - output vector

218:    Notes:
219:    The transformation implemented in this code is just OP=A^-1 and
220:    therefore it is of little use, merely as an example of working with
221:    a STSHELL.
222: */
223: PetscErrorCode STApply_User(ST st,Vec x,Vec y)
224: {
225:   SampleShellST  *shell;

229:   STShellGetContext(st,(void**)&shell);
230:   KSPSolve(shell->ksp,x,y);
231:   return(0);
232: }
233: /* ------------------------------------------------------------------- */
234: /*
235:    STApplyTranspose_User - This is not required unless using a two-sided
236:    eigensolver.

238:    Input Parameters:
239: +  st - spectral transformation context
240: -  x - input vector

242:    Output Parameter:
243: .  y - output vector
244: */
245: PetscErrorCode STApplyTranspose_User(ST st,Vec x,Vec y)
246: {
247:   SampleShellST  *shell;

251:   STShellGetContext(st,(void**)&shell);
252:   KSPSolveTranspose(shell->ksp,x,y);
253:   return(0);
254: }
255: /* ------------------------------------------------------------------- */
256: /*
257:    STBackTransform_User - This routine demonstrates the use of a
258:    user-provided spectral transformation.

260:    Input Parameters:
261: +  st - spectral transformation context
262: -  n  - number of eigenvalues to transform

264:    Input/Output Parameters:
265: +  eigr - pointer to real part of eigenvalues
266: -  eigi - pointer to imaginary part of eigenvalues

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

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

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

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

301: /*TEST

303:    testset:
304:       args: -eps_nev 5 -eps_non_hermitian -terse
305:       requires: !single
306:       output_file: output/ex10_1.out
307:       test:
308:          suffix: 1_sinvert
309:          args: -st_type sinvert
310:       test:
311:          suffix: 1_sinvert_twoside
312:          args: -st_type sinvert -eps_balance twoside
313:       test:
314:          suffix: 1_shell
315:          args: -st_type shell
316:       test:
317:          suffix: 1_shell_twoside
318:          args: -st_type shell -eps_balance twoside

320: TEST*/