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 : #pragma once
12 :
13 : #include <slepceps.h>
14 : #include <slepc/private/slepcimpl.h>
15 :
16 : /* SUBMANSEC = EPS */
17 :
18 : SLEPC_EXTERN PetscBool EPSRegisterAllCalled;
19 : SLEPC_EXTERN PetscBool EPSMonitorRegisterAllCalled;
20 : SLEPC_EXTERN PetscErrorCode EPSRegisterAll(void);
21 : SLEPC_EXTERN PetscErrorCode EPSMonitorRegisterAll(void);
22 : SLEPC_EXTERN PetscLogEvent EPS_SetUp,EPS_Solve,EPS_CISS_SVD;
23 :
24 : typedef struct _EPSOps *EPSOps;
25 :
26 : struct _EPSOps {
27 : PetscErrorCode (*solve)(EPS);
28 : PetscErrorCode (*setup)(EPS);
29 : PetscErrorCode (*setupsort)(EPS);
30 : PetscErrorCode (*setfromoptions)(EPS,PetscOptionItems*);
31 : PetscErrorCode (*publishoptions)(EPS);
32 : PetscErrorCode (*destroy)(EPS);
33 : PetscErrorCode (*reset)(EPS);
34 : PetscErrorCode (*view)(EPS,PetscViewer);
35 : PetscErrorCode (*backtransform)(EPS);
36 : PetscErrorCode (*computevectors)(EPS);
37 : PetscErrorCode (*setdefaultst)(EPS);
38 : PetscErrorCode (*setdstype)(EPS);
39 : };
40 :
41 : /*
42 : Maximum number of monitors you can run with a single EPS
43 : */
44 : #define MAXEPSMONITORS 5
45 :
46 : /*
47 : The solution process goes through several states
48 : */
49 : typedef enum { EPS_STATE_INITIAL,
50 : EPS_STATE_SETUP,
51 : EPS_STATE_SOLVED,
52 : EPS_STATE_EIGENVECTORS } EPSStateType;
53 :
54 : /*
55 : To classify the different solvers into categories
56 : */
57 : typedef enum { EPS_CATEGORY_KRYLOV, /* Krylov solver: relies on STApply and STBackTransform (same as OTHER) */
58 : EPS_CATEGORY_PRECOND, /* Preconditioned solver: uses ST only to manage preconditioner */
59 : EPS_CATEGORY_CONTOUR, /* Contour integral: ST used to solve linear systems at integration points */
60 : EPS_CATEGORY_OTHER } EPSSolverType;
61 :
62 : /*
63 : To check for unsupported features at EPSSetUp_XXX()
64 : */
65 : typedef enum { EPS_FEATURE_BALANCE=1, /* balancing */
66 : EPS_FEATURE_ARBITRARY=2, /* arbitrary selection of eigepairs */
67 : EPS_FEATURE_REGION=4, /* nontrivial region for filtering */
68 : EPS_FEATURE_EXTRACTION=8, /* extraction technique different from Ritz */
69 : EPS_FEATURE_CONVERGENCE=16, /* convergence test selected by user */
70 : EPS_FEATURE_STOPPING=32, /* stopping test */
71 : EPS_FEATURE_TWOSIDED=64 /* two-sided variant */
72 : } EPSFeatureType;
73 :
74 : /*
75 : Defines the EPS data structure
76 : */
77 : struct _p_EPS {
78 : PETSCHEADER(struct _EPSOps);
79 : /*------------------------- User parameters ---------------------------*/
80 : PetscInt max_it; /* maximum number of iterations */
81 : PetscInt nev; /* number of eigenvalues to compute */
82 : PetscInt ncv; /* number of basis vectors */
83 : PetscInt mpd; /* maximum dimension of projected problem */
84 : PetscInt nini,ninil; /* number of initial vectors (negative means not copied yet) */
85 : PetscInt nds; /* number of basis vectors of deflation space */
86 : PetscScalar target; /* target value */
87 : PetscReal tol; /* tolerance */
88 : EPSConv conv; /* convergence test */
89 : EPSStop stop; /* stopping test */
90 : EPSWhich which; /* which part of the spectrum to be sought */
91 : PetscReal inta,intb; /* interval [a,b] for spectrum slicing */
92 : EPSProblemType problem_type; /* which kind of problem to be solved */
93 : EPSExtraction extraction; /* which kind of extraction to be applied */
94 : EPSBalance balance; /* the balancing method */
95 : PetscInt balance_its; /* number of iterations of the balancing method */
96 : PetscReal balance_cutoff; /* cutoff value for balancing */
97 : PetscBool trueres; /* whether the true residual norm must be computed */
98 : PetscBool trackall; /* whether all the residuals must be computed */
99 : PetscBool purify; /* whether eigenvectors need to be purified */
100 : PetscBool twosided; /* whether to compute left eigenvectors (two-sided solver) */
101 :
102 : /*-------------- User-provided functions and contexts -----------------*/
103 : PetscErrorCode (*converged)(EPS,PetscScalar,PetscScalar,PetscReal,PetscReal*,void*);
104 : PetscErrorCode (*convergeduser)(EPS,PetscScalar,PetscScalar,PetscReal,PetscReal*,void*);
105 : PetscErrorCode (*convergeddestroy)(void*);
106 : PetscErrorCode (*stopping)(EPS,PetscInt,PetscInt,PetscInt,PetscInt,EPSConvergedReason*,void*);
107 : PetscErrorCode (*stoppinguser)(EPS,PetscInt,PetscInt,PetscInt,PetscInt,EPSConvergedReason*,void*);
108 : PetscErrorCode (*stoppingdestroy)(void*);
109 : PetscErrorCode (*arbitrary)(PetscScalar,PetscScalar,Vec,Vec,PetscScalar*,PetscScalar*,void*);
110 : void *convergedctx;
111 : void *stoppingctx;
112 : void *arbitraryctx;
113 : PetscErrorCode (*monitor[MAXEPSMONITORS])(EPS,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,void*);
114 : PetscErrorCode (*monitordestroy[MAXEPSMONITORS])(void**);
115 : void *monitorcontext[MAXEPSMONITORS];
116 : PetscInt numbermonitors;
117 :
118 : /*----------------- Child objects and working data -------------------*/
119 : ST st; /* spectral transformation object */
120 : DS ds; /* direct solver object */
121 : BV V; /* set of basis vectors and computed eigenvectors */
122 : BV W; /* left basis vectors (if left eigenvectors requested) */
123 : RG rg; /* optional region for filtering */
124 : SlepcSC sc; /* sorting criterion data */
125 : Vec D; /* diagonal matrix for balancing */
126 : Vec *IS,*ISL; /* references to user-provided initial spaces */
127 : Vec *defl; /* references to user-provided deflation space */
128 : PetscScalar *eigr,*eigi; /* real and imaginary parts of eigenvalues */
129 : PetscReal *errest; /* error estimates */
130 : PetscScalar *rr,*ri; /* values computed by user's arbitrary selection function */
131 : PetscInt *perm; /* permutation for eigenvalue ordering */
132 : PetscInt nwork; /* number of work vectors */
133 : Vec *work; /* work vectors */
134 : void *data; /* placeholder for solver-specific stuff */
135 :
136 : /* ----------------------- Status variables --------------------------*/
137 : EPSStateType state; /* initial -> setup -> solved -> eigenvectors */
138 : EPSSolverType categ; /* solver category */
139 : PetscInt nconv; /* number of converged eigenvalues */
140 : PetscInt its; /* number of iterations so far computed */
141 : PetscInt n,nloc; /* problem dimensions (global, local) */
142 : PetscReal nrma,nrmb; /* computed matrix norms */
143 : PetscBool useds; /* whether the solver uses the DS object or not */
144 : PetscBool isgeneralized;
145 : PetscBool ispositive;
146 : PetscBool ishermitian;
147 : EPSConvergedReason reason;
148 : };
149 :
150 : /*
151 : Macros to test valid EPS arguments
152 : */
153 : #if !defined(PETSC_USE_DEBUG)
154 :
155 : #define EPSCheckSolved(h,arg) do {(void)(h);} while (0)
156 :
157 : #else
158 :
159 : #define EPSCheckSolved(h,arg) \
160 : do { \
161 : PetscCheck((h)->state>=EPS_STATE_SOLVED,PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSolve() first: Parameter #%d",arg); \
162 : } while (0)
163 :
164 : #endif
165 :
166 : /*
167 : Macros to check settings at EPSSetUp()
168 : */
169 :
170 : /* EPSCheckHermitianDefinite: the problem is HEP or GHEP */
171 : #define EPSCheckHermitianDefiniteCondition(eps,condition,msg) \
172 : do { \
173 : if (condition) { \
174 : PetscCheck((eps)->ishermitian,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for non-%s problems",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
175 : PetscCheck(!(eps)->isgeneralized || (eps)->ispositive,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires that the problem is %s-definite",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
176 : } \
177 : } while (0)
178 : #define EPSCheckHermitianDefinite(eps) EPSCheckHermitianDefiniteCondition(eps,PETSC_TRUE,"")
179 :
180 : /* EPSCheckHermitian: the problem is HEP, GHEP, or GHIEP */
181 : #define EPSCheckHermitianCondition(eps,condition,msg) \
182 : do { \
183 : if (condition) { \
184 : PetscCheck((eps)->ishermitian,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for non-%s problems",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
185 : } \
186 : } while (0)
187 : #define EPSCheckHermitian(eps) EPSCheckHermitianCondition(eps,PETSC_TRUE,"")
188 :
189 : /* EPSCheckDefinite: the problem is not GHIEP */
190 : #define EPSCheckDefiniteCondition(eps,condition,msg) \
191 : do { \
192 : if (condition) { \
193 : PetscCheck(!(eps)->isgeneralized || !(eps)->ishermitian || (eps)->ispositive,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for %s-indefinite problems",((PetscObject)(eps))->type_name,(msg),SLEPC_STRING_HERMITIAN); \
194 : } \
195 : } while (0)
196 : #define EPSCheckDefinite(eps) EPSCheckDefiniteCondition(eps,PETSC_TRUE,"")
197 :
198 : /* EPSCheckStandard: the problem is HEP or NHEP */
199 : #define EPSCheckStandardCondition(eps,condition,msg) \
200 : do { \
201 : if (condition) { \
202 : PetscCheck(!(eps)->isgeneralized,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for generalized problems",((PetscObject)(eps))->type_name,(msg)); \
203 : } \
204 : } while (0)
205 : #define EPSCheckStandard(eps) EPSCheckStandardCondition(eps,PETSC_TRUE,"")
206 :
207 : /* EPSCheckSinvert: shift-and-invert ST */
208 : #define EPSCheckSinvertCondition(eps,condition,msg) \
209 : do { \
210 : if (condition) { \
211 : PetscBool __flg; \
212 : PetscCall(PetscObjectTypeCompare((PetscObject)(eps)->st,STSINVERT,&__flg)); \
213 : PetscCheck(__flg,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires a shift-and-invert spectral transform",((PetscObject)(eps))->type_name,(msg)); \
214 : } \
215 : } while (0)
216 : #define EPSCheckSinvert(eps) EPSCheckSinvertCondition(eps,PETSC_TRUE,"")
217 :
218 : /* EPSCheckSinvertCayley: shift-and-invert or Cayley ST */
219 : #define EPSCheckSinvertCayleyCondition(eps,condition,msg) \
220 : do { \
221 : if (condition) { \
222 : PetscBool __flg; \
223 : PetscCall(PetscObjectTypeCompareAny((PetscObject)(eps)->st,&__flg,STSINVERT,STCAYLEY,"")); \
224 : PetscCheck(__flg,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires shift-and-invert or Cayley transform",((PetscObject)(eps))->type_name,(msg)); \
225 : } \
226 : } while (0)
227 : #define EPSCheckSinvertCayley(eps) EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE,"")
228 :
229 : /* Check for unsupported features */
230 : #define EPSCheckUnsupportedCondition(eps,mask,condition,msg) \
231 : do { \
232 : if (condition) { \
233 : PetscCheck(!((mask) & EPS_FEATURE_BALANCE) || (eps)->balance==EPS_BALANCE_NONE,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support balancing",((PetscObject)(eps))->type_name,(msg)); \
234 : PetscCheck(!((mask) & EPS_FEATURE_ARBITRARY) || !(eps)->arbitrary,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support arbitrary selection of eigenpairs",((PetscObject)(eps))->type_name,(msg)); \
235 : if ((mask) & EPS_FEATURE_REGION) { \
236 : PetscBool __istrivial; \
237 : PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
238 : PetscCheck(__istrivial,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support region filtering",((PetscObject)(eps))->type_name,(msg)); \
239 : } \
240 : PetscCheck(!((mask) & EPS_FEATURE_EXTRACTION) || (eps)->extraction==EPS_RITZ,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports Ritz extraction",((PetscObject)(eps))->type_name,(msg)); \
241 : PetscCheck(!((mask) & EPS_FEATURE_CONVERGENCE) || (eps)->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default convergence test",((PetscObject)(eps))->type_name,(msg)); \
242 : PetscCheck(!((mask) & EPS_FEATURE_STOPPING) || (eps)->stopping==EPSStoppingBasic,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default stopping test",((PetscObject)(eps))->type_name,(msg)); \
243 : PetscCheck(!((mask) & EPS_FEATURE_TWOSIDED) || !(eps)->twosided,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot compute left eigenvectors (no two-sided variant)",((PetscObject)(eps))->type_name,(msg)); \
244 : } \
245 : } while (0)
246 : #define EPSCheckUnsupported(eps,mask) EPSCheckUnsupportedCondition(eps,mask,PETSC_TRUE,"")
247 :
248 : /* Check for ignored features */
249 : #define EPSCheckIgnoredCondition(eps,mask,condition,msg) \
250 : do { \
251 : if (condition) { \
252 : if (((mask) & EPS_FEATURE_BALANCE) && (eps)->balance!=EPS_BALANCE_NONE) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the balancing settings\n",((PetscObject)(eps))->type_name,(msg))); \
253 : if (((mask) & EPS_FEATURE_ARBITRARY) && (eps)->arbitrary) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the settings for arbitrary selection of eigenpairs\n",((PetscObject)(eps))->type_name,(msg))); \
254 : if ((mask) & EPS_FEATURE_REGION) { \
255 : PetscBool __istrivial; \
256 : PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
257 : if (!__istrivial) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the specified region\n",((PetscObject)(eps))->type_name,(msg))); \
258 : } \
259 : if (((mask) & EPS_FEATURE_EXTRACTION) && (eps)->extraction!=EPS_RITZ) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the extraction settings\n",((PetscObject)(eps))->type_name,(msg))); \
260 : if (((mask) & EPS_FEATURE_CONVERGENCE) && (eps)->converged!=EPSConvergedRelative) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the convergence test settings\n",((PetscObject)(eps))->type_name,(msg))); \
261 : if (((mask) & EPS_FEATURE_STOPPING) && (eps)->stopping!=EPSStoppingBasic) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the stopping test settings\n",((PetscObject)(eps))->type_name,(msg))); \
262 : if (((mask) & EPS_FEATURE_TWOSIDED) && (eps)->twosided) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the two-sided flag\n",((PetscObject)(eps))->type_name,(msg))); \
263 : } \
264 : } while (0)
265 : #define EPSCheckIgnored(eps,mask) EPSCheckIgnoredCondition(eps,mask,PETSC_TRUE,"")
266 :
267 : /*
268 : EPS_SetInnerProduct - set B matrix for inner product if appropriate.
269 : */
270 786 : static inline PetscErrorCode EPS_SetInnerProduct(EPS eps)
271 : {
272 786 : Mat B;
273 :
274 786 : PetscFunctionBegin;
275 786 : if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
276 786 : if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
277 148 : PetscCall(STGetBilinearForm(eps->st,&B));
278 148 : PetscCall(BVSetMatrix(eps->V,B,PetscNot(eps->ispositive)));
279 148 : PetscCall(MatDestroy(&B));
280 638 : } else PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
281 786 : PetscFunctionReturn(PETSC_SUCCESS);
282 : }
283 :
284 : /*
285 : EPS_Purify - purify the first k vectors in the V basis
286 : */
287 89 : static inline PetscErrorCode EPS_Purify(EPS eps,PetscInt k)
288 : {
289 89 : PetscInt i;
290 89 : Vec v,z;
291 :
292 89 : PetscFunctionBegin;
293 89 : PetscCall(BVCreateVec(eps->V,&v));
294 1801 : for (i=0;i<k;i++) {
295 1712 : PetscCall(BVCopyVec(eps->V,i,v));
296 1712 : PetscCall(BVGetColumn(eps->V,i,&z));
297 1712 : PetscCall(STApply(eps->st,v,z));
298 1712 : PetscCall(BVRestoreColumn(eps->V,i,&z));
299 : }
300 89 : PetscCall(VecDestroy(&v));
301 89 : PetscFunctionReturn(PETSC_SUCCESS);
302 : }
303 :
304 : /*
305 : EPS_KSPSetOperators - Sets the KSP matrices, see also ST_KSPSetOperators()
306 : */
307 403 : static inline PetscErrorCode EPS_KSPSetOperators(KSP ksp,Mat A,Mat B)
308 : {
309 403 : const char *prefix;
310 :
311 403 : PetscFunctionBegin;
312 403 : PetscCall(KSPSetOperators(ksp,A,B));
313 403 : PetscCall(MatGetOptionsPrefix(B,&prefix));
314 403 : if (!prefix) {
315 : /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS)
316 : only applies if the Mat has no user-defined prefix */
317 403 : PetscCall(KSPGetOptionsPrefix(ksp,&prefix));
318 403 : PetscCall(MatSetOptionsPrefix(B,prefix));
319 : }
320 403 : PetscFunctionReturn(PETSC_SUCCESS);
321 : }
322 :
323 : SLEPC_INTERN PetscErrorCode EPSSetWhichEigenpairs_Default(EPS);
324 : SLEPC_INTERN PetscErrorCode EPSSetDimensions_Default(EPS,PetscInt,PetscInt*,PetscInt*);
325 : SLEPC_INTERN PetscErrorCode EPSBackTransform_Default(EPS);
326 : SLEPC_INTERN PetscErrorCode EPSComputeVectors(EPS);
327 : SLEPC_INTERN PetscErrorCode EPSComputeVectors_Hermitian(EPS);
328 : SLEPC_INTERN PetscErrorCode EPSComputeVectors_Schur(EPS);
329 : SLEPC_INTERN PetscErrorCode EPSComputeVectors_Indefinite(EPS);
330 : SLEPC_INTERN PetscErrorCode EPSComputeVectors_Twosided(EPS);
331 : SLEPC_INTERN PetscErrorCode EPSComputeVectors_Slice(EPS);
332 : SLEPC_INTERN PetscErrorCode EPSComputeResidualNorm_Private(EPS,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
333 : SLEPC_INTERN PetscErrorCode EPSComputeRitzVector(EPS,PetscScalar*,PetscScalar*,BV,Vec,Vec);
334 : SLEPC_INTERN PetscErrorCode EPSGetStartVector(EPS,PetscInt,PetscBool*);
335 : SLEPC_INTERN PetscErrorCode EPSGetLeftStartVector(EPS,PetscInt,PetscBool*);
336 : SLEPC_INTERN PetscErrorCode MatEstimateSpectralRange_EPS(Mat,PetscReal*,PetscReal*);
337 :
338 : /* Private functions of the solver implementations */
339 :
340 : SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
341 : SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi1(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
342 : SLEPC_INTERN PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscInt,PetscInt,PetscReal,PetscReal,PetscReal,PetscInt*);
343 : SLEPC_INTERN PetscErrorCode EPSPseudoLanczos(EPS,PetscReal*,PetscReal*,PetscReal*,PetscInt,PetscInt*,PetscBool*,PetscBool*,PetscReal*,Vec);
344 : SLEPC_INTERN PetscErrorCode EPSBuildBalance_Krylov(EPS);
345 : SLEPC_INTERN PetscErrorCode EPSSetDefaultST(EPS);
346 : SLEPC_INTERN PetscErrorCode EPSSetDefaultST_Precond(EPS);
347 : SLEPC_INTERN PetscErrorCode EPSSetDefaultST_GMRES(EPS);
348 : SLEPC_INTERN PetscErrorCode EPSSetDefaultST_NoFactor(EPS);
349 : SLEPC_INTERN PetscErrorCode EPSSetUpSort_Basic(EPS);
350 : SLEPC_INTERN PetscErrorCode EPSSetUpSort_Default(EPS);
|
