Actual source code: ex10.c

slepc-3.17.1 2022-04-11
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

120:   EPSSolve(eps);

122:   /*
123:      Optional: Get some information from the solver and display it
124:   */
125:   EPSGetType(eps,&type);
126:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
127:   EPSGetDimensions(eps,&nev,NULL,NULL);
128:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);

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

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

150: /***********************************************************************/
151: /*     Routines for a user-defined shell spectral transformation       */
152: /***********************************************************************/

154: /*
155:    STCreate_User - This routine creates a user-defined
156:    spectral transformation context.

158:    Output Parameter:
159: .  shell - user-defined spectral transformation context
160: */
161: PetscErrorCode STCreate_User(SampleShellST **shell)
162: {
163:   SampleShellST  *newctx;

166:   PetscNew(&newctx);
167:   KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
168:   KSPAppendOptionsPrefix(newctx->ksp,"st_");
169:   *shell = newctx;
170:   PetscFunctionReturn(0);
171: }
172: /* ------------------------------------------------------------------- */
173: /*
174:    STSetUp_User - This routine sets up a user-defined
175:    spectral transformation context.

177:    Input Parameters:
178: +  shell - user-defined spectral transformation context
179: -  st    - spectral transformation context containing the operator matrices

181:    Output Parameter:
182: .  shell - fully set up user-defined transformation context

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

195:   STGetMatrix(st,0,&A);
196:   KSPSetOperators(shell->ksp,A,A);
197:   KSPSetFromOptions(shell->ksp);
198:   PetscFunctionReturn(0);
199: }
200: /* ------------------------------------------------------------------- */
201: /*
202:    STApply_User - This routine demonstrates the use of a
203:    user-provided spectral transformation.

205:    Input Parameters:
206: +  st - spectral transformation context
207: -  x - input vector

209:    Output Parameter:
210: .  y - output vector

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

222:   STShellGetContext(st,&shell);
223:   KSPSolve(shell->ksp,x,y);
224:   PetscFunctionReturn(0);
225: }
226: /* ------------------------------------------------------------------- */
227: /*
228:    STApplyTranspose_User - This is not required unless using a two-sided
229:    eigensolver.

231:    Input Parameters:
232: +  st - spectral transformation context
233: -  x - input vector

235:    Output Parameter:
236: .  y - output vector
237: */
238: PetscErrorCode STApplyTranspose_User(ST st,Vec x,Vec y)
239: {
240:   SampleShellST  *shell;

243:   STShellGetContext(st,&shell);
244:   KSPSolveTranspose(shell->ksp,x,y);
245:   PetscFunctionReturn(0);
246: }
247: /* ------------------------------------------------------------------- */
248: /*
249:    STBackTransform_User - This routine demonstrates the use of a
250:    user-provided spectral transformation.

252:    Input Parameters:
253: +  st - spectral transformation context
254: -  n  - number of eigenvalues to transform

256:    Input/Output Parameters:
257: +  eigr - pointer to real part of eigenvalues
258: -  eigi - pointer to imaginary part of eigenvalues

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

270:   for (j=0;j<n;j++) {
271:     eigr[j] = 1.0 / eigr[j];
272:   }
273:   PetscFunctionReturn(0);
274: }
275: /* ------------------------------------------------------------------- */
276: /*
277:    STDestroy_User - This routine destroys a user-defined
278:    spectral transformation context.

280:    Input Parameter:
281: .  shell - user-defined spectral transformation context
282: */
283: PetscErrorCode STDestroy_User(SampleShellST *shell)
284: {
286:   KSPDestroy(&shell->ksp);
287:   PetscFree(shell);
288:   PetscFunctionReturn(0);
289: }

291: /*TEST

293:    testset:
294:       args: -eps_nev 5 -eps_non_hermitian -terse
295:       output_file: output/ex10_1.out
296:       test:
297:          suffix: 1_sinvert
298:          args: -st_type sinvert
299:       test:
300:          suffix: 1_sinvert_twoside
301:          args: -st_type sinvert -eps_balance twoside
302:          requires: !single
303:       test:
304:          suffix: 1_shell
305:          args: -st_type shell
306:          requires: !single
307:       test:
308:          suffix: 1_shell_twoside
309:          args: -st_type shell -eps_balance twoside

311: TEST*/