Actual source code: ex10.c
slepc-3.19.2 2023-09-05
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: test:
301: suffix: 1_sinvert_twoside
302: args: -st_type sinvert -eps_balance twoside
303: requires: !single
304: test:
305: suffix: 1_shell
306: args: -st_type shell
307: requires: !single
308: test:
309: suffix: 1_shell_twoside
310: args: -st_type shell -eps_balance twoside
312: TEST*/