Actual source code: davidson.h
slepc-3.11.2 2019-07-30
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*);