Actual source code: davidson.h
slepc-3.14.0 2020-09-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: PETSC_STATIC_INLINE PetscErrorCode EPSDavidsonFLAdd(dvdFunctionList **fl,dvdCallback f)
188: {
190: dvdFunctionList *l;
193: PetscNew(&l);
194: l->f = f;
195: l->next = *fl;
196: *fl = l;
197: return(0);
198: }
200: PETSC_STATIC_INLINE PetscErrorCode EPSDavidsonFLCall(dvdFunctionList *fl,dvdDashboard *d)
201: {
203: dvdFunctionList *l;
206: for (l=fl;l;l=l->next) { (l->f)(d); }
207: return(0);
208: }
210: PETSC_STATIC_INLINE PetscErrorCode EPSDavidsonFLDestroy(dvdFunctionList **fl)
211: {
212: PetscErrorCode ierr;
213: dvdFunctionList *l,*l0;
216: for (l=*fl;l;l=l0) {
217: l0 = l->next;
218: PetscFree(l);
219: }
220: *fl = NULL;
221: return(0);
222: }
224: /*
225: The blackboard configuration structure: saves information about the memory
226: and other requirements.
228: The starting memory structure:
230: V W? AV BV? tKZ
231: |-----------|-----------|-----------|------------|-------|
232: nev+mpd nev+mpd scP+mpd nev?+mpd sP+scP
233: scP+mpd scP+mpd
235: The final memory structure considering W_shift:
237: cX V cY? W? cAV AV BcX? BV? KZ tKZ
238: |---|-------|----|------|---|-------|----|-------|---|---|
239: nev mpd nev mpd scP mpd nev mpd scP sP <- shift
240: scP scP <- !shift
241: */
242: typedef struct {
243: PetscInt max_size_V; /* max size of the searching subspace (mpd) */
244: PetscInt max_size_X; /* max size of X (bs) */
245: PetscInt size_V; /* real size of V (nev+size_P+mpd) */
246: PetscInt max_size_oldX; /* max size of oldX */
247: PetscInt max_nev; /* max number of converged pairs */
248: PetscInt max_size_P; /* number of computed vectors for the projector */
249: PetscInt max_size_cP; /* number of converged vectors in the projectors */
250: PetscInt max_size_proj; /* max size projected problem */
251: PetscInt max_size_cX_proj; /* max converged vectors in the projected problem */
252: PetscInt state; /* method states:
253: 0: preconfiguring
254: 1: configuring
255: 2: running */
256: } dvdBlackboard;
258: #define DVD_STATE_PRECONF 0
259: #define DVD_STATE_CONF 1
260: #define DVD_STATE_RUN 2
262: /* Prototypes of non-static auxiliary functions */
263: SLEPC_INTERN PetscErrorCode dvd_calcpairs_qz(dvdDashboard*,dvdBlackboard*,PetscBool,PetscBool);
264: SLEPC_INTERN PetscErrorCode dvd_improvex_gd2(dvdDashboard*,dvdBlackboard*,KSP,PetscInt);
265: SLEPC_INTERN PetscErrorCode dvd_improvex_jd(dvdDashboard*,dvdBlackboard*,KSP,PetscInt,PetscBool);
266: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard*,dvdBlackboard*);
267: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard*,dvdBlackboard*,PetscInt,PetscReal,PetscReal);
268: SLEPC_INTERN PetscErrorCode dvd_improvex_compute_X(dvdDashboard*,PetscInt,PetscInt,Vec*,PetscScalar*,PetscInt);
269: SLEPC_INTERN PetscErrorCode dvd_initV(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscBool);
270: SLEPC_INTERN PetscErrorCode dvd_orthV(BV,PetscInt,PetscInt);
271: SLEPC_INTERN PetscErrorCode dvd_schm_basic_preconf(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,HarmType_t,KSP,InitType_t,PetscBool,PetscBool,PetscBool);
272: 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);
273: SLEPC_INTERN PetscErrorCode dvd_testconv_slepc(dvdDashboard*,dvdBlackboard*);
274: SLEPC_INTERN PetscErrorCode dvd_managementV_basic(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool,PetscBool);
275: SLEPC_INTERN PetscErrorCode dvd_static_precond_PC(dvdDashboard*,dvdBlackboard*,PC);
276: SLEPC_INTERN PetscErrorCode dvd_harm_updateproj(dvdDashboard*);
277: SLEPC_INTERN PetscErrorCode dvd_harm_conf(dvdDashboard*,dvdBlackboard*,HarmType_t,PetscBool,PetscScalar);
279: /* Internal interface routines */
280: SLEPC_INTERN PetscErrorCode EPSReset_XD(EPS);
281: SLEPC_INTERN PetscErrorCode EPSSetUp_XD(EPS);
282: SLEPC_INTERN PetscErrorCode EPSSolve_XD(EPS);
283: SLEPC_INTERN PetscErrorCode EPSComputeVectors_XD(EPS);
284: SLEPC_INTERN PetscErrorCode EPSXDSetKrylovStart_XD(EPS,PetscBool);
285: SLEPC_INTERN PetscErrorCode EPSXDGetKrylovStart_XD(EPS,PetscBool*);
286: SLEPC_INTERN PetscErrorCode EPSXDSetBlockSize_XD(EPS,PetscInt);
287: SLEPC_INTERN PetscErrorCode EPSXDGetBlockSize_XD(EPS,PetscInt*);
288: SLEPC_INTERN PetscErrorCode EPSXDSetRestart_XD(EPS,PetscInt,PetscInt);
289: SLEPC_INTERN PetscErrorCode EPSXDGetRestart_XD(EPS,PetscInt*,PetscInt*);
290: SLEPC_INTERN PetscErrorCode EPSXDGetInitialSize_XD(EPS,PetscInt*);
291: SLEPC_INTERN PetscErrorCode EPSXDSetInitialSize_XD(EPS,PetscInt);
292: SLEPC_INTERN PetscErrorCode EPSXDSetBOrth_XD(EPS,PetscBool);
293: SLEPC_INTERN PetscErrorCode EPSXDGetBOrth_XD(EPS,PetscBool*);
294: SLEPC_INTERN PetscErrorCode EPSJDGetFix_JD(EPS,PetscReal*);
295: SLEPC_INTERN PetscErrorCode EPSJDGetConstCorrectionTol_JD(EPS,PetscBool*);