Actual source code: epsimpl.h
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: #pragma once
13: #include <slepceps.h>
14: #include <slepc/private/bvimpl.h>
16: /* SUBMANSEC = EPS */
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;
24: typedef struct _EPSOps *EPSOps;
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: };
41: /*
42: Maximum number of monitors you can run with a single EPS
43: */
44: #define MAXEPSMONITORS 5
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;
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;
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_THRESHOLD=64, /* threshold stopping test */
72: EPS_FEATURE_TWOSIDED=128 /* two-sided variant */
73: } EPSFeatureType;
75: /*
76: Defines the EPS data structure
77: */
78: struct _p_EPS {
79: PETSCHEADER(struct _EPSOps);
80: /*------------------------- User parameters ---------------------------*/
81: PetscInt max_it; /* maximum number of iterations */
82: PetscInt nev; /* number of eigenvalues to compute */
83: PetscInt ncv; /* number of basis vectors */
84: PetscInt mpd; /* maximum dimension of projected problem */
85: PetscInt nini,ninil; /* number of initial vectors (negative means not copied yet) */
86: PetscInt nds; /* number of basis vectors of deflation space */
87: PetscScalar target; /* target value */
88: PetscReal tol; /* tolerance */
89: PetscReal thres; /* threshold */
90: PetscBool threlative; /* threshold is relative */
91: EPSConv conv; /* convergence test */
92: EPSStop stop; /* stopping test */
93: EPSWhich which; /* which part of the spectrum to be sought */
94: PetscReal inta,intb; /* interval [a,b] for spectrum slicing */
95: EPSProblemType problem_type; /* which kind of problem to be solved */
96: EPSExtraction extraction; /* which kind of extraction to be applied */
97: EPSBalance balance; /* the balancing method */
98: PetscInt balance_its; /* number of iterations of the balancing method */
99: PetscReal balance_cutoff; /* cutoff value for balancing */
100: PetscBool trueres; /* whether the true residual norm must be computed */
101: PetscBool trackall; /* whether all the residuals must be computed */
102: PetscBool purify; /* whether eigenvectors need to be purified */
103: PetscBool twosided; /* whether to compute left eigenvectors (two-sided solver) */
105: /*-------------- User-provided functions and contexts -----------------*/
106: EPSConvergenceTestFn *converged;
107: EPSConvergenceTestFn *convergeduser;
108: PetscCtxDestroyFn *convergeddestroy;
109: EPSStoppingTestFn *stopping;
110: EPSStoppingTestFn *stoppinguser;
111: PetscCtxDestroyFn *stoppingdestroy;
112: SlepcArbitrarySelectionFn *arbitrary;
113: PetscCtxDestroyFn *arbitrarydestroy;
114: void *convergedctx;
115: void *stoppingctx;
116: void *arbitraryctx;
117: EPSMonitorFn *monitor[MAXEPSMONITORS];
118: PetscCtxDestroyFn *monitordestroy[MAXEPSMONITORS];
119: void *monitorcontext[MAXEPSMONITORS];
120: PetscInt numbermonitors;
122: /*----------------- Child objects and working data -------------------*/
123: ST st; /* spectral transformation object */
124: DS ds; /* direct solver object */
125: BV V; /* set of basis vectors and computed eigenvectors */
126: BV W; /* left basis vectors (if left eigenvectors requested) */
127: RG rg; /* optional region for filtering */
128: SlepcSC sc; /* sorting criterion data */
129: Vec D; /* diagonal matrix for balancing */
130: Vec *IS,*ISL; /* references to user-provided initial spaces */
131: Vec *defl; /* references to user-provided deflation space */
132: PetscScalar *eigr,*eigi; /* real and imaginary parts of eigenvalues */
133: PetscReal *errest; /* error estimates */
134: PetscScalar *rr,*ri; /* values computed by user's arbitrary selection function */
135: PetscInt *perm; /* permutation for eigenvalue ordering */
136: PetscInt nwork; /* number of work vectors */
137: Vec *work; /* work vectors */
138: void *data; /* placeholder for solver-specific stuff */
140: /* ----------------------- Status variables --------------------------*/
141: EPSStateType state; /* initial -> setup -> solved -> eigenvectors */
142: EPSSolverType categ; /* solver category */
143: PetscInt nconv; /* number of converged eigenvalues */
144: PetscInt its; /* number of iterations so far computed */
145: PetscInt n,nloc; /* problem dimensions (global, local) */
146: PetscReal nrma,nrmb; /* computed matrix norms */
147: PetscBool useds; /* whether the solver uses the DS object or not */
148: PetscBool isgeneralized;
149: PetscBool ispositive;
150: PetscBool ishermitian;
151: PetscBool isstructured;
152: PetscInt setfromoptionscalled;
153: EPSConvergedReason reason;
154: };
156: /*
157: Macros to test valid EPS arguments
158: */
159: #if !PetscDefined(USE_DEBUG)
161: #define EPSCheckSolved(h,arg) do {(void)(h);} while (0)
163: #else
165: #define EPSCheckSolved(h,arg) \
166: do { \
167: PetscCheck((h)->state>=EPS_STATE_SOLVED,PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSolve() first: Parameter #%d",arg); \
168: } while (0)
170: #endif
172: /*
173: Macros to check settings at EPSSetUp()
174: */
176: /* EPSCheckHermitianDefinite: the problem is HEP or GHEP */
177: #define EPSCheckHermitianDefiniteCondition(eps,condition,msg) \
178: do { \
179: if (condition) { \
180: 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); \
181: 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); \
182: } \
183: } while (0)
184: #define EPSCheckHermitianDefinite(eps) EPSCheckHermitianDefiniteCondition(eps,PETSC_TRUE,"")
186: /* EPSCheckHermitian: the problem is HEP, GHEP, or GHIEP */
187: #define EPSCheckHermitianCondition(eps,condition,msg) \
188: do { \
189: if (condition) { \
190: 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); \
191: } \
192: } while (0)
193: #define EPSCheckHermitian(eps) EPSCheckHermitianCondition(eps,PETSC_TRUE,"")
195: /* EPSCheckDefinite: the problem is not GHIEP */
196: #define EPSCheckDefiniteCondition(eps,condition,msg) \
197: do { \
198: if (condition) { \
199: 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); \
200: } \
201: } while (0)
202: #define EPSCheckDefinite(eps) EPSCheckDefiniteCondition(eps,PETSC_TRUE,"")
204: /* EPSCheckStandard: the problem is HEP or NHEP */
205: #define EPSCheckStandardCondition(eps,condition,msg) \
206: do { \
207: if (condition) { \
208: PetscCheck(!(eps)->isgeneralized,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for generalized problems",((PetscObject)(eps))->type_name,(msg)); \
209: } \
210: } while (0)
211: #define EPSCheckStandard(eps) EPSCheckStandardCondition(eps,PETSC_TRUE,"")
213: /* EPSCheckNotStructured: the problem is not structured */
214: #define EPSCheckNotStructuredCondition(eps,condition,msg) \
215: do { \
216: if (condition) { \
217: PetscCheck(!(eps)->isstructured,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not provide support for structured eigenproblems",((PetscObject)(eps))->type_name,(msg)); \
218: } \
219: } while (0)
220: #define EPSCheckNotStructured(eps) EPSCheckNotStructuredCondition(eps,PETSC_TRUE,"")
222: /* EPSCheckSinvert: shift-and-invert ST */
223: #define EPSCheckSinvertCondition(eps,condition,msg) \
224: do { \
225: if (condition) { \
226: PetscBool __flg; \
227: PetscCall(PetscObjectTypeCompare((PetscObject)(eps)->st,STSINVERT,&__flg)); \
228: PetscCheck(__flg,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires a shift-and-invert spectral transform",((PetscObject)(eps))->type_name,(msg)); \
229: } \
230: } while (0)
231: #define EPSCheckSinvert(eps) EPSCheckSinvertCondition(eps,PETSC_TRUE,"")
233: /* EPSCheckSinvertCayley: shift-and-invert or Cayley ST */
234: #define EPSCheckSinvertCayleyCondition(eps,condition,msg) \
235: do { \
236: if (condition) { \
237: PetscBool __flg; \
238: PetscCall(PetscObjectTypeCompareAny((PetscObject)(eps)->st,&__flg,STSINVERT,STCAYLEY,"")); \
239: PetscCheck(__flg,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires shift-and-invert or Cayley transform",((PetscObject)(eps))->type_name,(msg)); \
240: } \
241: } while (0)
242: #define EPSCheckSinvertCayley(eps) EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE,"")
244: /* Check for unsupported features */
245: #define EPSCheckUnsupportedCondition(eps,mask,condition,msg) \
246: do { \
247: if (condition) { \
248: 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)); \
249: 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)); \
250: if ((mask) & EPS_FEATURE_REGION) { \
251: PetscBool __istrivial; \
252: PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
253: PetscCheck(__istrivial,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support region filtering",((PetscObject)(eps))->type_name,(msg)); \
254: } \
255: 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)); \
256: 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)); \
257: 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)); \
258: PetscCheck(!((mask) & EPS_FEATURE_THRESHOLD) || (eps)->stopping!=EPSStoppingThreshold,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support the threshold stopping test",((PetscObject)(eps))->type_name,(msg)); \
259: 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)); \
260: } \
261: } while (0)
262: #define EPSCheckUnsupported(eps,mask) EPSCheckUnsupportedCondition(eps,mask,PETSC_TRUE,"")
264: /* Check for ignored features */
265: #define EPSCheckIgnoredCondition(eps,mask,condition,msg) \
266: do { \
267: if (condition) { \
268: 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))); \
269: 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))); \
270: if ((mask) & EPS_FEATURE_REGION) { \
271: PetscBool __istrivial; \
272: PetscCall(RGIsTrivial((eps)->rg,&__istrivial)); \
273: if (!__istrivial) PetscCall(PetscInfo((eps),"The solver '%s'%s ignores the specified region\n",((PetscObject)(eps))->type_name,(msg))); \
274: } \
275: 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))); \
276: 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))); \
277: 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))); \
278: 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))); \
279: } \
280: } while (0)
281: #define EPSCheckIgnored(eps,mask) EPSCheckIgnoredCondition(eps,mask,PETSC_TRUE,"")
283: /*
284: EPSSetCtxThreshold - Fills EPSStoppingCtx with data needed for the threshold stopping test
286: k = number of converged approximations, n = total number of available approximations
287: */
288: #define EPSSetCtxThreshold(eps,eigr,eigi,err_est,k,n) \
289: do { \
290: if ((eps)->stop==EPS_STOP_THRESHOLD && (k)) { \
291: PetscScalar __kr=(eigr)[(k)-1],__ki=(eigi)[(k)-1],__krn=0.0,__kin=0.0,__kr0=(eigr)[0],__ki0=(eigi)[0]; \
292: PetscCall(STBackTransform((eps)->st,1,&__kr,&__ki)); \
293: PetscCall(STBackTransform((eps)->st,1,&__kr0,&__ki0)); \
294: if ((n)>(k)) { \
295: __krn=(eigr)[k];__kin=(eigi)[k]; \
296: PetscCall(STBackTransform((eps)->st,1,&__krn,&__kin)); \
297: } \
298: if ((eps)->which==EPS_LARGEST_MAGNITUDE || (eps)->which==EPS_SMALLEST_MAGNITUDE) { \
299: ((EPSStoppingCtx)(eps)->stoppingctx)->firstev = SlepcAbsEigenvalue(__kr0,__ki0); \
300: ((EPSStoppingCtx)(eps)->stoppingctx)->lastev = SlepcAbsEigenvalue(__kr,__ki); \
301: ((EPSStoppingCtx)(eps)->stoppingctx)->firstnc = SlepcAbsEigenvalue(__krn,__kin); \
302: } else { \
303: ((EPSStoppingCtx)(eps)->stoppingctx)->firstev = PetscRealPart(__kr0); \
304: ((EPSStoppingCtx)(eps)->stoppingctx)->lastev = PetscRealPart(__kr); \
305: ((EPSStoppingCtx)(eps)->stoppingctx)->firstnc = PetscRealPart(__krn); \
306: } \
307: ((EPSStoppingCtx)(eps)->stoppingctx)->errest = (err_est)[k]; \
308: ((EPSStoppingCtx)(eps)->stoppingctx)->napprox = (n); \
309: } \
310: } while (0)
312: /*
313: EPS_SetInnerProduct - set B matrix for inner product if appropriate.
314: */
315: static inline PetscErrorCode EPS_SetInnerProduct(EPS eps)
316: {
317: Mat B;
319: PetscFunctionBegin;
320: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
321: if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
322: PetscCall(STGetBilinearForm(eps->st,&B));
323: PetscCall(BVSetMatrix(eps->V,B,PetscNot(eps->ispositive)));
324: if (eps->twosided) PetscCall(BVSetMatrix(eps->W,B,PetscNot(eps->ispositive)));
325: PetscCall(MatDestroy(&B));
326: } else PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*
331: EPS_Purify - purify the first k vectors in the V basis
332: */
333: static inline PetscErrorCode EPS_Purify(EPS eps,PetscInt k)
334: {
335: PetscInt i;
336: Vec v,z;
338: PetscFunctionBegin;
339: PetscCall(BVCreateVec(eps->V,&v));
340: for (i=0;i<k;i++) {
341: PetscCall(BVCopyVec(eps->V,i,v));
342: PetscCall(BVGetColumn(eps->V,i,&z));
343: PetscCall(STApply(eps->st,v,z));
344: PetscCall(BVRestoreColumn(eps->V,i,&z));
345: }
346: PetscCall(VecDestroy(&v));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*
351: EPS_KSPSetOperators - Sets the KSP matrices, see also ST_KSPSetOperators()
352: */
353: static inline PetscErrorCode EPS_KSPSetOperators(KSP ksp,Mat A,Mat B)
354: {
355: const char *prefix;
357: PetscFunctionBegin;
358: PetscCall(KSPSetOperators(ksp,A,B));
359: PetscCall(MatGetOptionsPrefix(B,&prefix));
360: if (!prefix) {
361: /* set Mat prefix to be the same as KSP to enable setting command-line options (e.g. MUMPS)
362: only applies if the Mat has no user-defined prefix */
363: PetscCall(KSPGetOptionsPrefix(ksp,&prefix));
364: PetscCall(MatSetOptionsPrefix(B,prefix));
365: }
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /*
370: EPS_GetActualConverged - Gets the actual value of nconv; in special cases the
371: number of available eigenvalues is larger than the computed ones
372: */
373: static inline PetscErrorCode EPS_GetActualConverged(EPS eps,PetscInt *nconv)
374: {
375: PetscFunctionBegin;
376: *nconv = eps->nconv;
377: if (eps->isstructured) {
378: if (eps->problem_type == EPS_BSE && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
379: else if (eps->problem_type == EPS_HAMILT && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
380: else if (eps->problem_type == EPS_LREP && (eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE)) *nconv *= 2;
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: static inline PetscErrorCode EPS_GetEigenvector_BSE(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
386: {
387: PetscInt k;
388: Vec v0,v1,w,w0,w1;
389: Mat H;
390: IS is[2];
392: PetscFunctionBegin;
393: PetscCheck(eps->which == EPS_SMALLEST_MAGNITUDE || eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)(eps)),PETSC_ERR_PLIB,"Inconsistent state");
394: /* BSE problem, even index is +lambda, odd index is -lambda */
395: k = eps->perm[i/2];
396: if (i%2) {
397: /* eigenvector of -lambda is J*conj(x) where J=[0 I; I 0] and x is eigenvector of lambda */
398: PetscCall(VecDuplicate(Vr?Vr:Vi,&w));
399: PetscCall(STGetMatrix(eps->st,0,&H));
400: PetscCall(MatNestGetISs(H,is,NULL));
401: if (Vr) {
402: PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],w,NULL));
403: PetscCall(VecConjugate(w));
404: PetscCall(VecGetSubVector(w,is[0],&w0));
405: PetscCall(VecGetSubVector(w,is[1],&w1));
406: PetscCall(VecGetSubVector(Vr,is[0],&v0));
407: PetscCall(VecGetSubVector(Vr,is[1],&v1));
408: PetscCall(VecCopy(w1,v0));
409: PetscCall(VecCopy(w0,v1));
410: PetscCall(VecRestoreSubVector(w,is[0],&w0));
411: PetscCall(VecRestoreSubVector(w,is[1],&w1));
412: PetscCall(VecRestoreSubVector(Vr,is[0],&v0));
413: PetscCall(VecRestoreSubVector(Vr,is[1],&v1));
414: }
415: #if !PetscDefined(USE_COMPLEX)
416: if (Vi) {
417: PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],NULL,w));
418: PetscCall(VecScale(w,-1.0));
419: PetscCall(VecGetSubVector(w,is[0],&w0));
420: PetscCall(VecGetSubVector(w,is[1],&w1));
421: PetscCall(VecGetSubVector(Vi,is[0],&v0));
422: PetscCall(VecGetSubVector(Vi,is[1],&v1));
423: PetscCall(VecCopy(w1,v0));
424: PetscCall(VecCopy(w0,v1));
425: PetscCall(VecRestoreSubVector(w,is[0],&w0));
426: PetscCall(VecRestoreSubVector(w,is[1],&w1));
427: PetscCall(VecRestoreSubVector(Vi,is[0],&v0));
428: PetscCall(VecRestoreSubVector(Vi,is[1],&v1));
429: }
430: #endif
431: PetscCall(VecDestroy(&w));
432: } else {
433: PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
434: }
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static inline PetscErrorCode EPS_GetEigenvector_HAMILT(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
439: {
440: #if !PetscDefined(USE_COMPLEX)
441: PetscInt k;
442: Vec w;
443: PetscInt k0,k1,k2,iquad;
444: PetscReal nrm,nrmr=0.0,nrmi=0.0,sgn;
445: #endif
447: PetscFunctionBegin;
448: #if !PetscDefined(USE_COMPLEX)
449: k = eps->perm[i/2];
450: if (eps->eigi[k]==0.0) { /* real eigenvalue */
451: if (Vr) {
452: PetscCall(BVCopyVec(V,k+eps->ncv/2+1,Vr));
453: PetscCall(BVGetColumn(V,k,&w));
454: PetscCall(VecAXPY(Vr,(i%2)?-eps->eigr[k]:eps->eigr[k],w));
455: PetscCall(BVRestoreColumn(V,k,&w));
456: PetscCall(VecNorm(Vr,NORM_2,&nrmr));
457: }
458: if (Vi) PetscCall(VecZeroEntries(Vi));
459: nrm = nrmr;
460: } else if (eps->eigr[k]==0.0 ) { /* purely imaginary eigenvalue */
461: if (Vr) {
462: PetscCall(BVCopyVec(V,k+eps->ncv/2+1,Vr));
463: PetscCall(VecNorm(Vr,NORM_2,&nrmr));
464: }
465: if (Vi) {
466: PetscCall(BVCopyVec(V,k,Vi));
467: PetscCall(VecScale(Vi,(i%2)?-eps->eigi[k]:eps->eigi[k]));
468: PetscCall(VecNorm(Vi,NORM_2,&nrmi));
469: }
470: nrm = SlepcAbs(nrmr,nrmi);
471: } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
472: iquad = i%2; /* index within the 4 values */
473: if (i>=2) {
474: k2 = eps->perm[(i-2)/2];
475: if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
476: }
477: k0 = (iquad<2)? k: k2;
478: k1 = k0+1;
479: /* Vr+Vi*i obtained as eig*u+v where u=ur+ui*i is stored in cols k0 (ur) and k1 (ui) and
480: v=vr+vi*i is in cols shifted by ncv/2+1.
481: For lambda=eigr+eigi*i:
482: Vr+Vi*i = (eigr+eigi*i)(ur+ui*i) + vr+vi*i
483: Vr+Vi*i = eigr*(ur+ui*i) - eigi*ui+eigi*ur*i + vr+vi*i
484: Vr+Vi*i = eigr*ur-eigi*ui+vr + (eigi*ur+eigr*ui+vi)*i
485: For -conj(lambda): eigr, ui and vi have the signs changed
486: For -lambda: eigr and eigi have the signs changed
487: For conj(lambda): eigi, ui and vi have the signs changed */
488: if (Vr) {
489: sgn = (iquad<2)? -1.0: 1.0;
490: PetscCall(BVCopyVec(V,k0,Vr)); /* ur */
491: PetscCall(VecScale(Vr,sgn*eps->eigr[k0]));
492: PetscCall(BVGetColumn(V,k1,&w)); /* ui */
493: PetscCall(VecAXPY(Vr,-sgn*eps->eigi[k0],w));
494: PetscCall(BVRestoreColumn(V,k1,&w));
495: PetscCall(BVGetColumn(V,k0+eps->ncv/2+1,&w)); /* vr */
496: PetscCall(VecAXPY(Vr,1.0,w));
497: PetscCall(BVRestoreColumn(V,k0+eps->ncv/2+1,&w));
498: PetscCall(VecNorm(Vr,NORM_2,&nrmr));
499: }
500: if (Vi) {
501: sgn = (iquad%2)? -1.0: 1.0;
502: PetscCall(BVCopyVec(V,k0,Vi)); /* ur */
503: PetscCall(VecScale(Vi,sgn*eps->eigi[k0]));
504: PetscCall(BVGetColumn(V,k1,&w)); /* ui */
505: PetscCall(VecAXPY(Vi,sgn*eps->eigr[k0],w));
506: PetscCall(BVRestoreColumn(V,k1,&w));
507: PetscCall(BVGetColumn(V,k1+eps->ncv/2+1,&w)); /* vi */
508: sgn = (iquad%3)? 1.0: -1.0;
509: PetscCall(VecAXPY(Vi,sgn,w));
510: PetscCall(BVRestoreColumn(V,k1+eps->ncv/2+1,&w));
511: PetscCall(VecNorm(Vi,NORM_2,&nrmi));
512: }
513: nrm = SlepcAbs(nrmr,nrmi);
514: }
515: if (Vr) PetscCall(VecScale(Vr,1.0/nrm));
516: if (Vi) PetscCall(VecScale(Vi,1.0/nrm));
517: #endif
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: static inline PetscErrorCode EPS_GetEigenvector_LREP(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
522: {
523: PetscInt k;
524: Vec v;
525: Mat H;
526: IS is[2];
528: PetscFunctionBegin;
529: /* LREP problem, even index is +lambda, odd index is -lambda */
530: k = eps->perm[i/2];
531: PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
532: if (i%2) {
533: /* eigenvector of -lambda is S*x where S=[I 0; 0 -I] and x is eigenvector of lambda */
534: PetscCall(STGetMatrix(eps->st,0,&H));
535: PetscCall(MatNestGetISs(H,is,NULL));
536: PetscCall(VecGetSubVector(Vr,is[1],&v));
537: PetscCall(VecScale(v,-1.0));
538: PetscCall(VecRestoreSubVector(Vr,is[1],&v));
539: }
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: /*
544: EPS_GetEigenvector - Gets the i-th eigenvector taking into account the case
545: where i exceeds the number of computed vectors (structure-preserving solver).
546: The argument V should be eps->V for right eigenvectors, eps->W for left ones.
547: */
548: static inline PetscErrorCode EPS_GetEigenvector(EPS eps,BV V,PetscInt i,Vec Vr,Vec Vi)
549: {
550: PetscInt k;
551: PetscBool reduced;
552: Mat H;
554: PetscFunctionBegin;
555: if (!eps->isstructured) {
556: k = eps->perm[i];
557: PetscCall(BV_GetEigenvector(V,k,eps->eigi[k],Vr,Vi));
558: } else {
559: switch (eps->problem_type) {
560: case EPS_BSE:
561: PetscCall(EPS_GetEigenvector_BSE(eps,V,i,Vr,Vi));
562: break;
563: case EPS_HAMILT:
564: PetscCall(EPS_GetEigenvector_HAMILT(eps,V,i,Vr,Vi));
565: break;
566: case EPS_LREP:
567: PetscCall(STGetMatrix(eps->st,0,&H));
568: PetscCall(SlepcCheckMatLREPReduced(H,&reduced));
569: if (reduced) PetscCall(EPS_GetEigenvector_LREP(eps,V,i,Vr,Vi));
570: else PetscCall(EPS_GetEigenvector_BSE(eps,V,i,Vr,Vi));
571: break;
572: default:
573: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Inconsistent state");
574: }
575: }
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: SLEPC_INTERN PetscErrorCode EPSSetWhichEigenpairs_Default(EPS);
580: SLEPC_INTERN PetscErrorCode EPSSetDimensions_Default(EPS,PetscInt*,PetscInt*,PetscInt*);
581: SLEPC_INTERN PetscErrorCode EPSBackTransform_Default(EPS);
582: SLEPC_INTERN PetscErrorCode EPSComputeVectors(EPS);
583: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Hermitian(EPS);
584: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Schur(EPS);
585: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Indefinite(EPS);
586: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Twosided(EPS);
587: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Slice(EPS);
588: SLEPC_INTERN PetscErrorCode EPSComputeResidualNorm_Private(EPS,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
589: SLEPC_INTERN PetscErrorCode EPSComputeRitzVector(EPS,PetscScalar*,PetscScalar*,BV,Vec,Vec);
590: SLEPC_INTERN PetscErrorCode EPSGetStartVector(EPS,PetscInt,PetscBool*);
591: SLEPC_INTERN PetscErrorCode EPSGetLeftStartVector(EPS,PetscInt,PetscBool*);
592: SLEPC_INTERN PetscErrorCode MatEstimateSpectralRange_EPS(Mat,PetscReal*,PetscReal*);
594: /* Private functions of the solver implementations */
596: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
597: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi1(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
598: SLEPC_INTERN PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscInt,PetscInt,PetscReal,PetscReal,PetscReal,PetscInt*);
599: SLEPC_INTERN PetscErrorCode EPSPseudoLanczos(EPS,PetscReal*,PetscReal*,PetscReal*,PetscInt,PetscInt*,PetscBool*,PetscBool*,PetscReal*,Vec);
600: SLEPC_INTERN PetscErrorCode EPSBuildBalance_Krylov(EPS);
601: SLEPC_INTERN PetscErrorCode EPSSetDefaultST(EPS);
602: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_Precond(EPS);
603: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_GMRES(EPS);
604: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_NoFactor(EPS);
605: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Basic(EPS);
606: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Default(EPS);