Actual source code: davidson.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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.
 12: */

 14: #include "davidson.h"

 16: static PetscBool  cited = PETSC_FALSE;
 17: static const char citation[] =
 18:   "@Article{slepc-davidson,\n"
 19:   "   author = \"E. Romero and J. E. Roman\",\n"
 20:   "   title = \"A parallel implementation of {Davidson} methods for large-scale eigenvalue problems in {SLEPc}\",\n"
 21:   "   journal = \"{ACM} Trans. Math. Software\",\n"
 22:   "   volume = \"40\",\n"
 23:   "   number = \"2\",\n"
 24:   "   pages = \"13:1--13:29\",\n"
 25:   "   year = \"2014,\"\n"
 26:   "   doi = \"https://doi.org/10.1145/2543696\"\n"
 27:   "}\n";

 29: PetscErrorCode EPSSetUp_XD(EPS eps)
 30: {
 32:   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;
 33:   dvdDashboard   *dvd = &data->ddb;
 34:   dvdBlackboard  b;
 35:   PetscInt       min_size_V,bs,initv,nmat;
 36:   Mat            A,B;
 37:   KSP            ksp;
 38:   PetscBool      ipB,ispositive;
 39:   HarmType_t     harm;
 40:   InitType_t     init;
 41:   PetscScalar    target;

 44:   /* Setup EPS options and get the problem specification */
 45:   bs = data->blocksize;
 46:   if (bs <= 0) bs = 1;
 47:   if (eps->ncv) {
 48:     if (eps->ncv<eps->nev) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of ncv must be at least nev");
 49:   } else if (eps->mpd) eps->ncv = eps->mpd + eps->nev + bs;
 50:   else if (eps->nev<500) eps->ncv = PetscMin(eps->n-bs,PetscMax(2*eps->nev,eps->nev+15))+bs;
 51:   else eps->ncv = PetscMin(eps->n-bs,eps->nev+500)+bs;
 52:   if (!eps->mpd) eps->mpd = eps->ncv;
 53:   if (eps->mpd > eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd has to be less or equal than ncv");
 54:   if (eps->mpd < 2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd has to be greater than 2");
 55:   if (!eps->max_it) eps->max_it = PetscMax(100*eps->ncv,2*eps->n);
 56:   if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
 57:   if (eps->ishermitian && (eps->which==EPS_LARGEST_IMAGINARY || eps->which==EPS_SMALLEST_IMAGINARY)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Wrong value of eps->which");
 58:   if (!(eps->nev + bs <= eps->ncv)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The ncv has to be greater than nev plus blocksize");
 59:   if (eps->trueres) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"-eps_true_residual is temporally disable in this solver.");

 61:   EPSXDSetRestart_XD(eps,data->minv,data->plusk);
 62:   min_size_V = data->minv;
 63:   if (!min_size_V) min_size_V = PetscMin(PetscMax(bs,5),eps->mpd/2);
 64:   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");
 65:   initv = data->initialsize;
 66:   if (eps->mpd < initv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The initv has to be less or equal than mpd");

 68:   /* Change the default sigma to inf if necessary */
 69:   if (eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_LARGEST_REAL || eps->which == EPS_LARGEST_IMAGINARY) {
 70:     STSetDefaultShift(eps->st,PETSC_MAX_REAL);
 71:   }

 73:   /* Set up preconditioner */
 74:   STSetUp(eps->st);

 76:   /* Setup problem specification in dvd */
 77:   STGetNumMatrices(eps->st,&nmat);
 78:   STGetMatrix(eps->st,0,&A);
 79:   if (nmat>1) { STGetMatrix(eps->st,1,&B); }
 80:   EPSReset_XD(eps);
 81:   PetscMemzero(dvd,sizeof(dvdDashboard));
 82:   dvd->A = A; dvd->B = eps->isgeneralized? B: NULL;
 83:   ispositive = eps->ispositive;
 84:   dvd->sA = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | ((ispositive && !eps->isgeneralized) ? DVD_MAT_POS_DEF: 0);
 85:   /* Asume -eps_hermitian means hermitian-definite in generalized problems */
 86:   if (!ispositive && !eps->isgeneralized && eps->ishermitian) ispositive = PETSC_TRUE;
 87:   if (!eps->isgeneralized) dvd->sB = DVD_MAT_IMPLICIT | DVD_MAT_HERMITIAN | DVD_MAT_IDENTITY | DVD_MAT_UNITARY | DVD_MAT_POS_DEF;
 88:   else dvd->sB = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | (ispositive? DVD_MAT_POS_DEF: 0);
 89:   ipB = (dvd->B && data->ipB && DVD_IS(dvd->sB,DVD_MAT_HERMITIAN))?PETSC_TRUE:PETSC_FALSE;
 90:   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);
 91:   if (data->ipB && !ipB) data->ipB = PETSC_FALSE;
 92:   dvd->correctXnorm = (dvd->B && (DVD_IS(dvd->sB,DVD_MAT_HERMITIAN)||DVD_IS(dvd->sEP,DVD_EP_INDEFINITE)))?PETSC_TRUE:PETSC_FALSE;
 93:   dvd->nev        = eps->nev;
 94:   dvd->which      = eps->which;
 95:   dvd->withTarget = PETSC_TRUE;
 96:   switch (eps->which) {
 97:     case EPS_TARGET_MAGNITUDE:
 98:     case EPS_TARGET_IMAGINARY:
 99:       dvd->target[0] = target = eps->target;
100:       dvd->target[1] = 1.0;
101:       break;
102:     case EPS_TARGET_REAL:
103:       dvd->target[0] = PetscRealPart(target = eps->target);
104:       dvd->target[1] = 1.0;
105:       break;
106:     case EPS_LARGEST_REAL:
107:     case EPS_LARGEST_MAGNITUDE:
108:     case EPS_LARGEST_IMAGINARY: /* TODO: think about this case */
109:       dvd->target[0] = 1.0;
110:       dvd->target[1] = target = 0.0;
111:       break;
112:     case EPS_SMALLEST_MAGNITUDE:
113:     case EPS_SMALLEST_REAL:
114:     case EPS_SMALLEST_IMAGINARY: /* TODO: think about this case */
115:       dvd->target[0] = target = 0.0;
116:       dvd->target[1] = 1.0;
117:       break;
118:     case EPS_WHICH_USER:
119:       STGetShift(eps->st,&target);
120:       dvd->target[0] = target;
121:       dvd->target[1] = 1.0;
122:       break;
123:     case EPS_ALL:
124:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported option: which == EPS_ALL");
125:       break;
126:     default:
127:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported value of option 'which'");
128:   }
129:   dvd->tol = (eps->tol==PETSC_DEFAULT)? SLEPC_DEFAULT_TOL: eps->tol;
130:   dvd->eps = eps;

132:   /* Setup the extraction technique */
133:   if (!eps->extraction) {
134:     if (ipB || ispositive) eps->extraction = EPS_RITZ;
135:     else {
136:       switch (eps->which) {
137:         case EPS_TARGET_REAL:
138:         case EPS_TARGET_MAGNITUDE:
139:         case EPS_TARGET_IMAGINARY:
140:         case EPS_SMALLEST_MAGNITUDE:
141:         case EPS_SMALLEST_REAL:
142:         case EPS_SMALLEST_IMAGINARY:
143:           eps->extraction = EPS_HARMONIC;
144:           break;
145:         case EPS_LARGEST_REAL:
146:         case EPS_LARGEST_MAGNITUDE:
147:         case EPS_LARGEST_IMAGINARY:
148:           eps->extraction = EPS_HARMONIC_LARGEST;
149:           break;
150:         default:
151:           eps->extraction = EPS_RITZ;
152:       }
153:     }
154:   }
155:   switch (eps->extraction) {
156:     case EPS_RITZ:              harm = DVD_HARM_NONE; break;
157:     case EPS_HARMONIC:          harm = DVD_HARM_RR; break;
158:     case EPS_HARMONIC_RELATIVE: harm = DVD_HARM_RRR; break;
159:     case EPS_HARMONIC_RIGHT:    harm = DVD_HARM_REIGS; break;
160:     case EPS_HARMONIC_LARGEST:  harm = DVD_HARM_LEIGS; break;
161:     default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
162:   }

164:   /* Setup the type of starting subspace */
165:   init = data->krylovstart? DVD_INITV_KRYLOV: DVD_INITV_CLASSIC;

167:   /* Preconfigure dvd */
168:   STGetKSP(eps->st,&ksp);
169:   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);

171:   /* Allocate memory */
172:   EPSAllocateSolution(eps,0);

174:   /* Setup orthogonalization */
175:   EPS_SetInnerProduct(eps);
176:   if (!(ipB && dvd->B)) {
177:     BVSetMatrix(eps->V,NULL,PETSC_FALSE);
178:   }

180:   /* Configure dvd for a basic GD */
181:   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);
182:   return(0);
183: }

185: PetscErrorCode EPSSolve_XD(EPS eps)
186: {
187:   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;
188:   dvdDashboard   *d = &data->ddb;
189:   PetscInt       l,k;

193:   PetscCitationsRegister(citation,&cited);
194:   /* Call the starting routines */
195:   EPSDavidsonFLCall(d->startList,d);

197:   while (eps->reason == EPS_CONVERGED_ITERATING) {

199:     /* Initialize V, if it is needed */
200:     BVGetActiveColumns(d->eps->V,&l,&k);
201:     if (l == k) { d->initV(d); }

203:     /* Find the best approximated eigenpairs in V, X */
204:     d->calcPairs(d);

206:     /* Test for convergence */
207:     (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);
208:     if (eps->reason != EPS_CONVERGED_ITERATING) break;

210:     /* Expand the subspace */
211:     d->updateV(d);

213:     /* Monitor progress */
214:     eps->nconv = d->nconv;
215:     eps->its++;
216:     BVGetActiveColumns(d->eps->V,NULL,&k);
217:     EPSMonitor(eps,eps->its,eps->nconv+d->npreconv,eps->eigr,eps->eigi,eps->errest,PetscMin(k,eps->nev));
218:   }

220:   /* Call the ending routines */
221:   EPSDavidsonFLCall(d->endList,d);
222:   return(0);
223: }

225: PetscErrorCode EPSReset_XD(EPS eps)
226: {
227:   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;
228:   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;

245:   data->krylovstart = krylovstart;
246:   return(0);
247: }

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

254:   *krylovstart = data->krylovstart;
255:   return(0);
256: }

258: PetscErrorCode EPSXDSetBlockSize_XD(EPS eps,PetscInt blocksize)
259: {
260:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

263:   if (blocksize == PETSC_DEFAULT || blocksize == PETSC_DECIDE) blocksize = 1;
264:   if (blocksize <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid blocksize value");
265:   data->blocksize = blocksize;
266:   return(0);
267: }

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

274:   *blocksize = data->blocksize;
275:   return(0);
276: }

278: PetscErrorCode EPSXDSetRestart_XD(EPS eps,PetscInt minv,PetscInt plusk)
279: {
280:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

283:   if (minv == PETSC_DEFAULT || minv == PETSC_DECIDE) minv = 5;
284:   if (minv <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid minv value");
285:   if (plusk == PETSC_DEFAULT || plusk == PETSC_DECIDE) plusk = eps->problem_type == EPS_GHIEP?0:1;
286:   if (plusk < 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid plusk value");
287:   data->minv = minv;
288:   data->plusk = plusk;
289:   return(0);
290: }

292: PetscErrorCode EPSXDGetRestart_XD(EPS eps,PetscInt *minv,PetscInt *plusk)
293: {
294:   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;

307:   *initialsize = data->initialsize;
308:   return(0);
309: }

311: PetscErrorCode EPSXDSetInitialSize_XD(EPS eps,PetscInt initialsize)
312: {
313:   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;

316:   if (initialsize == PETSC_DEFAULT || initialsize == PETSC_DECIDE) initialsize = 5;
317:   if (initialsize <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial size value");
318:   data->initialsize = initialsize;
319:   return(0);
320: }

322: PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
323: {
324:   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;

336:   *borth = data->ipB;
337:   return(0);
338: }

340: /*
341:   EPSComputeVectors_XD - Compute eigenvectors from the vectors
342:   provided by the eigensolver. This version is intended for solvers
343:   that provide Schur vectors from the QZ decomposition. Given the partial
344:   Schur decomposition OP*V=V*T, the following steps are performed:
345:       1) compute eigenvectors of (S,T): S*Z=T*Z*D
346:       2) compute eigenvectors of OP: X=V*Z
347:  */
348: PetscErrorCode EPSComputeVectors_XD(EPS eps)
349: {
351:   Mat            X;
352:   PetscBool      symm;

355:   PetscObjectTypeCompare((PetscObject)eps->ds,DSHEP,&symm);
356:   if (symm) return(0);
357:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);

359:   /* V <- V * X */
360:   DSGetMat(eps->ds,DS_MAT_X,&X);
361:   BVSetActiveColumns(eps->V,0,eps->nconv);
362:   BVMultInPlace(eps->V,X,0,eps->nconv);
363:   MatDestroy(&X);
364:   return(0);
365: }