Actual source code: epsimpl.h

slepc-3.16.1 2021-11-17
Report Typos and Errors
  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: #if !defined(SLEPCEPSIMPL_H)
 12: #define SLEPCEPSIMPL_H

 14: #include <slepceps.h>
 15: #include <slepc/private/slepcimpl.h>

 17: SLEPC_EXTERN PetscBool EPSRegisterAllCalled;
 18: SLEPC_EXTERN PetscBool EPSMonitorRegisterAllCalled;
 19: SLEPC_EXTERN PetscErrorCode EPSRegisterAll(void);
 20: SLEPC_EXTERN PetscErrorCode EPSMonitorRegisterAll(void);
 21: SLEPC_EXTERN PetscLogEvent EPS_SetUp,EPS_Solve,EPS_CISS_SVD;

 23: typedef struct _EPSOps *EPSOps;

 25: struct _EPSOps {
 26:   PetscErrorCode (*solve)(EPS);
 27:   PetscErrorCode (*setup)(EPS);
 28:   PetscErrorCode (*setupsort)(EPS);
 29:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,EPS);
 30:   PetscErrorCode (*publishoptions)(EPS);
 31:   PetscErrorCode (*destroy)(EPS);
 32:   PetscErrorCode (*reset)(EPS);
 33:   PetscErrorCode (*view)(EPS,PetscViewer);
 34:   PetscErrorCode (*backtransform)(EPS);
 35:   PetscErrorCode (*computevectors)(EPS);
 36:   PetscErrorCode (*setdefaultst)(EPS);
 37: };

 39: /*
 40:    Maximum number of monitors you can run with a single EPS
 41: */
 42: #define MAXEPSMONITORS 5

 44: /*
 45:    The solution process goes through several states
 46: */
 47: typedef enum { EPS_STATE_INITIAL,
 48:                EPS_STATE_SETUP,
 49:                EPS_STATE_SOLVED,
 50:                EPS_STATE_EIGENVECTORS } EPSStateType;

 52: /*
 53:    To classify the different solvers into categories
 54: */
 55: typedef enum { EPS_CATEGORY_KRYLOV,      /* Krylov solver: relies on STApply and STBackTransform (same as OTHER) */
 56:                EPS_CATEGORY_PRECOND,     /* Preconditioned solver: uses ST only to manage preconditioner */
 57:                EPS_CATEGORY_CONTOUR,     /* Contour integral: ST used to solve linear systems at integration points */
 58:                EPS_CATEGORY_OTHER } EPSSolverType;

 60: /*
 61:    To check for unsupported features at EPSSetUp_XXX()
 62: */
 63: typedef enum { EPS_FEATURE_BALANCE=1,       /* balancing */
 64:                EPS_FEATURE_ARBITRARY=2,     /* arbitrary selection of eigepairs */
 65:                EPS_FEATURE_REGION=4,        /* nontrivial region for filtering */
 66:                EPS_FEATURE_EXTRACTION=8,    /* extraction technique different from Ritz */
 67:                EPS_FEATURE_CONVERGENCE=16,  /* convergence test selected by user */
 68:                EPS_FEATURE_STOPPING=32,     /* stopping test */
 69:                EPS_FEATURE_TWOSIDED=64      /* two-sided variant */
 70:              } EPSFeatureType;

 72: /*
 73:    Defines the EPS data structure
 74: */
 75: struct _p_EPS {
 76:   PETSCHEADER(struct _EPSOps);
 77:   /*------------------------- User parameters ---------------------------*/
 78:   PetscInt       max_it;           /* maximum number of iterations */
 79:   PetscInt       nev;              /* number of eigenvalues to compute */
 80:   PetscInt       ncv;              /* number of basis vectors */
 81:   PetscInt       mpd;              /* maximum dimension of projected problem */
 82:   PetscInt       nini,ninil;       /* number of initial vectors (negative means not copied yet) */
 83:   PetscInt       nds;              /* number of basis vectors of deflation space */
 84:   PetscScalar    target;           /* target value */
 85:   PetscReal      tol;              /* tolerance */
 86:   EPSConv        conv;             /* convergence test */
 87:   EPSStop        stop;             /* stopping test */
 88:   EPSWhich       which;            /* which part of the spectrum to be sought */
 89:   PetscReal      inta,intb;        /* interval [a,b] for spectrum slicing */
 90:   EPSProblemType problem_type;     /* which kind of problem to be solved */
 91:   EPSExtraction  extraction;       /* which kind of extraction to be applied */
 92:   EPSBalance     balance;          /* the balancing method */
 93:   PetscInt       balance_its;      /* number of iterations of the balancing method */
 94:   PetscReal      balance_cutoff;   /* cutoff value for balancing */
 95:   PetscBool      trueres;          /* whether the true residual norm must be computed */
 96:   PetscBool      trackall;         /* whether all the residuals must be computed */
 97:   PetscBool      purify;           /* whether eigenvectors need to be purified */
 98:   PetscBool      twosided;         /* whether to compute left eigenvectors (two-sided solver) */

100:   /*-------------- User-provided functions and contexts -----------------*/
101:   PetscErrorCode (*converged)(EPS,PetscScalar,PetscScalar,PetscReal,PetscReal*,void*);
102:   PetscErrorCode (*convergeduser)(EPS,PetscScalar,PetscScalar,PetscReal,PetscReal*,void*);
103:   PetscErrorCode (*convergeddestroy)(void*);
104:   PetscErrorCode (*stopping)(EPS,PetscInt,PetscInt,PetscInt,PetscInt,EPSConvergedReason*,void*);
105:   PetscErrorCode (*stoppinguser)(EPS,PetscInt,PetscInt,PetscInt,PetscInt,EPSConvergedReason*,void*);
106:   PetscErrorCode (*stoppingdestroy)(void*);
107:   PetscErrorCode (*arbitrary)(PetscScalar,PetscScalar,Vec,Vec,PetscScalar*,PetscScalar*,void*);
108:   void           *convergedctx;
109:   void           *stoppingctx;
110:   void           *arbitraryctx;
111:   PetscErrorCode (*monitor[MAXEPSMONITORS])(EPS,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,void*);
112:   PetscErrorCode (*monitordestroy[MAXEPSMONITORS])(void**);
113:   void           *monitorcontext[MAXEPSMONITORS];
114:   PetscInt       numbermonitors;

116:   /*----------------- Child objects and working data -------------------*/
117:   ST             st;               /* spectral transformation object */
118:   DS             ds;               /* direct solver object */
119:   BV             V;                /* set of basis vectors and computed eigenvectors */
120:   BV             W;                /* left basis vectors (if left eigenvectors requested) */
121:   RG             rg;               /* optional region for filtering */
122:   SlepcSC        sc;               /* sorting criterion data */
123:   Vec            D;                /* diagonal matrix for balancing */
124:   Vec            *IS,*ISL;         /* references to user-provided initial spaces */
125:   Vec            *defl;            /* references to user-provided deflation space */
126:   PetscScalar    *eigr,*eigi;      /* real and imaginary parts of eigenvalues */
127:   PetscReal      *errest;          /* error estimates */
128:   PetscScalar    *rr,*ri;          /* values computed by user's arbitrary selection function */
129:   PetscInt       *perm;            /* permutation for eigenvalue ordering */
130:   PetscInt       nwork;            /* number of work vectors */
131:   Vec            *work;            /* work vectors */
132:   void           *data;            /* placeholder for solver-specific stuff */

134:   /* ----------------------- Status variables --------------------------*/
135:   EPSStateType   state;            /* initial -> setup -> solved -> eigenvectors */
136:   EPSSolverType  categ;            /* solver category */
137:   PetscInt       nconv;            /* number of converged eigenvalues */
138:   PetscInt       its;              /* number of iterations so far computed */
139:   PetscInt       n,nloc;           /* problem dimensions (global, local) */
140:   PetscReal      nrma,nrmb;        /* computed matrix norms */
141:   PetscBool      useds;            /* whether the solver uses the DS object or not */
142:   PetscBool      isgeneralized;
143:   PetscBool      ispositive;
144:   PetscBool      ishermitian;
145:   EPSConvergedReason reason;
146: };

148: /*
149:     Macros to test valid EPS arguments
150: */
151: #if !defined(PETSC_USE_DEBUG)

153: #define EPSCheckSolved(h,arg) do {(void)(h);} while (0)

155: #else

157: #define EPSCheckSolved(h,arg) \
158:   do { \
159:     if ((h)->state<EPS_STATE_SOLVED) SETERRQ1(PetscObjectComm((PetscObject)(h)),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSolve() first: Parameter #%d",arg); \
160:   } while (0)

162: #endif

164: /*
165:     Macros to check settings at EPSSetUp()
166: */

168: /* EPSCheckHermitianDefinite: the problem is HEP or GHEP */
169: #define EPSCheckHermitianDefiniteCondition(eps,condition,msg) \
170:   do { \
171:     if (condition) { \
172:       if (!(eps)->ishermitian) SETERRQ3(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); \
173:       else if ((eps)->isgeneralized && !(eps)->ispositive) SETERRQ3(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); \
174:     } \
175:   } while (0)
176: #define EPSCheckHermitianDefinite(eps) EPSCheckHermitianDefiniteCondition(eps,PETSC_TRUE,"")

178: /* EPSCheckHermitian: the problem is HEP, GHEP, or GHIEP */
179: #define EPSCheckHermitianCondition(eps,condition,msg) \
180:   do { \
181:     if (condition) { \
182:       if (!(eps)->ishermitian) SETERRQ3(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); \
183:     } \
184:   } while (0)
185: #define EPSCheckHermitian(eps) EPSCheckHermitianCondition(eps,PETSC_TRUE,"")

187: /* EPSCheckDefinite: the problem is not GHIEP */
188: #define EPSCheckDefiniteCondition(eps,condition,msg) \
189:   do { \
190:     if (condition) { \
191:       if ((eps)->isgeneralized && (eps)->ishermitian && !(eps)->ispositive) SETERRQ3(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); \
192:     } \
193:   } while (0)
194: #define EPSCheckDefinite(eps) EPSCheckDefiniteCondition(eps,PETSC_TRUE,"")

196: /* EPSCheckStandard: the problem is HEP or NHEP */
197: #define EPSCheckStandardCondition(eps,condition,msg) \
198:   do { \
199:     if (condition) { \
200:       if ((eps)->isgeneralized) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot be used for generalized problems",((PetscObject)(eps))->type_name,(msg)); \
201:     } \
202:   } while (0)
203: #define EPSCheckStandard(eps) EPSCheckStandardCondition(eps,PETSC_TRUE,"")

205: /* EPSCheckSinvert: shift-and-invert ST */
206: #define EPSCheckSinvertCondition(eps,condition,msg) \
207:   do { \
208:     if (condition) { \
209:       PetscBool __flg; \
210:       PetscObjectTypeCompare((PetscObject)(eps)->st,STSINVERT,&__flg); \
211:       if (!__flg) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires a shift-and-invert spectral transform",((PetscObject)(eps))->type_name,(msg)); \
212:     } \
213:   } while (0)
214: #define EPSCheckSinvert(eps) EPSCheckSinvertCondition(eps,PETSC_TRUE,"")

216: /* EPSCheckSinvertCayley: shift-and-invert or Cayley ST */
217: #define EPSCheckSinvertCayleyCondition(eps,condition,msg) \
218:   do { \
219:     if (condition) { \
220:       PetscBool __flg; \
221:       PetscObjectTypeCompareAny((PetscObject)(eps)->st,&__flg,STSINVERT,STCAYLEY,""); \
222:       if (!__flg) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s requires shift-and-invert or Cayley transform",((PetscObject)(eps))->type_name,(msg)); \
223:     } \
224:   } while (0)
225: #define EPSCheckSinvertCayley(eps) EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE,"")

227: /* Check for unsupported features */
228: #define EPSCheckUnsupportedCondition(eps,mask,condition,msg) \
229:   do { \
230:     if (condition) { \
231:       if (((mask) & EPS_FEATURE_BALANCE) && (eps)->balance!=EPS_BALANCE_NONE) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support balancing",((PetscObject)(eps))->type_name,(msg)); \
232:       if (((mask) & EPS_FEATURE_ARBITRARY) && (eps)->arbitrary) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support arbitrary selection of eigenpairs",((PetscObject)(eps))->type_name,(msg)); \
233:       if ((mask) & EPS_FEATURE_REGION) { \
234:         PetscBool      __istrivial; \
235:         PetscErrorCode __RGIsTrivial((eps)->rg,&__istrivial);CHKERRQ(__ierr); \
236:         if (!__istrivial) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s does not support region filtering",((PetscObject)(eps))->type_name,(msg)); \
237:       } \
238:       if (((mask) & EPS_FEATURE_EXTRACTION) && (eps)->extraction!=EPS_RITZ) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports Ritz extraction",((PetscObject)(eps))->type_name,(msg)); \
239:       if (((mask) & EPS_FEATURE_CONVERGENCE) && (eps)->converged!=EPSConvergedRelative) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default convergence test",((PetscObject)(eps))->type_name,(msg)); \
240:       if (((mask) & EPS_FEATURE_STOPPING) && (eps)->stopping!=EPSStoppingBasic) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s only supports the default stopping test",((PetscObject)(eps))->type_name,(msg)); \
241:       if (((mask) & EPS_FEATURE_TWOSIDED) && (eps)->twosided) SETERRQ2(PetscObjectComm((PetscObject)(eps)),PETSC_ERR_SUP,"The solver '%s'%s cannot compute left eigenvectors (no two-sided variant)",((PetscObject)(eps))->type_name,(msg)); \
242:     } \
243:   } while (0)
244: #define EPSCheckUnsupported(eps,mask) EPSCheckUnsupportedCondition(eps,mask,PETSC_TRUE,"")

246: /* Check for ignored features */
247: #define EPSCheckIgnoredCondition(eps,mask,condition,msg) \
248:   do { \
249:     PetscErrorCode __ierr; \
250:     if (condition) { \
251:       if (((mask) & EPS_FEATURE_BALANCE) && (eps)->balance!=EPS_BALANCE_NONE) { __PetscInfo2((eps),"The solver '%s'%s ignores the balancing settings\n",((PetscObject)(eps))->type_name,(msg)); } \
252:       if (((mask) & EPS_FEATURE_ARBITRARY) && (eps)->arbitrary) { __PetscInfo2((eps),"The solver '%s'%s ignores the settings for arbitrary selection of eigenpairs\n",((PetscObject)(eps))->type_name,(msg)); } \
253:       if ((mask) & EPS_FEATURE_REGION) { \
254:         PetscBool __istrivial; \
255:         __RGIsTrivial((eps)->rg,&__istrivial);CHKERRQ(__ierr); \
256:         if (!__istrivial) { __PetscInfo2((eps),"The solver '%s'%s ignores the specified region\n",((PetscObject)(eps))->type_name,(msg)); } \
257:       } \
258:       if (((mask) & EPS_FEATURE_EXTRACTION) && (eps)->extraction!=EPS_RITZ) { __PetscInfo2((eps),"The solver '%s'%s ignores the extraction settings\n",((PetscObject)(eps))->type_name,(msg)); } \
259:       if (((mask) & EPS_FEATURE_CONVERGENCE) && (eps)->converged!=EPSConvergedRelative) { __PetscInfo2((eps),"The solver '%s'%s ignores the convergence test settings\n",((PetscObject)(eps))->type_name,(msg)); } \
260:       if (((mask) & EPS_FEATURE_STOPPING) && (eps)->stopping!=EPSStoppingBasic) { __PetscInfo2((eps),"The solver '%s'%s ignores the stopping test settings\n",((PetscObject)(eps))->type_name,(msg)); } \
261:       if (((mask) & EPS_FEATURE_TWOSIDED) && (eps)->twosided) { __PetscInfo2((eps),"The solver '%s'%s ignores the two-sided flag\n",((PetscObject)(eps))->type_name,(msg)); } \
262:     } \
263:   } while (0)
264: #define EPSCheckIgnored(eps,mask) EPSCheckIgnoredCondition(eps,mask,PETSC_TRUE,"")

266: /*
267:   EPS_SetInnerProduct - set B matrix for inner product if appropriate.
268: */
269: PETSC_STATIC_INLINE PetscErrorCode EPS_SetInnerProduct(EPS eps)
270: {
272:   Mat            B;

275:   if (!eps->V) { EPSGetBV(eps,&eps->V); }
276:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
277:     STGetBilinearForm(eps->st,&B);
278:     BVSetMatrix(eps->V,B,PetscNot(eps->ispositive));
279:     MatDestroy(&B);
280:   } else {
281:     BVSetMatrix(eps->V,NULL,PETSC_FALSE);
282:   }
283:   return(0);
284: }

286: /*
287:   EPS_Purify - purify the first k vectors in the V basis
288: */
289: PETSC_STATIC_INLINE PetscErrorCode EPS_Purify(EPS eps,PetscInt k)
290: {
292:   PetscInt       i;
293:   Vec            v,z;

296:   BVCreateVec(eps->V,&v);
297:   for (i=0;i<k;i++) {
298:     BVCopyVec(eps->V,i,v);
299:     BVGetColumn(eps->V,i,&z);
300:     STApply(eps->st,v,z);
301:     BVRestoreColumn(eps->V,i,&z);
302:   }
303:   VecDestroy(&v);
304:   return(0);
305: }

307: SLEPC_INTERN PetscErrorCode EPSSetWhichEigenpairs_Default(EPS);
308: SLEPC_INTERN PetscErrorCode EPSSetDimensions_Default(EPS,PetscInt,PetscInt*,PetscInt*);
309: SLEPC_INTERN PetscErrorCode EPSBackTransform_Default(EPS);
310: SLEPC_INTERN PetscErrorCode EPSComputeVectors(EPS);
311: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Hermitian(EPS);
312: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Schur(EPS);
313: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Indefinite(EPS);
314: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Twosided(EPS);
315: SLEPC_INTERN PetscErrorCode EPSComputeVectors_Slice(EPS);
316: SLEPC_INTERN PetscErrorCode EPSComputeResidualNorm_Private(EPS,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
317: SLEPC_INTERN PetscErrorCode EPSComputeRitzVector(EPS,PetscScalar*,PetscScalar*,BV,Vec,Vec);
318: SLEPC_INTERN PetscErrorCode EPSGetStartVector(EPS,PetscInt,PetscBool*);
319: SLEPC_INTERN PetscErrorCode EPSGetLeftStartVector(EPS,PetscInt,PetscBool*);

321: /* Private functions of the solver implementations */

323: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
324: SLEPC_INTERN PetscErrorCode EPSDelayedArnoldi1(EPS,PetscScalar*,PetscInt,PetscInt,PetscInt*,PetscReal*,PetscBool*);
325: SLEPC_INTERN PetscErrorCode EPSKrylovConvergence(EPS,PetscBool,PetscInt,PetscInt,PetscReal,PetscReal,PetscReal,PetscInt*);
326: SLEPC_INTERN PetscErrorCode EPSPseudoLanczos(EPS,PetscReal*,PetscReal*,PetscReal*,PetscInt,PetscInt*,PetscBool*,PetscBool*,PetscReal*,Vec);
327: SLEPC_INTERN PetscErrorCode EPSBuildBalance_Krylov(EPS);
328: SLEPC_INTERN PetscErrorCode EPSSetDefaultST(EPS);
329: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_Precond(EPS);
330: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_GMRES(EPS);
331: SLEPC_INTERN PetscErrorCode EPSSetDefaultST_NoFactor(EPS);
332: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Basic(EPS);
333: SLEPC_INTERN PetscErrorCode EPSSetUpSort_Default(EPS);

335: #endif