Actual source code: davidson.c
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: Skeleton of Davidson solver. Actual solvers are GD and JD.
13: References:
15: [1] E. Romero and J.E. Roman, "A parallel implementation of Davidson
16: methods for large-scale eigenvalue problems in SLEPc", ACM Trans.
17: Math. Software 40(2):13, 2014.
18: */
20: #include "davidson.h"
22: static PetscBool cited = PETSC_FALSE;
23: static const char citation[] =
24: "@Article{slepc-davidson,\n"
25: " author = \"E. Romero and J. E. Roman\",\n"
26: " title = \"A parallel implementation of {Davidson} methods for large-scale eigenvalue problems in {SLEPc}\",\n"
27: " journal = \"{ACM} Trans. Math. Software\",\n"
28: " volume = \"40\",\n"
29: " number = \"2\",\n"
30: " pages = \"13:1--13:29\",\n"
31: " year = \"2014,\"\n"
32: " doi = \"https://doi.org/10.1145/2543696\"\n"
33: "}\n";
35: PetscErrorCode EPSSetUp_XD(EPS eps)
36: {
37: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
38: dvdDashboard *dvd = &data->ddb;
39: dvdBlackboard b;
40: PetscInt min_size_V,bs,initv,nmat;
41: Mat A,B;
42: KSP ksp;
43: PetscBool ipB,ispositive;
44: HarmType_t harm;
45: InitType_t init;
46: PetscScalar target;
48: /* Setup EPS options and get the problem specification */
49: bs = data->blocksize;
50: if (bs <= 0) bs = 1;
51: if (eps->ncv!=PETSC_DEFAULT && eps->ncv!=PETSC_DECIDE) {
53: } else if (eps->mpd!=PETSC_DEFAULT && eps->mpd!=PETSC_DECIDE) eps->ncv = eps->mpd + eps->nev + bs;
54: else if (eps->n < 10) eps->ncv = eps->n+eps->nev+bs;
55: else if (eps->nev < 500) eps->ncv = PetscMax(eps->nev,PetscMin(eps->n-bs,PetscMax(2*eps->nev,eps->nev+15))+bs);
56: else eps->ncv = PetscMax(eps->nev,PetscMin(eps->n-bs,eps->nev+500)+bs);
57: if (eps->mpd==PETSC_DEFAULT || eps->mpd==PETSC_DECIDE) eps->mpd = PetscMin(eps->n,eps->ncv);
60: if (eps->max_it == PETSC_DEFAULT || eps->max_it == PETSC_DECIDE) eps->max_it = PetscMax(100*eps->ncv,2*eps->n);
61: if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
64: EPSCheckUnsupported(eps,EPS_FEATURE_REGION | EPS_FEATURE_TWOSIDED);
66: if (!data->minv) data->minv = (eps->n && eps->n<10)? 1: PetscMin(PetscMax(bs,6),eps->mpd/2);
67: min_size_V = data->minv;
69: if (data->plusk == PETSC_DEFAULT) {
70: if (eps->problem_type == EPS_GHIEP || eps->nev+bs>eps->ncv) data->plusk = 0;
71: else data->plusk = 1;
72: }
73: if (!data->initialsize) data->initialsize = (eps->n && eps->n<10)? 1: 6;
74: initv = data->initialsize;
77: /* Change the default sigma to inf if necessary */
78: if (eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_LARGEST_REAL || eps->which == EPS_LARGEST_IMAGINARY) STSetDefaultShift(eps->st,PETSC_MAX_REAL);
80: /* Set up preconditioner */
81: STSetUp(eps->st);
83: /* Setup problem specification in dvd */
84: STGetNumMatrices(eps->st,&nmat);
85: STGetMatrix(eps->st,0,&A);
86: if (nmat>1) STGetMatrix(eps->st,1,&B);
87: EPSReset_XD(eps);
88: PetscMemzero(dvd,sizeof(dvdDashboard));
89: dvd->A = A; dvd->B = eps->isgeneralized? B: NULL;
90: ispositive = eps->ispositive;
91: dvd->sA = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | ((ispositive && !eps->isgeneralized) ? DVD_MAT_POS_DEF: 0);
92: /* Assume -eps_hermitian means hermitian-definite in generalized problems */
93: if (!ispositive && !eps->isgeneralized && eps->ishermitian) ispositive = PETSC_TRUE;
94: if (!eps->isgeneralized) dvd->sB = DVD_MAT_IMPLICIT | DVD_MAT_HERMITIAN | DVD_MAT_IDENTITY | DVD_MAT_UNITARY | DVD_MAT_POS_DEF;
95: else dvd->sB = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | (ispositive? DVD_MAT_POS_DEF: 0);
96: ipB = (dvd->B && data->ipB && DVD_IS(dvd->sB,DVD_MAT_HERMITIAN))?PETSC_TRUE:PETSC_FALSE;
97: dvd->sEP = ((!eps->isgeneralized || (eps->isgeneralized && ipB))? DVD_EP_STD: 0) | (ispositive? DVD_EP_HERMITIAN: 0) | ((eps->problem_type == EPS_GHIEP && ipB) ? DVD_EP_INDEFINITE : 0);
98: if (data->ipB && !ipB) data->ipB = PETSC_FALSE;
99: dvd->correctXnorm = (dvd->B && (DVD_IS(dvd->sB,DVD_MAT_HERMITIAN)||DVD_IS(dvd->sEP,DVD_EP_INDEFINITE)))?PETSC_TRUE:PETSC_FALSE;
100: dvd->nev = eps->nev;
101: dvd->which = eps->which;
102: dvd->withTarget = PETSC_TRUE;
103: switch (eps->which) {
104: case EPS_TARGET_MAGNITUDE:
105: case EPS_TARGET_IMAGINARY:
106: dvd->target[0] = target = eps->target;
107: dvd->target[1] = 1.0;
108: break;
109: case EPS_TARGET_REAL:
110: dvd->target[0] = PetscRealPart(target = eps->target);
111: dvd->target[1] = 1.0;
112: break;
113: case EPS_LARGEST_REAL:
114: case EPS_LARGEST_MAGNITUDE:
115: case EPS_LARGEST_IMAGINARY: /* TODO: think about this case */
116: dvd->target[0] = 1.0;
117: dvd->target[1] = target = 0.0;
118: break;
119: case EPS_SMALLEST_MAGNITUDE:
120: case EPS_SMALLEST_REAL:
121: case EPS_SMALLEST_IMAGINARY: /* TODO: think about this case */
122: dvd->target[0] = target = 0.0;
123: dvd->target[1] = 1.0;
124: break;
125: case EPS_WHICH_USER:
126: STGetShift(eps->st,&target);
127: dvd->target[0] = target;
128: dvd->target[1] = 1.0;
129: break;
130: case EPS_ALL:
131: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
132: default:
133: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported value of option 'which'");
134: }
135: dvd->tol = SlepcDefaultTol(eps->tol);
136: dvd->eps = eps;
138: /* Setup the extraction technique */
139: if (!eps->extraction) {
140: if (ipB || ispositive) eps->extraction = EPS_RITZ;
141: else {
142: switch (eps->which) {
143: case EPS_TARGET_REAL:
144: case EPS_TARGET_MAGNITUDE:
145: case EPS_TARGET_IMAGINARY:
146: case EPS_SMALLEST_MAGNITUDE:
147: case EPS_SMALLEST_REAL:
148: case EPS_SMALLEST_IMAGINARY:
149: eps->extraction = EPS_HARMONIC;
150: break;
151: case EPS_LARGEST_REAL:
152: case EPS_LARGEST_MAGNITUDE:
153: case EPS_LARGEST_IMAGINARY:
154: eps->extraction = EPS_HARMONIC_LARGEST;
155: break;
156: default:
157: eps->extraction = EPS_RITZ;
158: }
159: }
160: }
161: switch (eps->extraction) {
162: case EPS_RITZ: harm = DVD_HARM_NONE; break;
163: case EPS_HARMONIC: harm = DVD_HARM_RR; break;
164: case EPS_HARMONIC_RELATIVE: harm = DVD_HARM_RRR; break;
165: case EPS_HARMONIC_RIGHT: harm = DVD_HARM_REIGS; break;
166: case EPS_HARMONIC_LARGEST: harm = DVD_HARM_LEIGS; break;
167: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
168: }
170: /* Setup the type of starting subspace */
171: init = data->krylovstart? DVD_INITV_KRYLOV: DVD_INITV_CLASSIC;
173: /* Preconfigure dvd */
174: STGetKSP(eps->st,&ksp);
175: dvd_schm_basic_preconf(dvd,&b,eps->mpd,min_size_V,bs,initv,PetscAbs(eps->nini),data->plusk,harm,ksp,init,eps->trackall,data->ipB,data->doubleexp);
177: /* Allocate memory */
178: EPSAllocateSolution(eps,0);
180: /* Setup orthogonalization */
181: EPS_SetInnerProduct(eps);
182: if (!(ipB && dvd->B)) BVSetMatrix(eps->V,NULL,PETSC_FALSE);
184: /* Configure dvd for a basic GD */
185: dvd_schm_basic_conf(dvd,&b,eps->mpd,min_size_V,bs,initv,PetscAbs(eps->nini),data->plusk,harm,dvd->withTarget,target,ksp,data->fix,init,eps->trackall,data->ipB,data->dynamic,data->doubleexp);
186: return 0;
187: }
189: PetscErrorCode EPSSolve_XD(EPS eps)
190: {
191: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
192: dvdDashboard *d = &data->ddb;
193: PetscInt l,k;
195: PetscCitationsRegister(citation,&cited);
196: /* Call the starting routines */
197: EPSDavidsonFLCall(d->startList,d);
199: while (eps->reason == EPS_CONVERGED_ITERATING) {
201: /* Initialize V, if it is needed */
202: BVGetActiveColumns(d->eps->V,&l,&k);
203: if (PetscUnlikely(l == k)) d->initV(d);
205: /* Find the best approximated eigenpairs in V, X */
206: d->calcPairs(d);
208: /* Test for convergence */
209: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
210: if (eps->reason != EPS_CONVERGED_ITERATING) break;
212: /* Expand the subspace */
213: d->updateV(d);
215: /* Monitor progress */
216: eps->nconv = d->nconv;
217: eps->its++;
218: BVGetActiveColumns(d->eps->V,NULL,&k);
219: EPSMonitor(eps,eps->its,eps->nconv+d->npreconv,eps->eigr,eps->eigi,eps->errest,PetscMin(k,eps->nev));
220: }
222: /* Call the ending routines */
223: EPSDavidsonFLCall(d->endList,d);
224: return 0;
225: }
227: PetscErrorCode EPSReset_XD(EPS eps)
228: {
229: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
230: dvdDashboard *dvd = &data->ddb;
232: /* Call step destructors and destroys the list */
233: EPSDavidsonFLCall(dvd->destroyList,dvd);
234: EPSDavidsonFLDestroy(&dvd->destroyList);
235: EPSDavidsonFLDestroy(&dvd->startList);
236: EPSDavidsonFLDestroy(&dvd->endList);
237: return 0;
238: }
240: PetscErrorCode EPSXDSetKrylovStart_XD(EPS eps,PetscBool krylovstart)
241: {
242: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
244: data->krylovstart = krylovstart;
245: return 0;
246: }
248: PetscErrorCode EPSXDGetKrylovStart_XD(EPS eps,PetscBool *krylovstart)
249: {
250: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
252: *krylovstart = data->krylovstart;
253: return 0;
254: }
256: PetscErrorCode EPSXDSetBlockSize_XD(EPS eps,PetscInt blocksize)
257: {
258: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
260: if (blocksize == PETSC_DEFAULT || blocksize == PETSC_DECIDE) blocksize = 1;
262: if (data->blocksize != blocksize) {
263: data->blocksize = blocksize;
264: eps->state = EPS_STATE_INITIAL;
265: }
266: return 0;
267: }
269: PetscErrorCode EPSXDGetBlockSize_XD(EPS eps,PetscInt *blocksize)
270: {
271: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
273: *blocksize = data->blocksize;
274: return 0;
275: }
277: PetscErrorCode EPSXDSetRestart_XD(EPS eps,PetscInt minv,PetscInt plusk)
278: {
279: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
281: if (minv == PETSC_DEFAULT || minv == PETSC_DECIDE) minv = 0;
283: if (plusk == PETSC_DEFAULT || plusk == PETSC_DECIDE) plusk = PETSC_DEFAULT;
285: if (data->minv != minv || data->plusk != plusk) {
286: data->minv = minv;
287: data->plusk = plusk;
288: eps->state = EPS_STATE_INITIAL;
289: }
290: return 0;
291: }
293: PetscErrorCode EPSXDGetRestart_XD(EPS eps,PetscInt *minv,PetscInt *plusk)
294: {
295: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
297: if (minv) *minv = data->minv;
298: if (plusk) *plusk = data->plusk;
299: return 0;
300: }
302: PetscErrorCode EPSXDGetInitialSize_XD(EPS eps,PetscInt *initialsize)
303: {
304: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
306: *initialsize = data->initialsize;
307: return 0;
308: }
310: PetscErrorCode EPSXDSetInitialSize_XD(EPS eps,PetscInt initialsize)
311: {
312: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
314: if (initialsize == PETSC_DEFAULT || initialsize == PETSC_DECIDE) initialsize = 0;
316: if (data->initialsize != initialsize) {
317: data->initialsize = initialsize;
318: eps->state = EPS_STATE_INITIAL;
319: }
320: return 0;
321: }
323: PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
324: {
325: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
327: data->ipB = borth;
328: return 0;
329: }
331: PetscErrorCode EPSXDGetBOrth_XD(EPS eps,PetscBool *borth)
332: {
333: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
335: *borth = data->ipB;
336: return 0;
337: }
339: /*
340: EPSComputeVectors_XD - Compute eigenvectors from the vectors
341: provided by the eigensolver. This version is intended for solvers
342: that provide Schur vectors from the QZ decomposition. Given the partial
343: Schur decomposition OP*V=V*T, the following steps are performed:
344: 1) compute eigenvectors of (S,T): S*Z=T*Z*D
345: 2) compute eigenvectors of OP: X=V*Z
346: */
347: PetscErrorCode EPSComputeVectors_XD(EPS eps)
348: {
349: Mat X;
350: PetscBool symm;
352: PetscObjectTypeCompare((PetscObject)eps->ds,DSHEP,&symm);
353: if (symm) return 0;
354: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
356: /* V <- V * X */
357: DSGetMat(eps->ds,DS_MAT_X,&X);
358: BVSetActiveColumns(eps->V,0,eps->nconv);
359: BVMultInPlace(eps->V,X,0,eps->nconv);
360: DSRestoreMat(eps->ds,DS_MAT_X,&X);
361: return 0;
362: }