Actual source code: ex10.c
slepc-3.16.1 2021-11-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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;
43: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
45: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
46: PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%D\n\n",n);
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Compute the operator matrix that defines the eigensystem, Ax=kx
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: MatCreate(PETSC_COMM_WORLD,&A);
53: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);
54: MatSetFromOptions(A);
55: MatSetUp(A);
57: MatGetOwnershipRange(A,&Istart,&Iend);
58: for (i=Istart;i<Iend;i++) {
59: if (i>0) { MatSetValue(A,i,i-1,-1.0,INSERT_VALUES); }
60: if (i<n-1) { MatSetValue(A,i,i+1,-1.0,INSERT_VALUES); }
61: MatSetValue(A,i,i,2.0,INSERT_VALUES);
62: }
63: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
64: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create the eigensolver and set various options
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: /*
71: Create eigensolver context
72: */
73: EPSCreate(PETSC_COMM_WORLD,&eps);
75: /*
76: Set operators. In this case, it is a standard eigenvalue problem
77: */
78: EPSSetOperators(eps,A,NULL);
79: EPSSetProblemType(eps,EPS_HEP);
81: /*
82: Set solver parameters at runtime
83: */
84: EPSSetFromOptions(eps);
86: /*
87: Initialize shell spectral transformation if selected by user
88: */
89: EPSGetST(eps,&st);
90: PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell);
91: if (isShell) {
92: /* Change sorting criterion since this ST example computes values
93: closest to 0 */
94: 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: STCreate_User(&shell);
99: STShellSetContext(st,shell);
101: /* (Required) Set the user-defined routine for applying the operator */
102: STShellSetApply(st,STApply_User);
104: /* (Optional) Set the user-defined routine for applying the transposed operator */
105: STShellSetApplyTranspose(st,STApplyTranspose_User);
107: /* (Optional) Set the user-defined routine for back-transformation */
108: STShellSetBackTransform(st,STBackTransform_User);
110: /* (Optional) Set a name for the transformation, used for STView() */
111: PetscObjectSetName((PetscObject)st,"MyTransformation");
113: /* (Optional) Do any setup required for the new transformation */
114: STSetUp_User(shell,st);
115: }
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Solve the eigensystem
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: EPSSolve(eps);
123: /*
124: Optional: Get some information from the solver and display it
125: */
126: EPSGetType(eps,&type);
127: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
128: EPSGetDimensions(eps,&nev,NULL,NULL);
129: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Display solution and clean up
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135: /* show detailed info unless -terse option is given by user */
136: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
137: if (terse) {
138: EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
139: } else {
140: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
141: EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
142: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
143: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
144: }
145: if (isShell) {
146: STDestroy_User(shell);
147: }
148: EPSDestroy(&eps);
149: MatDestroy(&A);
150: SlepcFinalize();
151: return ierr;
152: }
154: /***********************************************************************/
155: /* Routines for a user-defined shell spectral transformation */
156: /***********************************************************************/
158: /*
159: STCreate_User - This routine creates a user-defined
160: spectral transformation context.
162: Output Parameter:
163: . shell - user-defined spectral transformation context
164: */
165: PetscErrorCode STCreate_User(SampleShellST **shell)
166: {
167: SampleShellST *newctx;
171: PetscNew(&newctx);
172: KSPCreate(PETSC_COMM_WORLD,&newctx->ksp);
173: KSPAppendOptionsPrefix(newctx->ksp,"st_");
174: *shell = newctx;
175: return(0);
176: }
177: /* ------------------------------------------------------------------- */
178: /*
179: STSetUp_User - This routine sets up a user-defined
180: spectral transformation context.
182: Input Parameters:
183: + shell - user-defined spectral transformation context
184: - st - spectral transformation context containing the operator matrices
186: Output Parameter:
187: . shell - fully set up user-defined transformation context
189: Notes:
190: In this example, the user-defined transformation is simply OP=A^-1.
191: Therefore, the eigenpairs converge in reversed order. The KSP object
192: used for the solution of linear systems with A is handled via the
193: user-defined context SampleShellST.
194: */
195: PetscErrorCode STSetUp_User(SampleShellST *shell,ST st)
196: {
197: Mat A;
201: STGetMatrix(st,0,&A);
202: KSPSetOperators(shell->ksp,A,A);
203: KSPSetFromOptions(shell->ksp);
204: return(0);
205: }
206: /* ------------------------------------------------------------------- */
207: /*
208: STApply_User - This routine demonstrates the use of a
209: user-provided spectral transformation.
211: Input Parameters:
212: + st - spectral transformation context
213: - x - input vector
215: Output Parameter:
216: . y - output vector
218: Notes:
219: The transformation implemented in this code is just OP=A^-1 and
220: therefore it is of little use, merely as an example of working with
221: a STSHELL.
222: */
223: PetscErrorCode STApply_User(ST st,Vec x,Vec y)
224: {
225: SampleShellST *shell;
229: STShellGetContext(st,&shell);
230: KSPSolve(shell->ksp,x,y);
231: return(0);
232: }
233: /* ------------------------------------------------------------------- */
234: /*
235: STApplyTranspose_User - This is not required unless using a two-sided
236: eigensolver.
238: Input Parameters:
239: + st - spectral transformation context
240: - x - input vector
242: Output Parameter:
243: . y - output vector
244: */
245: PetscErrorCode STApplyTranspose_User(ST st,Vec x,Vec y)
246: {
247: SampleShellST *shell;
251: STShellGetContext(st,&shell);
252: KSPSolveTranspose(shell->ksp,x,y);
253: return(0);
254: }
255: /* ------------------------------------------------------------------- */
256: /*
257: STBackTransform_User - This routine demonstrates the use of a
258: user-provided spectral transformation.
260: Input Parameters:
261: + st - spectral transformation context
262: - n - number of eigenvalues to transform
264: Input/Output Parameters:
265: + eigr - pointer to real part of eigenvalues
266: - eigi - pointer to imaginary part of eigenvalues
268: Notes:
269: This code implements the back transformation of eigenvalues in
270: order to retrieve the eigenvalues of the original problem. In this
271: example, simply set k_i = 1/k_i.
272: */
273: PetscErrorCode STBackTransform_User(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
274: {
275: PetscInt j;
278: for (j=0;j<n;j++) {
279: eigr[j] = 1.0 / eigr[j];
280: }
281: return(0);
282: }
283: /* ------------------------------------------------------------------- */
284: /*
285: STDestroy_User - This routine destroys a user-defined
286: spectral transformation context.
288: Input Parameter:
289: . shell - user-defined spectral transformation context
290: */
291: PetscErrorCode STDestroy_User(SampleShellST *shell)
292: {
296: KSPDestroy(&shell->ksp);
297: PetscFree(shell);
298: return(0);
299: }
301: /*TEST
303: testset:
304: args: -eps_nev 5 -eps_non_hermitian -terse
305: output_file: output/ex10_1.out
306: test:
307: suffix: 1_sinvert
308: args: -st_type sinvert
309: test:
310: suffix: 1_sinvert_twoside
311: args: -st_type sinvert -eps_balance twoside
312: requires: !single
313: test:
314: suffix: 1_shell
315: args: -st_type shell
316: requires: !single
317: test:
318: suffix: 1_shell_twoside
319: args: -st_type shell -eps_balance twoside
321: TEST*/