Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
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";
16 :
17 : #include <slepceps.h>
18 :
19 : /* Define context for user-provided spectral transformation */
20 : typedef struct {
21 : KSP ksp;
22 : } SampleShellST;
23 :
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*);
34 :
35 6 : int main (int argc,char **argv)
36 : {
37 6 : Mat A; /* operator matrix */
38 6 : EPS eps; /* eigenproblem solver context */
39 6 : ST st; /* spectral transformation context */
40 6 : SampleShellST *shell; /* user-defined spectral transform context */
41 6 : EPSType type;
42 6 : PetscInt n=30,i,Istart,Iend,nev;
43 6 : PetscBool isShell,terse;
44 6 : PetscBool set_ht=PETSC_FALSE;
45 :
46 6 : PetscFunctionBeginUser;
47 6 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
48 :
49 6 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
50 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%" PetscInt_FMT "\n\n",n));
51 6 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-set_ht",&set_ht,NULL));
52 :
53 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54 : Compute the operator matrix that defines the eigensystem, Ax=kx
55 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56 :
57 6 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
58 6 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
59 6 : PetscCall(MatSetFromOptions(A));
60 :
61 6 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
62 186 : for (i=Istart;i<Iend;i++) {
63 180 : if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
64 180 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
65 180 : PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
66 : }
67 6 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
68 6 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
69 :
70 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 : Create the eigensolver and set various options
72 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73 :
74 : /*
75 : Create eigensolver context
76 : */
77 6 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
78 :
79 : /*
80 : Set operators. In this case, it is a standard eigenvalue problem
81 : */
82 6 : PetscCall(EPSSetOperators(eps,A,NULL));
83 6 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
84 :
85 : /*
86 : Set solver parameters at runtime
87 : */
88 6 : PetscCall(EPSSetFromOptions(eps));
89 :
90 : /*
91 : Initialize shell spectral transformation if selected by user
92 : */
93 6 : PetscCall(EPSGetST(eps,&st));
94 6 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell));
95 6 : if (isShell) {
96 : /* Change sorting criterion since this ST example computes values
97 : closest to 0 */
98 3 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
99 :
100 : /* (Required) Create a context for the user-defined spectral transform;
101 : this context can be defined to contain any application-specific data. */
102 3 : PetscCall(STCreate_User(&shell));
103 3 : PetscCall(STShellSetContext(st,shell));
104 :
105 : /* (Required) Set the user-defined routine for applying the operator */
106 3 : PetscCall(STShellSetApply(st,STApply_User));
107 :
108 : /* (Optional) Set the user-defined routine for applying the transposed operator */
109 3 : PetscCall(STShellSetApplyTranspose(st,STApplyTranspose_User));
110 :
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
115 :
116 : /* (Optional) Set the user-defined routine for back-transformation */
117 3 : PetscCall(STShellSetBackTransform(st,STBackTransform_User));
118 :
119 : /* (Optional) Set a name for the transformation, used for STView() */
120 3 : PetscCall(PetscObjectSetName((PetscObject)st,"MyTransformation"));
121 :
122 : /* (Optional) Do any setup required for the new transformation */
123 3 : PetscCall(STSetUp_User(shell,st));
124 : }
125 :
126 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 : Solve the eigensystem
128 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129 :
130 6 : PetscCall(EPSSolve(eps));
131 :
132 : /*
133 : Optional: Get some information from the solver and display it
134 : */
135 6 : PetscCall(EPSGetType(eps,&type));
136 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
137 6 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
138 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
139 :
140 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141 : Display solution and clean up
142 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143 :
144 : /* show detailed info unless -terse option is given by user */
145 6 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
146 6 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
147 : else {
148 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
149 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
150 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
151 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
152 : }
153 6 : if (isShell) PetscCall(STDestroy_User(shell));
154 6 : PetscCall(EPSDestroy(&eps));
155 6 : PetscCall(MatDestroy(&A));
156 6 : PetscCall(SlepcFinalize());
157 : return 0;
158 : }
159 :
160 : /***********************************************************************/
161 : /* Routines for a user-defined shell spectral transformation */
162 : /***********************************************************************/
163 :
164 : /*
165 : STCreate_User - This routine creates a user-defined
166 : spectral transformation context.
167 :
168 : Output Parameter:
169 : . shell - user-defined spectral transformation context
170 : */
171 3 : PetscErrorCode STCreate_User(SampleShellST **shell)
172 : {
173 3 : SampleShellST *newctx;
174 :
175 3 : PetscFunctionBeginUser;
176 3 : PetscCall(PetscNew(&newctx));
177 3 : PetscCall(KSPCreate(PETSC_COMM_WORLD,&newctx->ksp));
178 3 : PetscCall(KSPAppendOptionsPrefix(newctx->ksp,"st_"));
179 3 : *shell = newctx;
180 3 : PetscFunctionReturn(PETSC_SUCCESS);
181 : }
182 : /* ------------------------------------------------------------------- */
183 : /*
184 : STSetUp_User - This routine sets up a user-defined
185 : spectral transformation context.
186 :
187 : Input Parameters:
188 : + shell - user-defined spectral transformation context
189 : - st - spectral transformation context containing the operator matrices
190 :
191 : Output Parameter:
192 : . shell - fully set up user-defined transformation context
193 :
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 3 : PetscErrorCode STSetUp_User(SampleShellST *shell,ST st)
201 : {
202 3 : Mat A;
203 :
204 3 : PetscFunctionBeginUser;
205 3 : PetscCall(STGetMatrix(st,0,&A));
206 3 : PetscCall(KSPSetOperators(shell->ksp,A,A));
207 3 : PetscCall(KSPSetFromOptions(shell->ksp));
208 3 : PetscFunctionReturn(PETSC_SUCCESS);
209 : }
210 : /* ------------------------------------------------------------------- */
211 : /*
212 : STApply_User - This routine demonstrates the use of a
213 : user-provided spectral transformation.
214 :
215 : Input Parameters:
216 : + st - spectral transformation context
217 : - x - input vector
218 :
219 : Output Parameter:
220 : . y - output vector
221 :
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 70 : PetscErrorCode STApply_User(ST st,Vec x,Vec y)
228 : {
229 70 : SampleShellST *shell;
230 :
231 70 : PetscFunctionBeginUser;
232 70 : PetscCall(STShellGetContext(st,&shell));
233 70 : PetscCall(KSPSolve(shell->ksp,x,y));
234 70 : PetscFunctionReturn(PETSC_SUCCESS);
235 : }
236 : /* ------------------------------------------------------------------- */
237 : /*
238 : STApplyTranspose_User - This is not required unless using a two-sided
239 : eigensolver.
240 :
241 : Input Parameters:
242 : + st - spectral transformation context
243 : - x - input vector
244 :
245 : Output Parameter:
246 : . y - output vector
247 : */
248 10 : PetscErrorCode STApplyTranspose_User(ST st,Vec x,Vec y)
249 : {
250 10 : SampleShellST *shell;
251 :
252 10 : PetscFunctionBeginUser;
253 10 : PetscCall(STShellGetContext(st,&shell));
254 10 : PetscCall(KSPSolveTranspose(shell->ksp,x,y));
255 10 : 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.
262 :
263 : Input Parameters:
264 : + st - spectral transformation context
265 : - x - input vector
266 :
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;
274 :
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.
290 :
291 : Input Parameters:
292 : + st - spectral transformation context
293 : - n - number of eigenvalues to transform
294 :
295 : Input/Output Parameters:
296 : + eigr - pointer to real part of eigenvalues
297 : - eigi - pointer to imaginary part of eigenvalues
298 :
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 573 : PetscErrorCode STBackTransform_User(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
305 : {
306 573 : PetscInt j;
307 :
308 573 : PetscFunctionBeginUser;
309 1734 : for (j=0;j<n;j++) {
310 1161 : eigr[j] = 1.0 / eigr[j];
311 : }
312 573 : PetscFunctionReturn(PETSC_SUCCESS);
313 : }
314 : /* ------------------------------------------------------------------- */
315 : /*
316 : STDestroy_User - This routine destroys a user-defined
317 : spectral transformation context.
318 :
319 : Input Parameter:
320 : . shell - user-defined spectral transformation context
321 : */
322 3 : PetscErrorCode STDestroy_User(SampleShellST *shell)
323 : {
324 3 : PetscFunctionBeginUser;
325 3 : PetscCall(KSPDestroy(&shell->ksp));
326 3 : PetscCall(PetscFree(shell));
327 3 : PetscFunctionReturn(PETSC_SUCCESS);
328 : }
329 :
330 : /*TEST
331 :
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}}
350 :
351 : TEST*/
|
