Actual source code: ex10.c
slepc-3.17.1 2022-04-11
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*/