Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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.
12 :
13 : References:
14 :
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 : */
19 :
20 : #include "davidson.h"
21 :
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";
34 :
35 95 : PetscErrorCode EPSSetUp_XD(EPS eps)
36 : {
37 95 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
38 95 : dvdDashboard *dvd = &data->ddb;
39 95 : dvdBlackboard b;
40 95 : PetscInt min_size_V,bs,initv,nmat;
41 95 : Mat A,B;
42 95 : KSP ksp;
43 95 : PetscBool ipB,ispositive;
44 95 : HarmType_t harm;
45 95 : InitType_t init;
46 95 : PetscScalar target;
47 :
48 95 : PetscFunctionBegin;
49 95 : EPSCheckNotStructured(eps);
50 : /* Setup EPS options and get the problem specification */
51 95 : if (eps->nev==0) eps->nev = 1;
52 95 : bs = data->blocksize;
53 95 : if (bs <= 0) bs = 1;
54 95 : if (eps->ncv!=PETSC_DETERMINE) {
55 41 : PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of ncv must be at least nev");
56 54 : } else if (eps->mpd!=PETSC_DETERMINE) eps->ncv = eps->mpd + eps->nev + bs;
57 53 : else if (eps->n < 10) eps->ncv = eps->n+eps->nev+bs;
58 53 : else if (eps->nev < 500) eps->ncv = PetscMax(eps->nev,PetscMin(eps->n-bs,PetscMax(2*eps->nev,eps->nev+15))+bs);
59 0 : else eps->ncv = PetscMax(eps->nev,PetscMin(eps->n-bs,eps->nev+500)+bs);
60 95 : if (eps->mpd==PETSC_DETERMINE) eps->mpd = PetscMin(eps->n,eps->ncv);
61 95 : PetscCheck(eps->mpd<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd parameter has to be less than or equal to ncv");
62 95 : PetscCheck(eps->mpd>=2,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd parameter has to be greater than 2");
63 95 : if (eps->max_it == PETSC_DETERMINE) eps->max_it = PetscMax(100*eps->ncv,2*eps->n);
64 95 : if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
65 95 : 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 95 : PetscCheck(!eps->trueres,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"-eps_true_residual is disabled in this solver.");
67 95 : EPSCheckUnsupported(eps,EPS_FEATURE_REGION | EPS_FEATURE_TWOSIDED | EPS_FEATURE_THRESHOLD);
68 :
69 95 : if (!data->minv) data->minv = (eps->n && eps->n<10)? 1: PetscMin(PetscMax(bs,6),eps->mpd/2);
70 95 : min_size_V = data->minv;
71 95 : 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 95 : if (data->plusk == PETSC_DETERMINE) {
73 62 : if (eps->problem_type == EPS_GHIEP || eps->nev+bs>eps->ncv) data->plusk = 0;
74 60 : else data->plusk = 1;
75 : }
76 157 : if (!data->initialsize) data->initialsize = (eps->n && eps->n<10)? 1: 6;
77 95 : initv = data->initialsize;
78 95 : PetscCheck(eps->mpd>=initv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The initv parameter has to be less than or equal to mpd");
79 :
80 : /* Change the default sigma to inf if necessary */
81 95 : if (eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_LARGEST_REAL || eps->which == EPS_LARGEST_IMAGINARY) PetscCall(STSetDefaultShift(eps->st,PETSC_MAX_REAL));
82 :
83 : /* Set up preconditioner */
84 95 : PetscCall(STSetUp(eps->st));
85 :
86 : /* Setup problem specification in dvd */
87 95 : PetscCall(STGetNumMatrices(eps->st,&nmat));
88 95 : PetscCall(STGetMatrix(eps->st,0,&A));
89 95 : if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
90 95 : PetscCall(EPSReset_XD(eps));
91 95 : PetscCall(PetscMemzero(dvd,sizeof(dvdDashboard)));
92 95 : dvd->A = A; dvd->B = eps->isgeneralized? B: NULL;
93 95 : ispositive = eps->ispositive;
94 118 : 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 95 : if (!ispositive && !eps->isgeneralized && eps->ishermitian) ispositive = PETSC_TRUE;
97 95 : if (!eps->isgeneralized) dvd->sB = DVD_MAT_IMPLICIT | DVD_MAT_HERMITIAN | DVD_MAT_IDENTITY | DVD_MAT_UNITARY | DVD_MAT_POS_DEF;
98 24 : else dvd->sB = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | (ispositive? DVD_MAT_POS_DEF: 0);
99 95 : ipB = (dvd->B && data->ipB && DVD_IS(dvd->sB,DVD_MAT_HERMITIAN))?PETSC_TRUE:PETSC_FALSE;
100 124 : 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 95 : if (data->ipB && !ipB) data->ipB = PETSC_FALSE;
102 95 : dvd->correctXnorm = (dvd->B && (DVD_IS(dvd->sB,DVD_MAT_HERMITIAN)||DVD_IS(dvd->sEP,DVD_EP_INDEFINITE)))?PETSC_TRUE:PETSC_FALSE;
103 95 : dvd->nev = eps->nev;
104 95 : dvd->which = eps->which;
105 95 : dvd->withTarget = PETSC_TRUE;
106 95 : switch (eps->which) {
107 15 : case EPS_TARGET_MAGNITUDE:
108 : case EPS_TARGET_IMAGINARY:
109 15 : dvd->target[0] = target = eps->target;
110 15 : dvd->target[1] = 1.0;
111 15 : break;
112 0 : case EPS_TARGET_REAL:
113 0 : dvd->target[0] = PetscRealPart(target = eps->target);
114 0 : dvd->target[1] = 1.0;
115 0 : break;
116 46 : case EPS_LARGEST_REAL:
117 : case EPS_LARGEST_MAGNITUDE:
118 : case EPS_LARGEST_IMAGINARY: /* TODO: think about this case */
119 46 : dvd->target[0] = 1.0;
120 46 : dvd->target[1] = target = 0.0;
121 46 : break;
122 30 : case EPS_SMALLEST_MAGNITUDE:
123 : case EPS_SMALLEST_REAL:
124 : case EPS_SMALLEST_IMAGINARY: /* TODO: think about this case */
125 30 : dvd->target[0] = target = 0.0;
126 30 : dvd->target[1] = 1.0;
127 30 : break;
128 4 : case EPS_WHICH_USER:
129 4 : PetscCall(STGetShift(eps->st,&target));
130 4 : dvd->target[0] = target;
131 4 : dvd->target[1] = 1.0;
132 4 : break;
133 0 : case EPS_ALL:
134 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
135 0 : default:
136 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported value of option 'which'");
137 : }
138 95 : dvd->tol = SlepcDefaultTol(eps->tol);
139 95 : dvd->eps = eps;
140 :
141 : /* Setup the extraction technique */
142 95 : if (!eps->extraction) {
143 86 : if (ipB || ispositive) eps->extraction = EPS_RITZ;
144 : else {
145 17 : switch (eps->which) {
146 6 : 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 6 : eps->extraction = EPS_HARMONIC;
153 6 : break;
154 8 : case EPS_LARGEST_REAL:
155 : case EPS_LARGEST_MAGNITUDE:
156 : case EPS_LARGEST_IMAGINARY:
157 8 : eps->extraction = EPS_HARMONIC_LARGEST;
158 8 : break;
159 3 : default:
160 3 : eps->extraction = EPS_RITZ;
161 : }
162 9 : }
163 : }
164 95 : switch (eps->extraction) {
165 : case EPS_RITZ: harm = DVD_HARM_NONE; break;
166 13 : case EPS_HARMONIC: harm = DVD_HARM_RR; break;
167 0 : case EPS_HARMONIC_RELATIVE: harm = DVD_HARM_RRR; break;
168 0 : case EPS_HARMONIC_RIGHT: harm = DVD_HARM_REIGS; break;
169 10 : case EPS_HARMONIC_LARGEST: harm = DVD_HARM_LEIGS; break;
170 0 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
171 : }
172 :
173 : /* Setup the type of starting subspace */
174 95 : init = data->krylovstart? DVD_INITV_KRYLOV: DVD_INITV_CLASSIC;
175 :
176 : /* Preconfigure dvd */
177 95 : PetscCall(STGetKSP(eps->st,&ksp));
178 95 : 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));
179 :
180 : /* Allocate memory */
181 95 : PetscCall(EPSAllocateSolution(eps,0));
182 :
183 : /* Setup orthogonalization */
184 95 : PetscCall(EPS_SetInnerProduct(eps));
185 95 : if (!(ipB && dvd->B)) PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
186 :
187 : /* Configure dvd for a basic GD */
188 95 : 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 95 : PetscFunctionReturn(PETSC_SUCCESS);
190 : }
191 :
192 95 : PetscErrorCode EPSSolve_XD(EPS eps)
193 : {
194 95 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
195 95 : dvdDashboard *d = &data->ddb;
196 95 : PetscInt l,k;
197 :
198 95 : PetscFunctionBegin;
199 95 : PetscCall(PetscCitationsRegister(citation,&cited));
200 : /* Call the starting routines */
201 95 : PetscCall(EPSDavidsonFLCall(d->startList,d));
202 :
203 7946 : while (eps->reason == EPS_CONVERGED_ITERATING) {
204 :
205 : /* Initialize V, if it is needed */
206 7946 : PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
207 7946 : if (PetscUnlikely(l == k)) PetscCall(d->initV(d));
208 :
209 : /* Find the best approximated eigenpairs in V, X */
210 7946 : PetscCall(d->calcPairs(d));
211 :
212 : /* Test for convergence */
213 7946 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
214 7946 : if (eps->reason != EPS_CONVERGED_ITERATING) break;
215 :
216 : /* Expand the subspace */
217 7851 : PetscCall(d->updateV(d));
218 :
219 : /* Monitor progress */
220 7851 : eps->nconv = d->nconv;
221 7851 : eps->its++;
222 7851 : PetscCall(BVGetActiveColumns(d->eps->V,NULL,&k));
223 7946 : PetscCall(EPSMonitor(eps,eps->its,eps->nconv+d->npreconv,eps->eigr,eps->eigi,eps->errest,PetscMin(k,eps->nev)));
224 : }
225 :
226 : /* Call the ending routines */
227 95 : PetscCall(EPSDavidsonFLCall(d->endList,d));
228 95 : PetscFunctionReturn(PETSC_SUCCESS);
229 : }
230 :
231 165 : PetscErrorCode EPSReset_XD(EPS eps)
232 : {
233 165 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
234 165 : dvdDashboard *dvd = &data->ddb;
235 :
236 165 : PetscFunctionBegin;
237 : /* Call step destructors and destroys the list */
238 165 : PetscCall(EPSDavidsonFLCall(dvd->destroyList,dvd));
239 165 : PetscCall(EPSDavidsonFLDestroy(&dvd->destroyList));
240 165 : PetscCall(EPSDavidsonFLDestroy(&dvd->startList));
241 165 : PetscCall(EPSDavidsonFLDestroy(&dvd->endList));
242 165 : PetscFunctionReturn(PETSC_SUCCESS);
243 : }
244 :
245 4 : PetscErrorCode EPSXDSetKrylovStart_XD(EPS eps,PetscBool krylovstart)
246 : {
247 4 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
248 :
249 4 : PetscFunctionBegin;
250 4 : data->krylovstart = krylovstart;
251 4 : PetscFunctionReturn(PETSC_SUCCESS);
252 : }
253 :
254 60 : PetscErrorCode EPSXDGetKrylovStart_XD(EPS eps,PetscBool *krylovstart)
255 : {
256 60 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
257 :
258 60 : PetscFunctionBegin;
259 60 : *krylovstart = data->krylovstart;
260 60 : PetscFunctionReturn(PETSC_SUCCESS);
261 : }
262 :
263 3 : PetscErrorCode EPSXDSetBlockSize_XD(EPS eps,PetscInt blocksize)
264 : {
265 3 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
266 :
267 3 : PetscFunctionBegin;
268 3 : if (blocksize == PETSC_DEFAULT || blocksize == PETSC_DECIDE) blocksize = 1;
269 3 : PetscCheck(blocksize>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid blocksize value, must be >0");
270 3 : if (data->blocksize != blocksize) {
271 3 : data->blocksize = blocksize;
272 3 : eps->state = EPS_STATE_INITIAL;
273 : }
274 3 : PetscFunctionReturn(PETSC_SUCCESS);
275 : }
276 :
277 60 : PetscErrorCode EPSXDGetBlockSize_XD(EPS eps,PetscInt *blocksize)
278 : {
279 60 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
280 :
281 60 : PetscFunctionBegin;
282 60 : *blocksize = data->blocksize;
283 60 : PetscFunctionReturn(PETSC_SUCCESS);
284 : }
285 :
286 3 : PetscErrorCode EPSXDSetRestart_XD(EPS eps,PetscInt minv,PetscInt plusk)
287 : {
288 3 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
289 :
290 3 : PetscFunctionBegin;
291 3 : if (minv == PETSC_DETERMINE) {
292 0 : if (data->minv != 0) eps->state = EPS_STATE_INITIAL;
293 0 : data->minv = 0;
294 3 : } else if (minv != PETSC_CURRENT) {
295 3 : PetscCheck(minv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid minv value, must be >0");
296 3 : if (data->minv != minv) eps->state = EPS_STATE_INITIAL;
297 3 : data->minv = minv;
298 : }
299 3 : if (plusk == PETSC_DETERMINE) {
300 1 : if (data->plusk != PETSC_DETERMINE) eps->state = EPS_STATE_INITIAL;
301 1 : data->plusk = PETSC_DETERMINE;
302 2 : } else if (plusk != PETSC_CURRENT) {
303 2 : PetscCheck(plusk>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid plusk value, must be >0");
304 2 : if (data->plusk != plusk) eps->state = EPS_STATE_INITIAL;
305 2 : data->plusk = plusk;
306 : }
307 3 : PetscFunctionReturn(PETSC_SUCCESS);
308 : }
309 :
310 60 : PetscErrorCode EPSXDGetRestart_XD(EPS eps,PetscInt *minv,PetscInt *plusk)
311 : {
312 60 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
313 :
314 60 : PetscFunctionBegin;
315 60 : if (minv) *minv = data->minv;
316 60 : if (plusk) *plusk = data->plusk;
317 60 : PetscFunctionReturn(PETSC_SUCCESS);
318 : }
319 :
320 59 : PetscErrorCode EPSXDGetInitialSize_XD(EPS eps,PetscInt *initialsize)
321 : {
322 59 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
323 :
324 59 : PetscFunctionBegin;
325 59 : *initialsize = data->initialsize;
326 59 : PetscFunctionReturn(PETSC_SUCCESS);
327 : }
328 :
329 2 : PetscErrorCode EPSXDSetInitialSize_XD(EPS eps,PetscInt initialsize)
330 : {
331 2 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
332 :
333 2 : PetscFunctionBegin;
334 2 : if (initialsize == PETSC_DEFAULT || initialsize == PETSC_DECIDE) initialsize = 0;
335 2 : else PetscCheck(initialsize>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial size value, must be >0");
336 2 : if (data->initialsize != initialsize) {
337 2 : data->initialsize = initialsize;
338 2 : eps->state = EPS_STATE_INITIAL;
339 : }
340 2 : PetscFunctionReturn(PETSC_SUCCESS);
341 : }
342 :
343 2 : PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
344 : {
345 2 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
346 :
347 2 : PetscFunctionBegin;
348 2 : data->ipB = borth;
349 2 : PetscFunctionReturn(PETSC_SUCCESS);
350 : }
351 :
352 60 : PetscErrorCode EPSXDGetBOrth_XD(EPS eps,PetscBool *borth)
353 : {
354 60 : EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
355 :
356 60 : PetscFunctionBegin;
357 60 : *borth = data->ipB;
358 60 : PetscFunctionReturn(PETSC_SUCCESS);
359 : }
360 :
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 91 : PetscErrorCode EPSComputeVectors_XD(EPS eps)
370 : {
371 91 : Mat X;
372 91 : PetscBool symm;
373 :
374 91 : PetscFunctionBegin;
375 91 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->ds,DSHEP,&symm));
376 91 : if (symm) PetscFunctionReturn(PETSC_SUCCESS);
377 26 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
378 :
379 : /* V <- V * X */
380 26 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
381 26 : PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
382 26 : PetscCall(BVMultInPlace(eps->V,X,0,eps->nconv));
383 26 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
384 26 : PetscFunctionReturn(PETSC_SUCCESS);
385 : }
|
