Actual source code: davidson.c
slepc-3.23.1 2025-05-01
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: PetscFunctionBegin;
49: EPSCheckNotStructured(eps);
50: /* Setup EPS options and get the problem specification */
51: if (eps->nev==0) eps->nev = 1;
52: bs = data->blocksize;
53: if (bs <= 0) bs = 1;
54: if (eps->ncv!=PETSC_DETERMINE) {
55: PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of ncv must be at least nev");
56: } else if (eps->mpd!=PETSC_DETERMINE) eps->ncv = eps->mpd + eps->nev + bs;
57: else if (eps->n < 10) eps->ncv = eps->n+eps->nev+bs;
58: else if (eps->nev < 500) eps->ncv = PetscMax(eps->nev,PetscMin(eps->n-bs,PetscMax(2*eps->nev,eps->nev+15))+bs);
59: else eps->ncv = PetscMax(eps->nev,PetscMin(eps->n-bs,eps->nev+500)+bs);
60: if (eps->mpd==PETSC_DETERMINE) eps->mpd = PetscMin(eps->n,eps->ncv);
61: PetscCheck(eps->mpd<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd parameter has to be less than or equal to ncv");
62: PetscCheck(eps->mpd>=2,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd parameter has to be greater than 2");
63: if (eps->max_it == PETSC_DETERMINE) eps->max_it = PetscMax(100*eps->ncv,2*eps->n);
64: if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
65: PetscCheck(eps->nev+bs<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of ncv has to be greater than nev plus blocksize");
66: PetscCheck(!eps->trueres,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"-eps_true_residual is disabled in this solver.");
67: EPSCheckUnsupported(eps,EPS_FEATURE_REGION | EPS_FEATURE_TWOSIDED | EPS_FEATURE_THRESHOLD);
69: if (!data->minv) data->minv = (eps->n && eps->n<10)? 1: PetscMin(PetscMax(bs,6),eps->mpd/2);
70: min_size_V = data->minv;
71: PetscCheck(min_size_V+bs<=eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of minv must be less than mpd minus blocksize");
72: if (data->plusk == PETSC_DETERMINE) {
73: if (eps->problem_type == EPS_GHIEP || eps->nev+bs>eps->ncv) data->plusk = 0;
74: else data->plusk = 1;
75: }
76: if (!data->initialsize) data->initialsize = (eps->n && eps->n<10)? 1: 6;
77: initv = data->initialsize;
78: PetscCheck(eps->mpd>=initv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The initv parameter has to be less than or equal to mpd");
80: /* Change the default sigma to inf if necessary */
81: if (eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_LARGEST_REAL || eps->which == EPS_LARGEST_IMAGINARY) PetscCall(STSetDefaultShift(eps->st,PETSC_MAX_REAL));
83: /* Set up preconditioner */
84: PetscCall(STSetUp(eps->st));
86: /* Setup problem specification in dvd */
87: PetscCall(STGetNumMatrices(eps->st,&nmat));
88: PetscCall(STGetMatrix(eps->st,0,&A));
89: if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
90: PetscCall(EPSReset_XD(eps));
91: PetscCall(PetscMemzero(dvd,sizeof(dvdDashboard)));
92: dvd->A = A; dvd->B = eps->isgeneralized? B: NULL;
93: ispositive = eps->ispositive;
94: dvd->sA = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | ((ispositive && !eps->isgeneralized) ? DVD_MAT_POS_DEF: 0);
95: /* Assume -eps_hermitian means hermitian-definite in generalized problems */
96: if (!ispositive && !eps->isgeneralized && eps->ishermitian) ispositive = PETSC_TRUE;
97: if (!eps->isgeneralized) dvd->sB = DVD_MAT_IMPLICIT | DVD_MAT_HERMITIAN | DVD_MAT_IDENTITY | DVD_MAT_UNITARY | DVD_MAT_POS_DEF;
98: else dvd->sB = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | (ispositive? DVD_MAT_POS_DEF: 0);
99: ipB = (dvd->B && data->ipB && DVD_IS(dvd->sB,DVD_MAT_HERMITIAN))?PETSC_TRUE:PETSC_FALSE;
100: 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);
101: if (data->ipB && !ipB) data->ipB = PETSC_FALSE;
102: dvd->correctXnorm = (dvd->B && (DVD_IS(dvd->sB,DVD_MAT_HERMITIAN)||DVD_IS(dvd->sEP,DVD_EP_INDEFINITE)))?PETSC_TRUE:PETSC_FALSE;
103: dvd->nev = eps->nev;
104: dvd->which = eps->which;
105: dvd->withTarget = PETSC_TRUE;
106: switch (eps->which) {
107: case EPS_TARGET_MAGNITUDE:
108: case EPS_TARGET_IMAGINARY:
109: dvd->target[0] = target = eps->target;
110: dvd->target[1] = 1.0;
111: break;
112: case EPS_TARGET_REAL:
113: dvd->target[0] = PetscRealPart(target = eps->target);
114: dvd->target[1] = 1.0;
115: break;
116: case EPS_LARGEST_REAL:
117: case EPS_LARGEST_MAGNITUDE:
118: case EPS_LARGEST_IMAGINARY: /* TODO: think about this case */
119: dvd->target[0] = 1.0;
120: dvd->target[1] = target = 0.0;
121: break;
122: case EPS_SMALLEST_MAGNITUDE:
123: case EPS_SMALLEST_REAL:
124: case EPS_SMALLEST_IMAGINARY: /* TODO: think about this case */
125: dvd->target[0] = target = 0.0;
126: dvd->target[1] = 1.0;
127: break;
128: case EPS_WHICH_USER:
129: PetscCall(STGetShift(eps->st,&target));
130: dvd->target[0] = target;
131: dvd->target[1] = 1.0;
132: break;
133: case EPS_ALL:
134: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
135: default:
136: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported value of option 'which'");
137: }
138: dvd->tol = SlepcDefaultTol(eps->tol);
139: dvd->eps = eps;
141: /* Setup the extraction technique */
142: if (!eps->extraction) {
143: if (ipB || ispositive) eps->extraction = EPS_RITZ;
144: else {
145: switch (eps->which) {
146: case EPS_TARGET_REAL:
147: case EPS_TARGET_MAGNITUDE:
148: case EPS_TARGET_IMAGINARY:
149: case EPS_SMALLEST_MAGNITUDE:
150: case EPS_SMALLEST_REAL:
151: case EPS_SMALLEST_IMAGINARY:
152: eps->extraction = EPS_HARMONIC;
153: break;
154: case EPS_LARGEST_REAL:
155: case EPS_LARGEST_MAGNITUDE:
156: case EPS_LARGEST_IMAGINARY:
157: eps->extraction = EPS_HARMONIC_LARGEST;
158: break;
159: default:
160: eps->extraction = EPS_RITZ;
161: }
162: }
163: }
164: switch (eps->extraction) {
165: case EPS_RITZ: harm = DVD_HARM_NONE; break;
166: case EPS_HARMONIC: harm = DVD_HARM_RR; break;
167: case EPS_HARMONIC_RELATIVE: harm = DVD_HARM_RRR; break;
168: case EPS_HARMONIC_RIGHT: harm = DVD_HARM_REIGS; break;
169: case EPS_HARMONIC_LARGEST: harm = DVD_HARM_LEIGS; break;
170: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
171: }
173: /* Setup the type of starting subspace */
174: init = data->krylovstart? DVD_INITV_KRYLOV: DVD_INITV_CLASSIC;
176: /* Preconfigure dvd */
177: PetscCall(STGetKSP(eps->st,&ksp));
178: PetscCall(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));
180: /* Allocate memory */
181: PetscCall(EPSAllocateSolution(eps,0));
183: /* Setup orthogonalization */
184: PetscCall(EPS_SetInnerProduct(eps));
185: if (!(ipB && dvd->B)) PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
187: /* Configure dvd for a basic GD */
188: PetscCall(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));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: PetscErrorCode EPSSolve_XD(EPS eps)
193: {
194: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
195: dvdDashboard *d = &data->ddb;
196: PetscInt l,k;
198: PetscFunctionBegin;
199: PetscCall(PetscCitationsRegister(citation,&cited));
200: /* Call the starting routines */
201: PetscCall(EPSDavidsonFLCall(d->startList,d));
203: while (eps->reason == EPS_CONVERGED_ITERATING) {
205: /* Initialize V, if it is needed */
206: PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
207: if (PetscUnlikely(l == k)) PetscCall(d->initV(d));
209: /* Find the best approximated eigenpairs in V, X */
210: PetscCall(d->calcPairs(d));
212: /* Test for convergence */
213: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
214: if (eps->reason != EPS_CONVERGED_ITERATING) break;
216: /* Expand the subspace */
217: PetscCall(d->updateV(d));
219: /* Monitor progress */
220: eps->nconv = d->nconv;
221: eps->its++;
222: PetscCall(BVGetActiveColumns(d->eps->V,NULL,&k));
223: PetscCall(EPSMonitor(eps,eps->its,eps->nconv+d->npreconv,eps->eigr,eps->eigi,eps->errest,PetscMin(k,eps->nev)));
224: }
226: /* Call the ending routines */
227: PetscCall(EPSDavidsonFLCall(d->endList,d));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: PetscErrorCode EPSReset_XD(EPS eps)
232: {
233: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
234: dvdDashboard *dvd = &data->ddb;
236: PetscFunctionBegin;
237: /* Call step destructors and destroys the list */
238: PetscCall(EPSDavidsonFLCall(dvd->destroyList,dvd));
239: PetscCall(EPSDavidsonFLDestroy(&dvd->destroyList));
240: PetscCall(EPSDavidsonFLDestroy(&dvd->startList));
241: PetscCall(EPSDavidsonFLDestroy(&dvd->endList));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: PetscErrorCode EPSXDSetKrylovStart_XD(EPS eps,PetscBool krylovstart)
246: {
247: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
249: PetscFunctionBegin;
250: data->krylovstart = krylovstart;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscErrorCode EPSXDGetKrylovStart_XD(EPS eps,PetscBool *krylovstart)
255: {
256: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
258: PetscFunctionBegin;
259: *krylovstart = data->krylovstart;
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: PetscErrorCode EPSXDSetBlockSize_XD(EPS eps,PetscInt blocksize)
264: {
265: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
267: PetscFunctionBegin;
268: if (blocksize == PETSC_DEFAULT || blocksize == PETSC_DECIDE) blocksize = 1;
269: PetscCheck(blocksize>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid blocksize value, must be >0");
270: if (data->blocksize != blocksize) {
271: data->blocksize = blocksize;
272: eps->state = EPS_STATE_INITIAL;
273: }
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode EPSXDGetBlockSize_XD(EPS eps,PetscInt *blocksize)
278: {
279: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
281: PetscFunctionBegin;
282: *blocksize = data->blocksize;
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: PetscErrorCode EPSXDSetRestart_XD(EPS eps,PetscInt minv,PetscInt plusk)
287: {
288: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
290: PetscFunctionBegin;
291: if (minv == PETSC_DETERMINE) {
292: if (data->minv != 0) eps->state = EPS_STATE_INITIAL;
293: data->minv = 0;
294: } else if (minv != PETSC_CURRENT) {
295: PetscCheck(minv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid minv value, must be >0");
296: if (data->minv != minv) eps->state = EPS_STATE_INITIAL;
297: data->minv = minv;
298: }
299: if (plusk == PETSC_DETERMINE) {
300: if (data->plusk != PETSC_DETERMINE) eps->state = EPS_STATE_INITIAL;
301: data->plusk = PETSC_DETERMINE;
302: } else if (plusk != PETSC_CURRENT) {
303: PetscCheck(plusk>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid plusk value, must be >0");
304: if (data->plusk != plusk) eps->state = EPS_STATE_INITIAL;
305: data->plusk = plusk;
306: }
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: PetscErrorCode EPSXDGetRestart_XD(EPS eps,PetscInt *minv,PetscInt *plusk)
311: {
312: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
314: PetscFunctionBegin;
315: if (minv) *minv = data->minv;
316: if (plusk) *plusk = data->plusk;
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: PetscErrorCode EPSXDGetInitialSize_XD(EPS eps,PetscInt *initialsize)
321: {
322: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
324: PetscFunctionBegin;
325: *initialsize = data->initialsize;
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: PetscErrorCode EPSXDSetInitialSize_XD(EPS eps,PetscInt initialsize)
330: {
331: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
333: PetscFunctionBegin;
334: if (initialsize == PETSC_DEFAULT || initialsize == PETSC_DECIDE) initialsize = 0;
335: else PetscCheck(initialsize>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial size value, must be >0");
336: if (data->initialsize != initialsize) {
337: data->initialsize = initialsize;
338: eps->state = EPS_STATE_INITIAL;
339: }
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
344: {
345: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
347: PetscFunctionBegin;
348: data->ipB = borth;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: PetscErrorCode EPSXDGetBOrth_XD(EPS eps,PetscBool *borth)
353: {
354: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
356: PetscFunctionBegin;
357: *borth = data->ipB;
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*
362: EPSComputeVectors_XD - Compute eigenvectors from the vectors
363: provided by the eigensolver. This version is intended for solvers
364: that provide Schur vectors from the QZ decomposition. Given the partial
365: Schur decomposition OP*V=V*T, the following steps are performed:
366: 1) compute eigenvectors of (S,T): S*Z=T*Z*D
367: 2) compute eigenvectors of OP: X=V*Z
368: */
369: PetscErrorCode EPSComputeVectors_XD(EPS eps)
370: {
371: Mat X;
372: PetscBool symm;
374: PetscFunctionBegin;
375: PetscCall(PetscObjectTypeCompare((PetscObject)eps->ds,DSHEP,&symm));
376: if (symm) PetscFunctionReturn(PETSC_SUCCESS);
377: PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
379: /* V <- V * X */
380: PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
381: PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
382: PetscCall(BVMultInPlace(eps->V,X,0,eps->nconv));
383: PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }