Actual source code: davidson.h

slepc-3.18.0 2022-10-01
Report Typos and Errors
  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: */
 10: /*
 11:    Method: General Davidson Method (includes GD and JD)

 13:    References:
 14:     - Ernest R. Davidson. Super-matrix methods. Computer Physics Communications,
 15:       53:49-60, May 1989.
 16: */

 18: #include <slepc/private/epsimpl.h>
 19: #include <slepc/private/vecimplslepc.h>

 21: typedef PetscInt MatType_t;
 22: #define DVD_MAT_HERMITIAN (1<<1)
 23: #define DVD_MAT_NEG_DEF (1<<2)
 24: #define DVD_MAT_POS_DEF (1<<3)
 25: #define DVD_MAT_SINGULAR (1<<4)
 26: #define DVD_MAT_COMPLEX (1<<5)
 27: #define DVD_MAT_IMPLICIT (1<<6)
 28: #define DVD_MAT_IDENTITY (1<<7)
 29: #define DVD_MAT_DIAG (1<<8)
 30: #define DVD_MAT_TRIANG (1<<9)
 31: #define DVD_MAT_UTRIANG (1<<9)
 32: #define DVD_MAT_LTRIANG (1<<10)
 33: #define DVD_MAT_UNITARY (1<<11)

 35: typedef PetscInt EPType_t;
 36: #define DVD_EP_STD (1<<1)
 37: #define DVD_EP_HERMITIAN (1<<2)
 38: #define DVD_EP_INDEFINITE (1<<3)

 40: #define DVD_IS(T,P) ((T) & (P))
 41: #define DVD_ISNOT(T,P) (((T) & (P)) ^ (P))

 43: struct _dvdDashboard;
 44: typedef PetscErrorCode (*dvdCallback)(struct _dvdDashboard*);
 45: typedef struct _dvdFunctionList {
 46:   dvdCallback f;
 47:   struct _dvdFunctionList *next;
 48: } dvdFunctionList;

 50: typedef enum {
 51:   DVD_HARM_NONE,
 52:   DVD_HARM_RR,
 53:   DVD_HARM_RRR,
 54:   DVD_HARM_REIGS,
 55:   DVD_HARM_LEIGS
 56: } HarmType_t;

 58: typedef enum {
 59:   DVD_INITV_CLASSIC,
 60:   DVD_INITV_KRYLOV
 61: } InitType_t;

 63: /*
 64:    Dashboard struct: contains the methods that will be employed and the tuning
 65:    options.
 66: */

 68: typedef struct _dvdDashboard {
 69:   /**** Function steps ****/
 70:   /* Initialize V */
 71:   PetscErrorCode (*initV)(struct _dvdDashboard*);
 72:   void *initV_data;

 74:   /* Find the approximate eigenpairs from V */
 75:   PetscErrorCode (*calcPairs)(struct _dvdDashboard*);
 76:   void *calcPairs_data;

 78:   /* Eigenpair test for convergence */
 79:   PetscBool (*testConv)(struct _dvdDashboard*,PetscScalar,PetscScalar,PetscReal,PetscReal*);
 80:   void *testConv_data;

 82:   /* Improve the selected eigenpairs */
 83:   PetscErrorCode (*improveX)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt*);
 84:   void *improveX_data;

 86:   /* Check for restarting */
 87:   PetscErrorCode (*isRestarting)(struct _dvdDashboard*,PetscBool*);
 88:   void *isRestarting_data;

 90:   /* Perform restarting */
 91:   PetscErrorCode (*restartV)(struct _dvdDashboard*);
 92:   void *restartV_data;

 94:   /* Update V */
 95:   PetscErrorCode (*updateV)(struct _dvdDashboard*);
 96:   void *updateV_data;

 98:   /**** Problem specification ****/
 99:   Mat         A,B;            /* problem matrices */
100:   MatType_t   sA,sB;          /* matrix specifications */
101:   EPType_t    sEP;            /* problem specifications */
102:   PetscInt    nev;            /* number of eigenpairs */
103:   EPSWhich    which;          /* spectrum selection */
104:   PetscBool   withTarget;     /* if there is a target */
105:   PetscScalar target[2];      /* target value */
106:   PetscReal   tol;            /* tolerance */
107:   PetscBool   correctXnorm;   /* if true, norm of X are computed */

109:   /**** Subspaces specification ****/
110:   PetscInt nconv;             /* number of converged eigenpairs */
111:   PetscInt npreconv;          /* number of pairs ready to converge */

113:   BV       W;                 /* left basis for harmonic case */
114:   BV       AX;                /* A*V */
115:   BV       BX;                /* B*V */
116:   PetscInt size_D;            /* active vectors */
117:   PetscInt max_size_proj;     /* max size projected problem */
118:   PetscInt max_size_P;        /* max unconverged vectors in the projector */
119:   PetscInt bs;                /* max vectors that expands the subspace every iteration */
120:   EPS      eps;               /* connection to SLEPc */

122:   /**** Auxiliary space ****/
123:   VecPool auxV;               /* auxiliary vectors */
124:   BV      auxBV;              /* auxiliary vectors */

126:   /**** Eigenvalues and errors ****/
127:   PetscScalar *ceigr,*ceigi;  /* converged eigenvalues */
128:   PetscScalar *eigr,*eigi;    /* current eigenvalues */
129:   PetscReal   *nR;            /* residual norm */
130:   PetscReal   *real_nR;       /* original nR */
131:   PetscReal   *nX;            /* X norm */
132:   PetscReal   *real_nX;       /* original nX */
133:   PetscReal   *errest;        /* relative error eigenpairs */
134:   PetscReal   *nBds;          /* B-norms of projected problem  */

136:   /**** Shared function and variables ****/
137:   PetscErrorCode (*e_Vchanged)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt,PetscInt);
138:   void *e_Vchanged_data;
139:   PetscErrorCode (*calcpairs_residual)(struct _dvdDashboard*,PetscInt,PetscInt);
140:   PetscErrorCode (*calcpairs_selectPairs)(struct _dvdDashboard*,PetscInt);
141:   void *calcpairs_residual_data;
142:   PetscErrorCode (*improvex_precond)(struct _dvdDashboard*,PetscInt,Vec,Vec);
143:   void *improvex_precond_data;
144:   PetscErrorCode (*improvex_jd_proj_uv)(struct _dvdDashboard*,PetscInt,PetscInt,Vec*,Vec*,Vec*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscInt);
145:   PetscErrorCode (*improvex_jd_lit)(struct _dvdDashboard*,PetscInt,PetscScalar*,PetscScalar*,PetscInt*,PetscReal*);
146:   PetscErrorCode (*calcpairs_W)(struct _dvdDashboard*);
147:   void *calcpairs_W_data;
148:   PetscErrorCode (*calcpairs_proj_trans)(struct _dvdDashboard*);
149:   PetscErrorCode (*calcpairs_eigs_trans)(struct _dvdDashboard*);
150:   PetscErrorCode (*calcpairs_eig_backtrans)(struct _dvdDashboard*,PetscScalar,PetscScalar,PetscScalar*,PetscScalar*);
151:   PetscErrorCode (*calcpairs_proj_res)(struct _dvdDashboard*,PetscInt,PetscInt,Vec*);
152:   PetscErrorCode (*preTestConv)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt,PetscInt*);
153:   PetscErrorCode (*e_newIteration)(struct _dvdDashboard*);
154:   void *e_newIteration_data;

156:   dvdFunctionList *startList;  /* starting list */
157:   dvdFunctionList *endList;    /* ending list */
158:   dvdFunctionList *destroyList;/* destructor list */

160:   Mat       H,G;               /* projected problem matrices */
161:   Mat       auxM;              /* auxiliary dense matrix */
162:   PetscInt  size_MT;           /* rows in MT */

164:   PetscInt  V_tra_s;
165:   PetscInt  V_tra_e;       /* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
166:   PetscInt  V_new_s;
167:   PetscInt  V_new_e;           /* added to V the columns V_new_s:V_new_e */
168:   PetscBool W_shift;           /* if true W is shifted when vectors converge */
169: } dvdDashboard;

171: typedef struct {
172:   /*------------------------- User parameters ---------------------------*/
173:   PetscInt  blocksize;     /* block size */
174:   PetscInt  initialsize;   /* initial size of V */
175:   PetscInt  minv;          /* size of V after restarting */
176:   PetscInt  plusk;         /* keep plusk eigenvectors from the last iteration */
177:   PetscBool ipB;           /* true if B-ortho is used */
178:   PetscReal fix;           /* the fix parameter */
179:   PetscBool krylovstart;   /* true if the starting subspace is a Krylov basis */
180:   PetscBool dynamic;       /* true if dynamic stopping criterion is used */
181:   PetscBool doubleexp;     /* double expansion in GD (GD2) */

183:   /*----------------- Child objects and working data -------------------*/
184:   dvdDashboard ddb;
185: } EPS_DAVIDSON;

187: static inline PetscErrorCode EPSDavidsonFLAdd(dvdFunctionList **fl,dvdCallback f)
188: {
189:   dvdFunctionList *l;

191:   PetscNew(&l);
192:   l->f = f;
193:   l->next = *fl;
194:   *fl = l;
195:   return 0;
196: }

198: static inline PetscErrorCode EPSDavidsonFLCall(dvdFunctionList *fl,dvdDashboard *d)
199: {
200:   dvdFunctionList *l;

202:   for (l=fl;l;l=l->next) (l->f)(d);
203:   return 0;
204: }

206: static inline PetscErrorCode EPSDavidsonFLDestroy(dvdFunctionList **fl)
207: {
208:   dvdFunctionList *l,*l0;

210:   for (l=*fl;l;l=l0) {
211:     l0 = l->next;
212:     PetscFree(l);
213:   }
214:   *fl = NULL;
215:   return 0;
216: }

218: /*
219:   The blackboard configuration structure: saves information about the memory
220:   and other requirements.

222:   The starting memory structure:

224:   V           W?          AV          BV?          tKZ
225:   |-----------|-----------|-----------|------------|-------|
226:   nev+mpd     nev+mpd     scP+mpd     nev?+mpd     sP+scP
227:               scP+mpd                 scP+mpd

229:   The final memory structure considering W_shift:

231:   cX  V       cY?  W?     cAV AV      BcX? BV?     KZ  tKZ
232:   |---|-------|----|------|---|-------|----|-------|---|---|
233:   nev mpd     nev  mpd    scP mpd     nev  mpd     scP sP    <- shift
234:               scP                     scP                    <- !shift
235: */
236: typedef struct {
237:   PetscInt max_size_V;         /* max size of the searching subspace (mpd) */
238:   PetscInt max_size_X;         /* max size of X (bs) */
239:   PetscInt size_V;             /* real size of V (nev+size_P+mpd) */
240:   PetscInt max_size_oldX;      /* max size of oldX */
241:   PetscInt max_nev;            /* max number of converged pairs */
242:   PetscInt max_size_P;         /* number of computed vectors for the projector */
243:   PetscInt max_size_cP;        /* number of converged vectors in the projectors */
244:   PetscInt max_size_proj;      /* max size projected problem */
245:   PetscInt max_size_cX_proj;   /* max converged vectors in the projected problem */
246:   PetscInt state;              /* method states:
247:                                    0: preconfiguring
248:                                    1: configuring
249:                                    2: running */
250: } dvdBlackboard;

252: #define DVD_STATE_PRECONF 0
253: #define DVD_STATE_CONF 1
254: #define DVD_STATE_RUN 2

256: /* Prototypes of non-static auxiliary functions */
257: SLEPC_INTERN PetscErrorCode dvd_calcpairs_qz(dvdDashboard*,dvdBlackboard*,PetscBool,PetscBool);
258: SLEPC_INTERN PetscErrorCode dvd_improvex_gd2(dvdDashboard*,dvdBlackboard*,KSP,PetscInt);
259: SLEPC_INTERN PetscErrorCode dvd_improvex_jd(dvdDashboard*,dvdBlackboard*,KSP,PetscInt,PetscBool);
260: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard*,dvdBlackboard*);
261: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard*,dvdBlackboard*,PetscInt,PetscReal,PetscReal);
262: SLEPC_INTERN PetscErrorCode dvd_improvex_compute_X(dvdDashboard*,PetscInt,PetscInt,Vec*,PetscScalar*,PetscInt);
263: SLEPC_INTERN PetscErrorCode dvd_initV(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscBool);
264: SLEPC_INTERN PetscErrorCode dvd_orthV(BV,PetscInt,PetscInt);
265: SLEPC_INTERN PetscErrorCode dvd_schm_basic_preconf(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,HarmType_t,KSP,InitType_t,PetscBool,PetscBool,PetscBool);
266: SLEPC_INTERN PetscErrorCode dvd_schm_basic_conf(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,HarmType_t,PetscBool,PetscScalar,KSP,PetscReal,InitType_t,PetscBool,PetscBool,PetscBool,PetscBool);
267: SLEPC_INTERN PetscErrorCode dvd_testconv_slepc(dvdDashboard*,dvdBlackboard*);
268: SLEPC_INTERN PetscErrorCode dvd_managementV_basic(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool,PetscBool);
269: SLEPC_INTERN PetscErrorCode dvd_static_precond_PC(dvdDashboard*,dvdBlackboard*,PC);
270: SLEPC_INTERN PetscErrorCode dvd_harm_updateproj(dvdDashboard*);
271: SLEPC_INTERN PetscErrorCode dvd_harm_conf(dvdDashboard*,dvdBlackboard*,HarmType_t,PetscBool,PetscScalar);

273: /* Internal interface routines */
274: SLEPC_INTERN PetscErrorCode EPSReset_XD(EPS);
275: SLEPC_INTERN PetscErrorCode EPSSetUp_XD(EPS);
276: SLEPC_INTERN PetscErrorCode EPSSolve_XD(EPS);
277: SLEPC_INTERN PetscErrorCode EPSComputeVectors_XD(EPS);
278: SLEPC_INTERN PetscErrorCode EPSXDSetKrylovStart_XD(EPS,PetscBool);
279: SLEPC_INTERN PetscErrorCode EPSXDGetKrylovStart_XD(EPS,PetscBool*);
280: SLEPC_INTERN PetscErrorCode EPSXDSetBlockSize_XD(EPS,PetscInt);
281: SLEPC_INTERN PetscErrorCode EPSXDGetBlockSize_XD(EPS,PetscInt*);
282: SLEPC_INTERN PetscErrorCode EPSXDSetRestart_XD(EPS,PetscInt,PetscInt);
283: SLEPC_INTERN PetscErrorCode EPSXDGetRestart_XD(EPS,PetscInt*,PetscInt*);
284: SLEPC_INTERN PetscErrorCode EPSXDGetInitialSize_XD(EPS,PetscInt*);
285: SLEPC_INTERN PetscErrorCode EPSXDSetInitialSize_XD(EPS,PetscInt);
286: SLEPC_INTERN PetscErrorCode EPSXDSetBOrth_XD(EPS,PetscBool);
287: SLEPC_INTERN PetscErrorCode EPSXDGetBOrth_XD(EPS,PetscBool*);
288: SLEPC_INTERN PetscErrorCode EPSJDGetFix_JD(EPS,PetscReal*);
289: SLEPC_INTERN PetscErrorCode EPSJDGetConstCorrectionTol_JD(EPS,PetscBool*);