Actual source code: ex10.c

slepc-3.20.1 2023-11-27
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:   PetscFunctionBeginUser;
 43:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

121:   PetscCall(EPSSolve(eps));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

292: /*TEST

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

313: TEST*/