Actual source code: davidson.h
slepc-3.18.2 2023-01-26
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*);