Actual source code: ex10.c

slepc-main 2024-11-15
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: #if defined(PETSC_USE_COMPLEX)
 30: PetscErrorCode STApplyHermitianTranspose_User(ST,Vec,Vec);
 31: #endif
 32: PetscErrorCode STBackTransform_User(ST,PetscInt,PetscScalar*,PetscScalar*);
 33: PetscErrorCode STDestroy_User(SampleShellST*);

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

 46:   PetscFunctionBeginUser;
 47:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

 49:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
 50:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%" PetscInt_FMT "\n\n",n));
 51:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-set_ht",&set_ht,NULL));

 53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:      Compute the operator matrix that defines the eigensystem, Ax=kx
 55:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 57:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
 58:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
 59:   PetscCall(MatSetFromOptions(A));

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

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

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

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

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

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

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

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

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

111:     /* (Optional) Set the user-defined routine for applying the conjugate-transposed operator */
112: #if defined(PETSC_USE_COMPLEX)
113:     if (set_ht) PetscCall(STShellSetApplyHermitianTranspose(st,STApplyHermitianTranspose_User));
114: #endif

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

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

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

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

130:   PetscCall(EPSSolve(eps));

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

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

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

160: /***********************************************************************/
161: /*     Routines for a user-defined shell spectral transformation       */
162: /***********************************************************************/

164: /*
165:    STCreate_User - This routine creates a user-defined
166:    spectral transformation context.

168:    Output Parameter:
169: .  shell - user-defined spectral transformation context
170: */
171: PetscErrorCode STCreate_User(SampleShellST **shell)
172: {
173:   SampleShellST  *newctx;

175:   PetscFunctionBeginUser;
176:   PetscCall(PetscNew(&newctx));
177:   PetscCall(KSPCreate(PETSC_COMM_WORLD,&newctx->ksp));
178:   PetscCall(KSPAppendOptionsPrefix(newctx->ksp,"st_"));
179:   *shell = newctx;
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }
182: /* ------------------------------------------------------------------- */
183: /*
184:    STSetUp_User - This routine sets up a user-defined
185:    spectral transformation context.

187:    Input Parameters:
188: +  shell - user-defined spectral transformation context
189: -  st    - spectral transformation context containing the operator matrices

191:    Output Parameter:
192: .  shell - fully set up user-defined transformation context

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

204:   PetscFunctionBeginUser;
205:   PetscCall(STGetMatrix(st,0,&A));
206:   PetscCall(KSPSetOperators(shell->ksp,A,A));
207:   PetscCall(KSPSetFromOptions(shell->ksp));
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }
210: /* ------------------------------------------------------------------- */
211: /*
212:    STApply_User - This routine demonstrates the use of a
213:    user-provided spectral transformation.

215:    Input Parameters:
216: +  st - spectral transformation context
217: -  x - input vector

219:    Output Parameter:
220: .  y - output vector

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

231:   PetscFunctionBeginUser;
232:   PetscCall(STShellGetContext(st,&shell));
233:   PetscCall(KSPSolve(shell->ksp,x,y));
234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }
236: /* ------------------------------------------------------------------- */
237: /*
238:    STApplyTranspose_User - This is not required unless using a two-sided
239:    eigensolver.

241:    Input Parameters:
242: +  st - spectral transformation context
243: -  x - input vector

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

252:   PetscFunctionBeginUser;
253:   PetscCall(STShellGetContext(st,&shell));
254:   PetscCall(KSPSolveTranspose(shell->ksp,x,y));
255:   PetscFunctionReturn(PETSC_SUCCESS);
256: }
257: #if defined(PETSC_USE_COMPLEX)
258: /* ------------------------------------------------------------------- */
259: /*
260:    STApplyHermitianTranspose_User - This is not required unless using a two-sided
261:    eigensolver in complex scalars.

263:    Input Parameters:
264: +  st - spectral transformation context
265: -  x - input vector

267:    Output Parameter:
268: .  y - output vector
269: */
270: PetscErrorCode STApplyHermitianTranspose_User(ST st,Vec x,Vec y)
271: {
272:   SampleShellST  *shell;
273:   Vec            w;

275:   PetscFunctionBeginUser;
276:   PetscCall(STShellGetContext(st,&shell));
277:   PetscCall(VecDuplicate(x,&w));
278:   PetscCall(VecCopy(x,w));
279:   PetscCall(VecConjugate(w));
280:   PetscCall(KSPSolveTranspose(shell->ksp,w,y));
281:   PetscCall(VecConjugate(y));
282:   PetscCall(VecDestroy(&w));
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }
285: #endif
286: /* ------------------------------------------------------------------- */
287: /*
288:    STBackTransform_User - This routine demonstrates the use of a
289:    user-provided spectral transformation.

291:    Input Parameters:
292: +  st - spectral transformation context
293: -  n  - number of eigenvalues to transform

295:    Input/Output Parameters:
296: +  eigr - pointer to real part of eigenvalues
297: -  eigi - pointer to imaginary part of eigenvalues

299:    Notes:
300:    This code implements the back transformation of eigenvalues in
301:    order to retrieve the eigenvalues of the original problem. In this
302:    example, simply set k_i = 1/k_i.
303: */
304: PetscErrorCode STBackTransform_User(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
305: {
306:   PetscInt j;

308:   PetscFunctionBeginUser;
309:   for (j=0;j<n;j++) {
310:     eigr[j] = 1.0 / eigr[j];
311:   }
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }
314: /* ------------------------------------------------------------------- */
315: /*
316:    STDestroy_User - This routine destroys a user-defined
317:    spectral transformation context.

319:    Input Parameter:
320: .  shell - user-defined spectral transformation context
321: */
322: PetscErrorCode STDestroy_User(SampleShellST *shell)
323: {
324:   PetscFunctionBeginUser;
325:   PetscCall(KSPDestroy(&shell->ksp));
326:   PetscCall(PetscFree(shell));
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

330: /*TEST

332:    testset:
333:       args: -eps_nev 5 -eps_non_hermitian -terse
334:       output_file: output/ex10_1.out
335:       test:
336:          suffix: 1_sinvert
337:          args: -st_type sinvert
338:          requires: !single
339:       test:
340:          suffix: 1_sinvert_twoside
341:          args: -st_type sinvert -eps_balance twoside -set_ht {{0 1}}
342:          requires: !single
343:       test:
344:          suffix: 1_shell
345:          args: -st_type shell
346:          requires: !single
347:       test:
348:          suffix: 1_shell_twoside
349:          args: -st_type shell -eps_balance twoside -set_ht {{0 1}}

351: TEST*/