Actual source code: ex10.c
slepc-3.21.0 2024-03-30
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,(char*)0,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*/