Actual source code: davidson.h
slepc-3.21.1 2024-04-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: #pragma once
20: #include <slepc/private/epsimpl.h>
21: #include <slepc/private/vecimplslepc.h>
23: typedef PetscInt MatType_t;
24: #define DVD_MAT_HERMITIAN (1<<1)
25: #define DVD_MAT_NEG_DEF (1<<2)
26: #define DVD_MAT_POS_DEF (1<<3)
27: #define DVD_MAT_SINGULAR (1<<4)
28: #define DVD_MAT_COMPLEX (1<<5)
29: #define DVD_MAT_IMPLICIT (1<<6)
30: #define DVD_MAT_IDENTITY (1<<7)
31: #define DVD_MAT_DIAG (1<<8)
32: #define DVD_MAT_TRIANG (1<<9)
33: #define DVD_MAT_UTRIANG (1<<9)
34: #define DVD_MAT_LTRIANG (1<<10)
35: #define DVD_MAT_UNITARY (1<<11)
37: typedef PetscInt EPType_t;
38: #define DVD_EP_STD (1<<1)
39: #define DVD_EP_HERMITIAN (1<<2)
40: #define DVD_EP_INDEFINITE (1<<3)
42: #define DVD_IS(T,P) ((T) & (P))
43: #define DVD_ISNOT(T,P) (((T) & (P)) ^ (P))
45: struct _dvdDashboard;
46: typedef PetscErrorCode (*dvdCallback)(struct _dvdDashboard*);
47: typedef struct _dvdFunctionList {
48: dvdCallback f;
49: struct _dvdFunctionList *next;
50: } dvdFunctionList;
52: typedef enum {
53: DVD_HARM_NONE,
54: DVD_HARM_RR,
55: DVD_HARM_RRR,
56: DVD_HARM_REIGS,
57: DVD_HARM_LEIGS
58: } HarmType_t;
60: typedef enum {
61: DVD_INITV_CLASSIC,
62: DVD_INITV_KRYLOV
63: } InitType_t;
65: /*
66: Dashboard struct: contains the methods that will be employed and the tuning
67: options.
68: */
70: typedef struct _dvdDashboard {
71: /**** Function steps ****/
72: /* Initialize V */
73: PetscErrorCode (*initV)(struct _dvdDashboard*);
74: void *initV_data;
76: /* Find the approximate eigenpairs from V */
77: PetscErrorCode (*calcPairs)(struct _dvdDashboard*);
78: void *calcPairs_data;
80: /* Eigenpair test for convergence */
81: PetscBool (*testConv)(struct _dvdDashboard*,PetscScalar,PetscScalar,PetscReal,PetscReal*);
82: void *testConv_data;
84: /* Improve the selected eigenpairs */
85: PetscErrorCode (*improveX)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt*);
86: void *improveX_data;
88: /* Check for restarting */
89: PetscErrorCode (*isRestarting)(struct _dvdDashboard*,PetscBool*);
90: void *isRestarting_data;
92: /* Perform restarting */
93: PetscErrorCode (*restartV)(struct _dvdDashboard*);
94: void *restartV_data;
96: /* Update V */
97: PetscErrorCode (*updateV)(struct _dvdDashboard*);
98: void *updateV_data;
100: /**** Problem specification ****/
101: Mat A,B; /* problem matrices */
102: MatType_t sA,sB; /* matrix specifications */
103: EPType_t sEP; /* problem specifications */
104: PetscInt nev; /* number of eigenpairs */
105: EPSWhich which; /* spectrum selection */
106: PetscBool withTarget; /* if there is a target */
107: PetscScalar target[2]; /* target value */
108: PetscReal tol; /* tolerance */
109: PetscBool correctXnorm; /* if true, norm of X are computed */
111: /**** Subspaces specification ****/
112: PetscInt nconv; /* number of converged eigenpairs */
113: PetscInt npreconv; /* number of pairs ready to converge */
115: BV W; /* left basis for harmonic case */
116: BV AX; /* A*V */
117: BV BX; /* B*V */
118: PetscInt size_D; /* active vectors */
119: PetscInt max_size_proj; /* max size projected problem */
120: PetscInt max_size_P; /* max unconverged vectors in the projector */
121: PetscInt bs; /* max vectors that expands the subspace every iteration */
122: EPS eps; /* connection to SLEPc */
124: /**** Auxiliary space ****/
125: VecPool auxV; /* auxiliary vectors */
126: BV auxBV; /* auxiliary vectors */
128: /**** Eigenvalues and errors ****/
129: PetscScalar *ceigr,*ceigi; /* converged eigenvalues */
130: PetscScalar *eigr,*eigi; /* current eigenvalues */
131: PetscReal *nR; /* residual norm */
132: PetscReal *real_nR; /* original nR */
133: PetscReal *nX; /* X norm */
134: PetscReal *real_nX; /* original nX */
135: PetscReal *errest; /* relative error eigenpairs */
136: PetscReal *nBds; /* B-norms of projected problem */
138: /**** Shared function and variables ****/
139: PetscErrorCode (*e_Vchanged)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt,PetscInt);
140: void *e_Vchanged_data;
141: PetscErrorCode (*calcpairs_residual)(struct _dvdDashboard*,PetscInt,PetscInt);
142: PetscErrorCode (*calcpairs_selectPairs)(struct _dvdDashboard*,PetscInt);
143: void *calcpairs_residual_data;
144: PetscErrorCode (*improvex_precond)(struct _dvdDashboard*,PetscInt,Vec,Vec);
145: void *improvex_precond_data;
146: PetscErrorCode (*improvex_jd_proj_uv)(struct _dvdDashboard*,PetscInt,PetscInt,Vec*,Vec*,Vec*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscInt);
147: PetscErrorCode (*improvex_jd_lit)(struct _dvdDashboard*,PetscInt,PetscScalar*,PetscScalar*,PetscInt*,PetscReal*);
148: PetscErrorCode (*calcpairs_W)(struct _dvdDashboard*);
149: void *calcpairs_W_data;
150: PetscErrorCode (*calcpairs_proj_trans)(struct _dvdDashboard*);
151: PetscErrorCode (*calcpairs_eigs_trans)(struct _dvdDashboard*);
152: PetscErrorCode (*calcpairs_eig_backtrans)(struct _dvdDashboard*,PetscScalar,PetscScalar,PetscScalar*,PetscScalar*);
153: PetscErrorCode (*calcpairs_proj_res)(struct _dvdDashboard*,PetscInt,PetscInt,Vec*);
154: PetscErrorCode (*preTestConv)(struct _dvdDashboard*,PetscInt,PetscInt,PetscInt,PetscInt*);
155: PetscErrorCode (*e_newIteration)(struct _dvdDashboard*);
156: void *e_newIteration_data;
158: dvdFunctionList *startList; /* starting list */
159: dvdFunctionList *endList; /* ending list */
160: dvdFunctionList *destroyList;/* destructor list */
162: Mat H,G; /* projected problem matrices */
163: Mat auxM; /* auxiliary dense matrix */
164: PetscInt size_MT; /* rows in MT */
166: PetscInt V_tra_s;
167: PetscInt V_tra_e; /* cX <- [cX V*MT(0:V_tra_s-1)], V <- V*MT(V_tra_s:V_tra_e) */
168: PetscInt V_new_s;
169: PetscInt V_new_e; /* added to V the columns V_new_s:V_new_e */
170: PetscBool W_shift; /* if true W is shifted when vectors converge */
171: } dvdDashboard;
173: typedef struct {
174: /*------------------------- User parameters ---------------------------*/
175: PetscInt blocksize; /* block size */
176: PetscInt initialsize; /* initial size of V */
177: PetscInt minv; /* size of V after restarting */
178: PetscInt plusk; /* keep plusk eigenvectors from the last iteration */
179: PetscBool ipB; /* true if B-ortho is used */
180: PetscReal fix; /* the fix parameter */
181: PetscBool krylovstart; /* true if the starting subspace is a Krylov basis */
182: PetscBool dynamic; /* true if dynamic stopping criterion is used */
183: PetscBool doubleexp; /* double expansion in GD (GD2) */
185: /*----------------- Child objects and working data -------------------*/
186: dvdDashboard ddb;
187: } EPS_DAVIDSON;
189: static inline PetscErrorCode EPSDavidsonFLAdd(dvdFunctionList **fl,dvdCallback f)
190: {
191: dvdFunctionList *l;
193: PetscFunctionBegin;
194: PetscCall(PetscNew(&l));
195: l->f = f;
196: l->next = *fl;
197: *fl = l;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static inline PetscErrorCode EPSDavidsonFLCall(dvdFunctionList *fl,dvdDashboard *d)
202: {
203: dvdFunctionList *l;
205: PetscFunctionBegin;
206: for (l=fl;l;l=l->next) PetscCall(l->f(d));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: static inline PetscErrorCode EPSDavidsonFLDestroy(dvdFunctionList **fl)
211: {
212: dvdFunctionList *l,*l0;
214: PetscFunctionBegin;
215: for (l=*fl;l;l=l0) {
216: l0 = l->next;
217: PetscCall(PetscFree(l));
218: }
219: *fl = NULL;
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*
224: The blackboard configuration structure: saves information about the memory
225: and other requirements.
227: The starting memory structure:
229: V W? AV BV? tKZ
230: |-----------|-----------|-----------|------------|-------|
231: nev+mpd nev+mpd scP+mpd nev?+mpd sP+scP
232: scP+mpd scP+mpd
234: The final memory structure considering W_shift:
236: cX V cY? W? cAV AV BcX? BV? KZ tKZ
237: |---|-------|----|------|---|-------|----|-------|---|---|
238: nev mpd nev mpd scP mpd nev mpd scP sP <- shift
239: scP scP <- !shift
240: */
241: typedef struct {
242: PetscInt max_size_V; /* max size of the searching subspace (mpd) */
243: PetscInt max_size_X; /* max size of X (bs) */
244: PetscInt size_V; /* real size of V (nev+size_P+mpd) */
245: PetscInt max_size_oldX; /* max size of oldX */
246: PetscInt max_nev; /* max number of converged pairs */
247: PetscInt max_size_P; /* number of computed vectors for the projector */
248: PetscInt max_size_cP; /* number of converged vectors in the projectors */
249: PetscInt max_size_proj; /* max size projected problem */
250: PetscInt max_size_cX_proj; /* max converged vectors in the projected problem */
251: PetscInt state; /* method states:
252: 0: preconfiguring
253: 1: configuring
254: 2: running */
255: } dvdBlackboard;
257: #define DVD_STATE_PRECONF 0
258: #define DVD_STATE_CONF 1
259: #define DVD_STATE_RUN 2
261: /* Prototypes of non-static auxiliary functions */
262: SLEPC_INTERN PetscErrorCode dvd_calcpairs_qz(dvdDashboard*,dvdBlackboard*,PetscBool,PetscBool);
263: SLEPC_INTERN PetscErrorCode dvd_improvex_gd2(dvdDashboard*,dvdBlackboard*,KSP,PetscInt);
264: SLEPC_INTERN PetscErrorCode dvd_improvex_jd(dvdDashboard*,dvdBlackboard*,KSP,PetscInt,PetscBool);
265: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_proj_uv(dvdDashboard*,dvdBlackboard*);
266: SLEPC_INTERN PetscErrorCode dvd_improvex_jd_lit_const(dvdDashboard*,dvdBlackboard*,PetscInt,PetscReal,PetscReal);
267: SLEPC_INTERN PetscErrorCode dvd_improvex_compute_X(dvdDashboard*,PetscInt,PetscInt,Vec*,PetscScalar*,PetscInt);
268: SLEPC_INTERN PetscErrorCode dvd_initV(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscBool);
269: SLEPC_INTERN PetscErrorCode dvd_orthV(BV,PetscInt,PetscInt);
270: SLEPC_INTERN PetscErrorCode dvd_schm_basic_preconf(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,PetscInt,HarmType_t,KSP,InitType_t,PetscBool,PetscBool,PetscBool);
271: 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);
272: SLEPC_INTERN PetscErrorCode dvd_testconv_slepc(dvdDashboard*,dvdBlackboard*);
273: SLEPC_INTERN PetscErrorCode dvd_managementV_basic(dvdDashboard*,dvdBlackboard*,PetscInt,PetscInt,PetscInt,PetscInt,PetscBool,PetscBool);
274: SLEPC_INTERN PetscErrorCode dvd_static_precond_PC(dvdDashboard*,dvdBlackboard*,PC);
275: SLEPC_INTERN PetscErrorCode dvd_harm_updateproj(dvdDashboard*);
276: SLEPC_INTERN PetscErrorCode dvd_harm_conf(dvdDashboard*,dvdBlackboard*,HarmType_t,PetscBool,PetscScalar);
278: /* Internal interface routines */
279: SLEPC_INTERN PetscErrorCode EPSReset_XD(EPS);
280: SLEPC_INTERN PetscErrorCode EPSSetUp_XD(EPS);
281: SLEPC_INTERN PetscErrorCode EPSSolve_XD(EPS);
282: SLEPC_INTERN PetscErrorCode EPSComputeVectors_XD(EPS);
283: SLEPC_INTERN PetscErrorCode EPSXDSetKrylovStart_XD(EPS,PetscBool);
284: SLEPC_INTERN PetscErrorCode EPSXDGetKrylovStart_XD(EPS,PetscBool*);
285: SLEPC_INTERN PetscErrorCode EPSXDSetBlockSize_XD(EPS,PetscInt);
286: SLEPC_INTERN PetscErrorCode EPSXDGetBlockSize_XD(EPS,PetscInt*);
287: SLEPC_INTERN PetscErrorCode EPSXDSetRestart_XD(EPS,PetscInt,PetscInt);
288: SLEPC_INTERN PetscErrorCode EPSXDGetRestart_XD(EPS,PetscInt*,PetscInt*);
289: SLEPC_INTERN PetscErrorCode EPSXDGetInitialSize_XD(EPS,PetscInt*);
290: SLEPC_INTERN PetscErrorCode EPSXDSetInitialSize_XD(EPS,PetscInt);
291: SLEPC_INTERN PetscErrorCode EPSXDSetBOrth_XD(EPS,PetscBool);
292: SLEPC_INTERN PetscErrorCode EPSXDGetBOrth_XD(EPS,PetscBool*);
293: SLEPC_INTERN PetscErrorCode EPSJDGetFix_JD(EPS,PetscReal*);
294: SLEPC_INTERN PetscErrorCode EPSJDGetConstCorrectionTol_JD(EPS,PetscBool*);