Actual source code: davidson.c

slepc-3.17.0 2022-03-31
Report Typos and Errors
  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:   PetscFunctionReturn(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:   PetscFunctionReturn(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:   PetscFunctionReturn(0);
238: }

240: PetscErrorCode EPSXDSetKrylovStart_XD(EPS eps,PetscBool krylovstart)
241: {
242:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

244:   data->krylovstart = krylovstart;
245:   PetscFunctionReturn(0);
246: }

248: PetscErrorCode EPSXDGetKrylovStart_XD(EPS eps,PetscBool *krylovstart)
249: {
250:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

252:   *krylovstart = data->krylovstart;
253:   PetscFunctionReturn(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:   PetscFunctionReturn(0);
267: }

269: PetscErrorCode EPSXDGetBlockSize_XD(EPS eps,PetscInt *blocksize)
270: {
271:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

273:   *blocksize = data->blocksize;
274:   PetscFunctionReturn(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:   PetscFunctionReturn(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:   PetscFunctionReturn(0);
300: }

302: PetscErrorCode EPSXDGetInitialSize_XD(EPS eps,PetscInt *initialsize)
303: {
304:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

306:   *initialsize = data->initialsize;
307:   PetscFunctionReturn(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:   PetscFunctionReturn(0);
321: }

323: PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
324: {
325:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

327:   data->ipB = borth;
328:   PetscFunctionReturn(0);
329: }

331: PetscErrorCode EPSXDGetBOrth_XD(EPS eps,PetscBool *borth)
332: {
333:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

335:   *borth = data->ipB;
336:   PetscFunctionReturn(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) PetscFunctionReturn(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:   MatDestroy(&X);
361:   PetscFunctionReturn(0);
362: }