Actual source code: davidson.h

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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: struct _dvdDashboard;
 22: typedef PetscErrorCode (*dvdCallback)(struct _dvdDashboard*);
 23: typedef struct _dvdFunctionList {
 24:   dvdCallback f;
 25:   struct _dvdFunctionList *next;
 26: } dvdFunctionList;

 28: typedef enum {
 29:   DVD_HARM_NONE,
 30:   DVD_HARM_RR,
 31:   DVD_HARM_RRR,
 32:   DVD_HARM_REIGS,
 33:   DVD_HARM_LEIGS
 34: } HarmType_t;

 36: typedef enum {
 37:   DVD_INITV_CLASSIC,
 38:   DVD_INITV_KRYLOV
 39: } InitType_t;

 41: /*
 42:    Dashboard struct: contains the methods that will be employed and the tuning
 43:    options.
 44: */

 46: typedef struct _dvdDashboard {
 47:   /**** Function steps ****/
 48:   /* Initialize V */
 49:   PetscErrorCode (*initV)(struct _dvdDashboard*);
 50:   void *initV_data;

 52:   /* Find the approximate eigenpairs from V */
 53:   PetscErrorCode (*calcPairs)(struct _dvdDashboard*);
 54:   void *calcPairs_data;

 56:   /* Eigenpair test for convergence */
 57:   PetscBool (*testConv)(struct _dvdDashboard*,PetscScalar,PetscScalar,PetscReal,PetscReal*);
 58:   void *testConv_data;

 60:   /* Improve the selected eigenpairs */
 61:   PetscErrorCode (*improveX)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt*);
 62:   void *improveX_data;

 64:   /* Check for restarting */
 65:   PetscErrorCode (*isRestarting)(struct _dvdDashboard*,PetscBool*);
 66:   void *isRestarting_data;

 68:   /* Perform restarting */
 69:   PetscErrorCode (*restartV)(struct _dvdDashboard*);
 70:   void *restartV_data;

 72:   /* Update V */
 73:   PetscErrorCode (*updateV)(struct _dvdDashboard*);
 74:   void *updateV_data;

 76:   /**** Problem specification ****/
 77:   Mat         A,B;            /* problem matrices */
 78:   MatType_t   sA,sB;          /* matrix specifications */
 79:   EPType_t    sEP;            /* problem specifications */
 80:   PetscInt    nev;            /* number of eigenpairs */
 81:   EPSWhich    which;          /* spectrum selection */
 82:   PetscBool   withTarget;     /* if there is a target */
 83:   PetscScalar target[2];      /* target value */
 84:   PetscReal   tol;            /* tolerance */
 85:   PetscBool   correctXnorm;   /* if true, norm of X are computed */

 87:   /**** Subspaces specification ****/
 88:   PetscInt nconv;             /* number of converged eigenpairs */
 89:   PetscInt npreconv;          /* number of pairs ready to converge */

 91:   BV       W;                 /* left basis for harmonic case */
 92:   BV       AX;                /* A*V */
 93:   BV       BX;                /* B*V */
 94:   PetscInt size_D;            /* active vectors */
 95:   PetscInt max_size_proj;     /* max size projected problem */
 96:   PetscInt max_size_P;        /* max unconverged vectors in the projector */
 97:   PetscInt bs;                /* max vectors that expands the subspace every iteration */
 98:   EPS      eps;               /* connection to SLEPc */

100:   /**** Auxiliary space ****/
101:   VecPool auxV;               /* auxiliary vectors */
102:   BV      auxBV;              /* auxiliary vectors */

104:   /**** Eigenvalues and errors ****/
105:   PetscScalar *ceigr,*ceigi;  /* converged eigenvalues */
106:   PetscScalar *eigr,*eigi;    /* current eigenvalues */
107:   PetscReal   *nR;            /* residual norm */
108:   PetscReal   *real_nR;       /* original nR */
109:   PetscReal   *nX;            /* X norm */
110:   PetscReal   *real_nX;       /* original nX */
111:   PetscReal   *errest;        /* relative error eigenpairs */
112:   PetscReal   *nBds;          /* B-norms of projected problem  */

114:   /**** Shared function and variables ****/
115:   PetscErrorCode (*e_Vchanged)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt,PetscInt);
116:   void *e_Vchanged_data;
117:   PetscErrorCode (*calcpairs_residual)(struct _dvdDashboard*,PetscInt,PetscInt);
118:   PetscErrorCode (*calcpairs_selectPairs)(struct _dvdDashboard*,PetscInt);
119:   void *calcpairs_residual_data;
120:   PetscErrorCode (*improvex_precond)(struct _dvdDashboard*,PetscInt,Vec,Vec);
121:   void *improvex_precond_data;
122:   PetscErrorCode (*improvex_jd_proj_uv)(struct _dvdDashboard*,PetscInt,PetscInt,Vec*,Vec*,Vec*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscInt);
123:   PetscErrorCode (*improvex_jd_lit)(struct _dvdDashboard*,PetscInt,PetscScalar*,PetscScalar*,PetscInt*,PetscReal*);
124:   PetscErrorCode (*calcpairs_W)(struct _dvdDashboard*);
125:   void *calcpairs_W_data;
126:   PetscErrorCode (*calcpairs_proj_trans)(struct _dvdDashboard*);
127:   PetscErrorCode (*calcpairs_eigs_trans)(struct _dvdDashboard*);
128:   PetscErrorCode (*calcpairs_eig_backtrans)(struct _dvdDashboard*,PetscScalar,PetscScalar,PetscScalar*,PetscScalar*);
129:   PetscErrorCode (*calcpairs_proj_res)(struct _dvdDashboard*,PetscInt,PetscInt,Vec*);
130:   PetscErrorCode (*preTestConv)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt,PetscInt*);
131:   PetscErrorCode (*e_newIteration)(struct _dvdDashboard*);
132:   void *e_newIteration_data;

134:   dvdFunctionList *startList;  /* starting list */
135:   dvdFunctionList *endList;    /* ending list */
136:   dvdFunctionList *destroyList;/* destructor list */

138:   Mat       H,G;               /* projected problem matrices */
139:   Mat       auxM;              /* auxiliary dense matrix */
140:   PetscInt  size_MT;           /* rows in MT */

142:   PetscInt  V_tra_s;
143:   PetscInt  V_tra_e;       /* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
144:   PetscInt  V_new_s;
145:   PetscInt  V_new_e;           /* added to V the columns V_new_s:V_new_e */
146:   PetscBool W_shift;           /* if true W is shifted when vectors converge */
147: } dvdDashboard;

149: typedef struct {
150:   /*------------------------- User parameters ---------------------------*/
151:   PetscInt  blocksize;     /* block size */
152:   PetscInt  initialsize;   /* initial size of V */
153:   PetscInt  minv;          /* size of V after restarting */
154:   PetscInt  plusk;         /* keep plusk eigenvectors from the last iteration */
155:   PetscBool ipB;           /* true if B-ortho is used */
156:   PetscReal fix;           /* the fix parameter */
157:   PetscBool krylovstart;   /* true if the starting subspace is a Krylov basis */
158:   PetscBool dynamic;       /* true if dynamic stopping criterion is used */
159:   PetscBool doubleexp;     /* double expansion in GD (GD2) */

161:   /*----------------- Child objects and working data -------------------*/
162:   dvdDashboard ddb;
163: } EPS_DAVIDSON;

165: PETSC_STATIC_INLINE PetscErrorCode EPSDavidsonFLAdd(dvdFunctionList **fl,dvdCallback f)
166: {
168:   dvdFunctionList *l;

171:   PetscNew(&l);
172:   l->f = f;
173:   l->next = *fl;
174:   *fl = l;
175:   return(0);
176: }

178: PETSC_STATIC_INLINE PetscErrorCode EPSDavidsonFLCall(dvdFunctionList *fl,dvdDashboard *d)
179: {
181:   dvdFunctionList *l;

184:   for (l=fl;l;l=l->next) { (l->f)(d); }
185:   return(0);
186: }

188: PETSC_STATIC_INLINE PetscErrorCode EPSDavidsonFLDestroy(dvdFunctionList **fl)
189: {
190:   PetscErrorCode  ierr;
191:   dvdFunctionList *l,*l0;

194:   for (l=*fl;l;l=l0) {
195:     l0 = l->next;
196:     PetscFree(l);
197:   }
198:   *fl = NULL;
199:   return(0);
200: }

202: /*
203:   The blackboard configuration structure: saves information about the memory
204:   and other requirements.

206:   The starting memory structure:

208:   V           W?          AV          BV?          tKZ
209:   |-----------|-----------|-----------|------------|-------|
210:   nev+mpd     nev+mpd     scP+mpd     nev?+mpd     sP+scP
211:               scP+mpd                 scP+mpd

213:   The final memory structure considering W_shift:

215:   cX  V       cY?  W?     cAV AV      BcX? BV?     KZ  tKZ
216:   |---|-------|----|------|---|-------|----|-------|---|---|
217:   nev mpd     nev  mpd    scP mpd     nev  mpd     scP sP    <- shift
218:               scP                     scP                    <- !shift
219: */
220: typedef struct {
221:   PetscInt max_size_V;         /* max size of the searching subspace (mpd) */
222:   PetscInt max_size_X;         /* max size of X (bs) */
223:   PetscInt size_V;             /* real size of V (nev+size_P+mpd) */
224:   PetscInt max_size_oldX;      /* max size of oldX */
225:   PetscInt max_nev;            /* max number of converged pairs */
226:   PetscInt max_size_P;         /* number of computed vectors for the projector */
227:   PetscInt max_size_cP;        /* number of converged vectors in the projectors */
228:   PetscInt max_size_proj;      /* max size projected problem */
229:   PetscInt max_size_cX_proj;   /* max converged vectors in the projected problem */
230:   PetscInt state;              /* method states:
231:                                    0: preconfiguring
232:                                    1: configuring
233:                                    2: running */
234: } dvdBlackboard;

236: #define DVD_STATE_PRECONF 0
237: #define DVD_STATE_CONF 1
238: #define DVD_STATE_RUN 2

240: /* Prototypes of non-static auxiliary functions */
241: SLEPC_INTERN PetscErrorCode dvd_calcpairs_qz(dvdDashboard*,dvdBlackboard*,PetscBool,PetscBool);
242: SLEPC_INTERN PetscErrorCode dvd_improvex_gd2(dvdDashboard*,dvdBlackboard*,KSP,PetscInt);
243: SLEPC_INTERN PetscErrorCode dvd_improvex_jd(dvdDashboard*,dvdBlackboard*,KSP,PetscInt,PetscBool);
244: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard*,dvdBlackboard*);
245: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard*,dvdBlackboard*,PetscInt,PetscReal,PetscReal);
246: SLEPC_INTERN PetscErrorCode dvd_improvex_compute_X(dvdDashboard*,PetscInt,PetscInt,Vec*,PetscScalar*,PetscInt);
247: SLEPC_INTERN PetscErrorCode dvd_initV(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscBool);
248: SLEPC_INTERN PetscErrorCode dvd_orthV(BV,PetscInt,PetscInt);
249: SLEPC_INTERN PetscErrorCode dvd_schm_basic_preconf(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,HarmType_t,KSP,InitType_t,PetscBool,PetscBool,PetscBool);
250: 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);
251: SLEPC_INTERN PetscErrorCode dvd_testconv_slepc(dvdDashboard*,dvdBlackboard*);
252: SLEPC_INTERN PetscErrorCode dvd_managementV_basic(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool,PetscBool);
253: SLEPC_INTERN PetscErrorCode dvd_static_precond_PC(dvdDashboard*,dvdBlackboard*,PC);
254: SLEPC_INTERN PetscErrorCode dvd_harm_updateproj(dvdDashboard*);
255: SLEPC_INTERN PetscErrorCode dvd_harm_conf(dvdDashboard*,dvdBlackboard*,HarmType_t,PetscBool,PetscScalar);

257: /* Internal interface routines */
258: SLEPC_INTERN PetscErrorCode EPSReset_XD(EPS);
259: SLEPC_INTERN PetscErrorCode EPSSetUp_XD(EPS);
260: SLEPC_INTERN PetscErrorCode EPSSolve_XD(EPS);
261: SLEPC_INTERN PetscErrorCode EPSComputeVectors_XD(EPS);
262: SLEPC_INTERN PetscErrorCode EPSXDSetKrylovStart_XD(EPS,PetscBool);
263: SLEPC_INTERN PetscErrorCode EPSXDGetKrylovStart_XD(EPS,PetscBool*);
264: SLEPC_INTERN PetscErrorCode EPSXDSetBlockSize_XD(EPS,PetscInt);
265: SLEPC_INTERN PetscErrorCode EPSXDGetBlockSize_XD(EPS,PetscInt*);
266: SLEPC_INTERN PetscErrorCode EPSXDSetRestart_XD(EPS,PetscInt,PetscInt);
267: SLEPC_INTERN PetscErrorCode EPSXDGetRestart_XD(EPS,PetscInt*,PetscInt*);
268: SLEPC_INTERN PetscErrorCode EPSXDGetInitialSize_XD(EPS,PetscInt*);
269: SLEPC_INTERN PetscErrorCode EPSXDSetInitialSize_XD(EPS,PetscInt);
270: SLEPC_INTERN PetscErrorCode EPSXDSetBOrth_XD(EPS,PetscBool);
271: SLEPC_INTERN PetscErrorCode EPSXDGetBOrth_XD(EPS,PetscBool*);
272: SLEPC_INTERN PetscErrorCode EPSJDGetFix_JD(EPS,PetscReal*);
273: SLEPC_INTERN PetscErrorCode EPSJDGetConstCorrectionTol_JD(EPS,PetscBool*);