Actual source code: davidson.c
slepc-3.14.1 2020-12-08
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: {
38: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
39: dvdDashboard *dvd = &data->ddb;
40: dvdBlackboard b;
41: PetscInt min_size_V,bs,initv,nmat;
42: Mat A,B;
43: KSP ksp;
44: PetscBool ipB,ispositive;
45: HarmType_t harm;
46: InitType_t init;
47: PetscScalar target;
50: /* Setup EPS options and get the problem specification */
51: bs = data->blocksize;
52: if (bs <= 0) bs = 1;
53: if (eps->ncv!=PETSC_DEFAULT) {
54: if (eps->ncv<eps->nev) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of ncv must be at least nev");
55: } else if (eps->mpd!=PETSC_DEFAULT) eps->ncv = eps->mpd + eps->nev + bs;
56: else if (eps->nev<500) eps->ncv = PetscMin(eps->n-bs,PetscMax(2*eps->nev,eps->nev+15))+bs;
57: else eps->ncv = PetscMin(eps->n-bs,eps->nev+500)+bs;
58: if (eps->mpd==PETSC_DEFAULT) eps->mpd = eps->ncv;
59: if (eps->mpd > eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd has to be less or equal than ncv");
60: if (eps->mpd < 2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd has to be greater than 2");
61: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100*eps->ncv,2*eps->n);
62: if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
63: if (!(eps->nev + bs <= eps->ncv)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The ncv has to be greater than nev plus blocksize");
64: if (eps->trueres) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"-eps_true_residual is temporally disable in this solver.");
65: EPSCheckUnsupported(eps,EPS_FEATURE_REGION | EPS_FEATURE_TWOSIDED);
67: EPSXDSetRestart_XD(eps,data->minv,data->plusk);
68: min_size_V = data->minv;
69: if (!min_size_V) min_size_V = PetscMin(PetscMax(bs,5),eps->mpd/2);
70: if (!(min_size_V+bs <= eps->mpd)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of minv must be less than mpd minus blocksize");
71: initv = data->initialsize;
72: if (eps->mpd < initv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The initv has to be less or equal than mpd");
74: /* Change the default sigma to inf if necessary */
75: if (eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_LARGEST_REAL || eps->which == EPS_LARGEST_IMAGINARY) {
76: STSetDefaultShift(eps->st,PETSC_MAX_REAL);
77: }
79: /* Set up preconditioner */
80: STSetUp(eps->st);
82: /* Setup problem specification in dvd */
83: STGetNumMatrices(eps->st,&nmat);
84: STGetMatrix(eps->st,0,&A);
85: if (nmat>1) { STGetMatrix(eps->st,1,&B); }
86: EPSReset_XD(eps);
87: PetscMemzero(dvd,sizeof(dvdDashboard));
88: dvd->A = A; dvd->B = eps->isgeneralized? B: NULL;
89: ispositive = eps->ispositive;
90: dvd->sA = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | ((ispositive && !eps->isgeneralized) ? DVD_MAT_POS_DEF: 0);
91: /* Asume -eps_hermitian means hermitian-definite in generalized problems */
92: if (!ispositive && !eps->isgeneralized && eps->ishermitian) ispositive = PETSC_TRUE;
93: if (!eps->isgeneralized) dvd->sB = DVD_MAT_IMPLICIT | DVD_MAT_HERMITIAN | DVD_MAT_IDENTITY | DVD_MAT_UNITARY | DVD_MAT_POS_DEF;
94: else dvd->sB = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | (ispositive? DVD_MAT_POS_DEF: 0);
95: ipB = (dvd->B && data->ipB && DVD_IS(dvd->sB,DVD_MAT_HERMITIAN))?PETSC_TRUE:PETSC_FALSE;
96: 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);
97: if (data->ipB && !ipB) data->ipB = PETSC_FALSE;
98: dvd->correctXnorm = (dvd->B && (DVD_IS(dvd->sB,DVD_MAT_HERMITIAN)||DVD_IS(dvd->sEP,DVD_EP_INDEFINITE)))?PETSC_TRUE:PETSC_FALSE;
99: dvd->nev = eps->nev;
100: dvd->which = eps->which;
101: dvd->withTarget = PETSC_TRUE;
102: switch (eps->which) {
103: case EPS_TARGET_MAGNITUDE:
104: case EPS_TARGET_IMAGINARY:
105: dvd->target[0] = target = eps->target;
106: dvd->target[1] = 1.0;
107: break;
108: case EPS_TARGET_REAL:
109: dvd->target[0] = PetscRealPart(target = eps->target);
110: dvd->target[1] = 1.0;
111: break;
112: case EPS_LARGEST_REAL:
113: case EPS_LARGEST_MAGNITUDE:
114: case EPS_LARGEST_IMAGINARY: /* TODO: think about this case */
115: dvd->target[0] = 1.0;
116: dvd->target[1] = target = 0.0;
117: break;
118: case EPS_SMALLEST_MAGNITUDE:
119: case EPS_SMALLEST_REAL:
120: case EPS_SMALLEST_IMAGINARY: /* TODO: think about this case */
121: dvd->target[0] = target = 0.0;
122: dvd->target[1] = 1.0;
123: break;
124: case EPS_WHICH_USER:
125: STGetShift(eps->st,&target);
126: dvd->target[0] = target;
127: dvd->target[1] = 1.0;
128: break;
129: case EPS_ALL:
130: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
131: default:
132: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported value of option 'which'");
133: }
134: dvd->tol = (eps->tol==PETSC_DEFAULT)? SLEPC_DEFAULT_TOL: eps->tol;
135: dvd->eps = eps;
137: /* Setup the extraction technique */
138: if (!eps->extraction) {
139: if (ipB || ispositive) eps->extraction = EPS_RITZ;
140: else {
141: switch (eps->which) {
142: case EPS_TARGET_REAL:
143: case EPS_TARGET_MAGNITUDE:
144: case EPS_TARGET_IMAGINARY:
145: case EPS_SMALLEST_MAGNITUDE:
146: case EPS_SMALLEST_REAL:
147: case EPS_SMALLEST_IMAGINARY:
148: eps->extraction = EPS_HARMONIC;
149: break;
150: case EPS_LARGEST_REAL:
151: case EPS_LARGEST_MAGNITUDE:
152: case EPS_LARGEST_IMAGINARY:
153: eps->extraction = EPS_HARMONIC_LARGEST;
154: break;
155: default:
156: eps->extraction = EPS_RITZ;
157: }
158: }
159: }
160: switch (eps->extraction) {
161: case EPS_RITZ: harm = DVD_HARM_NONE; break;
162: case EPS_HARMONIC: harm = DVD_HARM_RR; break;
163: case EPS_HARMONIC_RELATIVE: harm = DVD_HARM_RRR; break;
164: case EPS_HARMONIC_RIGHT: harm = DVD_HARM_REIGS; break;
165: case EPS_HARMONIC_LARGEST: harm = DVD_HARM_LEIGS; break;
166: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
167: }
169: /* Setup the type of starting subspace */
170: init = data->krylovstart? DVD_INITV_KRYLOV: DVD_INITV_CLASSIC;
172: /* Preconfigure dvd */
173: STGetKSP(eps->st,&ksp);
174: 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);
176: /* Allocate memory */
177: EPSAllocateSolution(eps,0);
179: /* Setup orthogonalization */
180: EPS_SetInnerProduct(eps);
181: if (!(ipB && dvd->B)) {
182: BVSetMatrix(eps->V,NULL,PETSC_FALSE);
183: }
185: /* Configure dvd for a basic GD */
186: 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);
187: return(0);
188: }
190: PetscErrorCode EPSSolve_XD(EPS eps)
191: {
192: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
193: dvdDashboard *d = &data->ddb;
194: PetscInt l,k;
198: PetscCitationsRegister(citation,&cited);
199: /* Call the starting routines */
200: EPSDavidsonFLCall(d->startList,d);
202: while (eps->reason == EPS_CONVERGED_ITERATING) {
204: /* Initialize V, if it is needed */
205: BVGetActiveColumns(d->eps->V,&l,&k);
206: if (l == k) { d->initV(d); }
208: /* Find the best approximated eigenpairs in V, X */
209: d->calcPairs(d);
211: /* Test for convergence */
212: (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
213: if (eps->reason != EPS_CONVERGED_ITERATING) break;
215: /* Expand the subspace */
216: d->updateV(d);
218: /* Monitor progress */
219: eps->nconv = d->nconv;
220: eps->its++;
221: BVGetActiveColumns(d->eps->V,NULL,&k);
222: EPSMonitor(eps,eps->its,eps->nconv+d->npreconv,eps->eigr,eps->eigi,eps->errest,PetscMin(k,eps->nev));
223: }
225: /* Call the ending routines */
226: EPSDavidsonFLCall(d->endList,d);
227: return(0);
228: }
230: PetscErrorCode EPSReset_XD(EPS eps)
231: {
232: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
233: dvdDashboard *dvd = &data->ddb;
237: /* Call step destructors and destroys the list */
238: EPSDavidsonFLCall(dvd->destroyList,dvd);
239: EPSDavidsonFLDestroy(&dvd->destroyList);
240: EPSDavidsonFLDestroy(&dvd->startList);
241: EPSDavidsonFLDestroy(&dvd->endList);
242: return(0);
243: }
245: PetscErrorCode EPSXDSetKrylovStart_XD(EPS eps,PetscBool krylovstart)
246: {
247: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
250: data->krylovstart = krylovstart;
251: return(0);
252: }
254: PetscErrorCode EPSXDGetKrylovStart_XD(EPS eps,PetscBool *krylovstart)
255: {
256: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
259: *krylovstart = data->krylovstart;
260: return(0);
261: }
263: PetscErrorCode EPSXDSetBlockSize_XD(EPS eps,PetscInt blocksize)
264: {
265: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
268: if (blocksize == PETSC_DEFAULT || blocksize == PETSC_DECIDE) blocksize = 1;
269: if (blocksize <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid blocksize value");
270: data->blocksize = blocksize;
271: return(0);
272: }
274: PetscErrorCode EPSXDGetBlockSize_XD(EPS eps,PetscInt *blocksize)
275: {
276: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
279: *blocksize = data->blocksize;
280: return(0);
281: }
283: PetscErrorCode EPSXDSetRestart_XD(EPS eps,PetscInt minv,PetscInt plusk)
284: {
285: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
288: if (minv == PETSC_DEFAULT || minv == PETSC_DECIDE) minv = 5;
289: if (minv <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid minv value");
290: if (plusk == PETSC_DEFAULT || plusk == PETSC_DECIDE) plusk = eps->problem_type == EPS_GHIEP?0:1;
291: if (plusk < 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid plusk value");
292: data->minv = minv;
293: data->plusk = plusk;
294: return(0);
295: }
297: PetscErrorCode EPSXDGetRestart_XD(EPS eps,PetscInt *minv,PetscInt *plusk)
298: {
299: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
302: if (minv) *minv = data->minv;
303: if (plusk) *plusk = data->plusk;
304: return(0);
305: }
307: PetscErrorCode EPSXDGetInitialSize_XD(EPS eps,PetscInt *initialsize)
308: {
309: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
312: *initialsize = data->initialsize;
313: return(0);
314: }
316: PetscErrorCode EPSXDSetInitialSize_XD(EPS eps,PetscInt initialsize)
317: {
318: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
321: if (initialsize == PETSC_DEFAULT || initialsize == PETSC_DECIDE) initialsize = 5;
322: if (initialsize <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial size value");
323: data->initialsize = initialsize;
324: return(0);
325: }
327: PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
328: {
329: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
332: data->ipB = borth;
333: return(0);
334: }
336: PetscErrorCode EPSXDGetBOrth_XD(EPS eps,PetscBool *borth)
337: {
338: EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
341: *borth = data->ipB;
342: return(0);
343: }
345: /*
346: EPSComputeVectors_XD - Compute eigenvectors from the vectors
347: provided by the eigensolver. This version is intended for solvers
348: that provide Schur vectors from the QZ decomposition. Given the partial
349: Schur decomposition OP*V=V*T, the following steps are performed:
350: 1) compute eigenvectors of (S,T): S*Z=T*Z*D
351: 2) compute eigenvectors of OP: X=V*Z
352: */
353: PetscErrorCode EPSComputeVectors_XD(EPS eps)
354: {
356: Mat X;
357: PetscBool symm;
360: PetscObjectTypeCompare((PetscObject)eps->ds,DSHEP,&symm);
361: if (symm) return(0);
362: DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
364: /* V <- V * X */
365: DSGetMat(eps->ds,DS_MAT_X,&X);
366: BVSetActiveColumns(eps->V,0,eps->nconv);
367: BVMultInPlace(eps->V,X,0,eps->nconv);
368: MatDestroy(&X);
369: return(0);
370: }