Actual source code: dvdupdatev.c

slepc-3.21.1 2024-04-26
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:    SLEPc eigensolver: "davidson"

 13:    Step: test for restarting, updateV, restartV
 14: */

 16: #include "davidson.h"

 18: typedef struct {
 19:   PetscInt          min_size_V;        /* restart with this number of eigenvectors */
 20:   PetscInt          plusk;             /* at restart, save plusk vectors from last iteration */
 21:   PetscInt          mpd;               /* max size of the searching subspace */
 22:   void              *old_updateV_data; /* old updateV data */
 23:   PetscErrorCode    (*old_isRestarting)(dvdDashboard*,PetscBool*);  /* old isRestarting */
 24:   Mat               oldU;              /* previous projected right igenvectors */
 25:   Mat               oldV;              /* previous projected left eigenvectors */
 26:   PetscInt          size_oldU;         /* size of oldU */
 27:   PetscBool         allResiduals;      /* if computing all the residuals */
 28: } dvdManagV_basic;

 30: static PetscErrorCode dvd_updateV_start(dvdDashboard *d)
 31: {
 32:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
 33:   PetscInt        i;

 35:   PetscFunctionBegin;
 36:   for (i=0;i<d->eps->ncv;i++) d->eigi[i] = 0.0;
 37:   d->nR = d->real_nR;
 38:   for (i=0;i<d->eps->ncv;i++) d->nR[i] = 1.0;
 39:   d->nX = d->real_nX;
 40:   for (i=0;i<d->eps->ncv;i++) d->errest[i] = 1.0;
 41:   data->size_oldU = 0;
 42:   d->nconv = 0;
 43:   d->npreconv = 0;
 44:   d->V_tra_s = d->V_tra_e = d->V_new_s = d->V_new_e = 0;
 45:   d->size_D = 0;
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

 49: static PetscErrorCode dvd_isrestarting_fullV(dvdDashboard *d,PetscBool *r)
 50: {
 51:   PetscInt        l,k;
 52:   PetscBool       restart;
 53:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;

 55:   PetscFunctionBegin;
 56:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
 57:   restart = (k+2 > d->eps->ncv)? PETSC_TRUE: PETSC_FALSE;

 59:   /* Check old isRestarting function */
 60:   if (PetscUnlikely(!restart && data->old_isRestarting)) PetscCall(data->old_isRestarting(d,&restart));
 61:   *r = restart;
 62:   PetscFunctionReturn(PETSC_SUCCESS);
 63: }

 65: static PetscErrorCode dvd_managementV_basic_d(dvdDashboard *d)
 66: {
 67:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;

 69:   PetscFunctionBegin;
 70:   /* Restore changes in dvdDashboard */
 71:   d->updateV_data = data->old_updateV_data;

 73:   /* Free local data */
 74:   PetscCall(MatDestroy(&data->oldU));
 75:   PetscCall(MatDestroy(&data->oldV));
 76:   PetscCall(PetscFree(d->real_nR));
 77:   PetscCall(PetscFree(d->real_nX));
 78:   PetscCall(PetscFree(data));
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: static PetscErrorCode dvd_updateV_conv_gen(dvdDashboard *d)
 83: {
 84:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
 85:   PetscInt        npreconv,cMT,cMTX,lV,kV,nV;
 86:   Mat             Z,Z0,Q,Q0;
 87:   PetscBool       t;
 88: #if !defined(PETSC_USE_COMPLEX)
 89:   PetscInt        i;
 90: #endif

 92:   PetscFunctionBegin;
 93:   npreconv = d->npreconv;
 94:   /* Constrains the converged pairs to nev */
 95: #if !defined(PETSC_USE_COMPLEX)
 96:   /* Tries to maintain together conjugate eigenpairs */
 97:   for (i=0; (i + (d->eigi[i]!=0.0?1:0) < npreconv) && (d->nconv + i < d->nev); i+= (d->eigi[i]!=0.0?2:1));
 98:   npreconv = i;
 99: #else
100:   npreconv = PetscMax(PetscMin(d->nev-d->nconv,npreconv),0);
101: #endif
102:   /* For GHEP without B-ortho, converge all of the requested pairs at once */
103:   PetscCall(PetscObjectTypeCompare((PetscObject)d->eps->ds,DSGHEP,&t));
104:   if (t && d->nconv+npreconv<d->nev) npreconv = 0;
105:   /* Quick exit */
106:   if (npreconv == 0) PetscFunctionReturn(PETSC_SUCCESS);

108:   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
109:   nV  = kV - lV;
110:   cMT = nV - npreconv;
111:   /* Harmonics restarts with right eigenvectors, and other with the left ones.
112:      If the problem is standard or hermitian, left and right vectors are the same */
113:   if (!(d->W||DVD_IS(d->sEP,DVD_EP_STD)||DVD_IS(d->sEP,DVD_EP_HERMITIAN))) {
114:     /* ps.Q <- [ps.Q(0:npreconv-1) ps.Z(npreconv:size_H-1)] */
115:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Q,&Q));
116:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Z,&Z));
117:     PetscCall(MatDenseGetSubMatrix(Q,0,npreconv,nV,npreconv+cMT,&Q0));
118:     PetscCall(MatDenseGetSubMatrix(Z,0,npreconv,nV,npreconv+cMT,&Z0));
119:     PetscCall(MatCopy(Z0,Q0,SAME_NONZERO_PATTERN));
120:     PetscCall(MatDenseRestoreSubMatrix(Q,&Q0));
121:     PetscCall(MatDenseRestoreSubMatrix(Z,&Z0));
122:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q));
123:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z));
124:   }
125:   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) PetscCall(DSPseudoOrthogonalize(d->eps->ds,DS_MAT_Q,nV,d->nBds,&cMTX,d->nBds));
126:   else PetscCall(DSOrthogonalize(d->eps->ds,DS_MAT_Q,nV,&cMTX));
127:   cMT = cMTX - npreconv;

129:   if (d->W) {
130:     PetscCall(DSOrthogonalize(d->eps->ds,DS_MAT_Z,nV,&cMTX));
131:     cMT = PetscMin(cMT,cMTX - npreconv);
132:   }

134:   /* Lock the converged pairs */
135:   d->eigr+= npreconv;
136: #if !defined(PETSC_USE_COMPLEX)
137:   if (d->eigi) d->eigi+= npreconv;
138: #endif
139:   d->nconv+= npreconv;
140:   d->errest+= npreconv;
141:   /* Notify the changes in V and update the other subspaces */
142:   d->V_tra_s = npreconv;          d->V_tra_e = nV;
143:   d->V_new_s = cMT;               d->V_new_e = d->V_new_s;
144:   /* Remove oldU */
145:   data->size_oldU = 0;

147:   d->npreconv-= npreconv;
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: static PetscErrorCode dvd_updateV_restart_gen(dvdDashboard *d)
152: {
153:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
154:   PetscInt        lV,kV,nV,size_plusk,size_X,cMTX,cMTY,max_restart_size;
155:   Mat             Q,Q0,Z,Z0,U,V;

157:   PetscFunctionBegin;
158:   /* Select size_X desired pairs from V */
159:   /* The restarted basis should:
160:      - have at least one spot to add a new direction;
161:      - keep converged vectors, npreconv;
162:      - keep at least 1 oldU direction if possible.
163:   */
164:   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
165:   nV = kV - lV;
166:   max_restart_size = PetscMax(0,PetscMin(d->eps->mpd - 1,d->eps->ncv - lV - 2));
167:   size_X = PetscMin(PetscMin(data->min_size_V+d->npreconv,max_restart_size - (max_restart_size - d->npreconv > 1 && data->plusk > 0 && data->size_oldU > 0 ? 1 : 0)), nV);

169:   /* Add plusk eigenvectors from the previous iteration */
170:   size_plusk = PetscMax(0,PetscMin(PetscMin(PetscMin(data->plusk,data->size_oldU),max_restart_size - size_X),nV - size_X));

172:   d->size_MT = nV;
173:   /* ps.Q <- orth([pX(0:size_X-1) [oldU(0:size_plusk-1); 0] ]) */
174:   /* Harmonics restarts with right eigenvectors, and other with the left ones.
175:      If the problem is standard or hermitian, left and right vectors are the same */
176:   if (!(d->W||DVD_IS(d->sEP,DVD_EP_STD)||DVD_IS(d->sEP,DVD_EP_HERMITIAN))) {
177:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Q,&Q));
178:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Z,&Z));
179:     PetscCall(MatDenseGetSubMatrix(Q,0,nV,0,size_X,&Q0));
180:     PetscCall(MatDenseGetSubMatrix(Z,0,nV,0,size_X,&Z0));
181:     PetscCall(MatCopy(Z0,Q0,SAME_NONZERO_PATTERN));
182:     PetscCall(MatDenseRestoreSubMatrix(Q,&Q0));
183:     PetscCall(MatDenseRestoreSubMatrix(Z,&Z0));
184:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q));
185:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z));
186:   }
187:   PetscCheck(size_plusk<=0 || !DVD_IS(d->sEP,DVD_EP_INDEFINITE),PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported plusk>0 in indefinite eigenvalue problems");
188:   if (size_plusk > 0) {
189:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Q,&Q));
190:     PetscCall(MatDenseGetSubMatrix(Q,0,nV,size_X,size_X+size_plusk,&Q0));
191:     PetscCall(MatDenseGetSubMatrix(data->oldU,0,nV,0,size_plusk,&U));
192:     PetscCall(MatCopy(U,Q0,SAME_NONZERO_PATTERN));
193:     PetscCall(MatDenseRestoreSubMatrix(Q,&Q0));
194:     PetscCall(MatDenseRestoreSubMatrix(data->oldU,&U));
195:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q));
196:   }
197:   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) PetscCall(DSPseudoOrthogonalize(d->eps->ds,DS_MAT_Q,size_X,d->nBds,&cMTX,d->nBds));
198:   else PetscCall(DSOrthogonalize(d->eps->ds,DS_MAT_Q,size_X+size_plusk,&cMTX));

200:   if (d->W && size_plusk > 0) {
201:     /* ps.Z <- orth([ps.Z(0:size_X-1) [oldV(0:size_plusk-1); 0] ]) */
202:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Z,&Z));
203:     PetscCall(MatDenseGetSubMatrix(Z,0,nV,size_X,size_X+size_plusk,&Z0));
204:     PetscCall(MatDenseGetSubMatrix(data->oldV,0,nV,0,size_plusk,&V));
205:     PetscCall(MatCopy(V,Z0,SAME_NONZERO_PATTERN));
206:     PetscCall(MatDenseRestoreSubMatrix(Z,&Z0));
207:     PetscCall(MatDenseRestoreSubMatrix(data->oldV,&V));
208:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z));
209:     PetscCall(DSOrthogonalize(d->eps->ds,DS_MAT_Z,size_X+size_plusk,&cMTY));
210:     cMTX = PetscMin(cMTX, cMTY);
211:   }
212:   PetscAssert(cMTX<=size_X+size_plusk,PETSC_COMM_SELF,PETSC_ERR_SUP,"Invalid number of columns to restart");

214:   /* Notify the changes in V and update the other subspaces */
215:   d->V_tra_s = 0;                     d->V_tra_e = cMTX;
216:   d->V_new_s = d->V_tra_e;            d->V_new_e = d->V_new_s;

218:   /* Remove oldU */
219:   data->size_oldU = 0;

221:   /* Remove npreconv */
222:   d->npreconv = 0;
223:   PetscFunctionReturn(PETSC_SUCCESS);
224: }

226: static PetscErrorCode dvd_updateV_testConv(dvdDashboard *d,PetscInt s,PetscInt pre,PetscInt e,PetscInt *nConv)
227: {
228:   PetscInt        i,j,b;
229:   PetscReal       norm;
230:   PetscBool       conv, c;
231:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;

233:   PetscFunctionBegin;
234:   if (nConv) *nConv = s;
235:   for (i=s,conv=PETSC_TRUE;(conv || data->allResiduals) && (i < e);i+=b) {
236: #if !defined(PETSC_USE_COMPLEX)
237:     b = d->eigi[i]!=0.0?2:1;
238: #else
239:     b = 1;
240: #endif
241:     if (i+b-1 >= pre) PetscCall(d->calcpairs_residual(d,i,i+b));
242:     /* Test the Schur vector */
243:     for (j=0,c=PETSC_TRUE;j<b && c;j++) {
244:       norm = d->nR[i+j]/d->nX[i+j];
245:       c = d->testConv(d,d->eigr[i+j],d->eigi[i+j],norm,&d->errest[i+j]);
246:     }
247:     if (conv && c) { if (nConv) *nConv = i+b; }
248:     else conv = PETSC_FALSE;
249:   }
250:   pre = PetscMax(pre,i);

252: #if !defined(PETSC_USE_COMPLEX)
253:   /* Enforce converged conjugate complex eigenpairs */
254:   if (nConv) {
255:     for (j=0;j<*nConv;j++) if (d->eigi[j] != 0.0) j++;
256:     if (j>*nConv) (*nConv)--;
257:   }
258: #endif
259:   for (i=pre;i<e;i++) d->errest[i] = d->nR[i] = 1.0;
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: static PetscErrorCode dvd_updateV_update_gen(dvdDashboard *d)
264: {
265:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
266:   PetscInt        size_D,s,lV,kV,nV;
267:   Mat             Q,Q0,Z,Z0,U,V;

269:   PetscFunctionBegin;
270:   /* Select the desired pairs */
271:   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
272:   nV = kV - lV;
273:   size_D = PetscMin(PetscMin(PetscMin(d->bs,nV),d->eps->ncv-nV),nV);
274:   if (size_D == 0) PetscFunctionReturn(PETSC_SUCCESS);

276:   /* Fill V with D */
277:   PetscCall(d->improveX(d,d->npreconv,d->npreconv+size_D,&size_D));

279:   /* If D is empty, exit */
280:   d->size_D = size_D;
281:   if (size_D == 0) PetscFunctionReturn(PETSC_SUCCESS);

283:   /* Get the residual of all pairs */
284: #if !defined(PETSC_USE_COMPLEX)
285:   s = (d->eigi[0]!=0.0)? 2: 1;
286: #else
287:   s = 1;
288: #endif
289:   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
290:   nV = kV - lV;
291:   PetscCall(dvd_updateV_testConv(d,s,s,data->allResiduals?nV:size_D,NULL));

293:   /* Notify the changes in V */
294:   d->V_tra_s = 0;                 d->V_tra_e = 0;
295:   d->V_new_s = nV;                d->V_new_e = nV+size_D;

297:   /* Save the projected eigenvectors */
298:   if (data->plusk > 0) {
299:     PetscCall(MatZeroEntries(data->oldU));
300:     data->size_oldU = nV;
301:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_Q,&Q));
302:     PetscCall(MatDenseGetSubMatrix(Q,0,nV,0,nV,&Q0));
303:     PetscCall(MatDenseGetSubMatrix(data->oldU,0,nV,0,nV,&U));
304:     PetscCall(MatCopy(Q0,U,SAME_NONZERO_PATTERN));
305:     PetscCall(MatDenseRestoreSubMatrix(Q,&Q0));
306:     PetscCall(MatDenseRestoreSubMatrix(data->oldU,&U));
307:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q));
308:     if (d->W) {
309:       PetscCall(MatZeroEntries(data->oldV));
310:       PetscCall(DSGetMat(d->eps->ds,DS_MAT_Z,&Z));
311:       PetscCall(MatDenseGetSubMatrix(Z,0,nV,0,nV,&Z0));
312:       PetscCall(MatDenseGetSubMatrix(data->oldV,0,nV,0,nV,&V));
313:       PetscCall(MatCopy(Z0,V,SAME_NONZERO_PATTERN));
314:       PetscCall(MatDenseRestoreSubMatrix(Z,&Z0));
315:       PetscCall(MatDenseRestoreSubMatrix(data->oldV,&V));
316:       PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z));
317:     }
318:   }
319:   PetscFunctionReturn(PETSC_SUCCESS);
320: }

322: static PetscErrorCode dvd_updateV_extrapol(dvdDashboard *d)
323: {
324:   dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
325:   PetscInt        i;
326:   PetscBool       restart,t;

328:   PetscFunctionBegin;
329:   /* TODO: restrict select pairs to each case */
330:   PetscCall(d->calcpairs_selectPairs(d, data->min_size_V+d->npreconv));

332:   /* If the subspaces doesn't need restart, add new vector */
333:   PetscCall(d->isRestarting(d,&restart));
334:   if (!restart) {
335:     d->size_D = 0;
336:     PetscCall(dvd_updateV_update_gen(d));

338:     /* If no vector were converged, exit */
339:     /* For GHEP without B-ortho, converge all of the requested pairs at once */
340:     PetscCall(PetscObjectTypeCompare((PetscObject)d->eps->ds,DSGHEP,&t));
341:     if (d->nconv+d->npreconv < d->nev && (t || d->npreconv == 0)) PetscFunctionReturn(PETSC_SUCCESS);
342:   }

344:   /* If some eigenpairs were converged, lock them  */
345:   if (d->npreconv > 0) {
346:     i = d->npreconv;
347:     PetscCall(dvd_updateV_conv_gen(d));

349:     /* If some eigenpair was locked, exit */
350:     if (i > d->npreconv) PetscFunctionReturn(PETSC_SUCCESS);
351:   }

353:   /* Else, a restarting is performed */
354:   PetscCall(dvd_updateV_restart_gen(d));
355:   PetscFunctionReturn(PETSC_SUCCESS);
356: }

358: PetscErrorCode dvd_managementV_basic(dvdDashboard *d,dvdBlackboard *b,PetscInt bs,PetscInt mpd,PetscInt min_size_V,PetscInt plusk,PetscBool harm,PetscBool allResiduals)
359: {
360:   dvdManagV_basic *data;
361: #if !defined(PETSC_USE_COMPLEX)
362:   PetscBool       her_probl,std_probl;
363: #endif

365:   PetscFunctionBegin;
366:   /* Setting configuration constrains */
367: #if !defined(PETSC_USE_COMPLEX)
368:   /* if the last converged eigenvalue is complex its conjugate pair is also
369:      converged */
370:   her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
371:   std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
372:   b->max_size_X = PetscMax(b->max_size_X,bs+((her_probl && std_probl)?0:1));
373: #else
374:   b->max_size_X = PetscMax(b->max_size_X,bs);
375: #endif

377:   b->max_size_V = PetscMax(b->max_size_V,mpd);
378:   min_size_V = PetscMin(min_size_V,mpd-bs);
379:   b->size_V = PetscMax(b->size_V,b->max_size_V+b->max_size_P+b->max_nev);
380:   b->max_size_oldX = plusk;

382:   /* Setup the step */
383:   if (b->state >= DVD_STATE_CONF) {
384:     PetscCall(PetscNew(&data));
385:     data->mpd = b->max_size_V;
386:     data->min_size_V = min_size_V;
387:     d->bs = bs;
388:     data->plusk = plusk;
389:     data->allResiduals = allResiduals;

391:     d->eigr = d->eps->eigr;
392:     d->eigi = d->eps->eigi;
393:     d->errest = d->eps->errest;
394:     PetscCall(PetscMalloc1(d->eps->ncv,&d->real_nR));
395:     PetscCall(PetscMalloc1(d->eps->ncv,&d->real_nX));
396:     if (plusk > 0) PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&data->oldU));
397:     else data->oldU = NULL;
398:     if (harm && plusk>0) PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&data->oldV));
399:     else data->oldV = NULL;

401:     data->old_updateV_data = d->updateV_data;
402:     d->updateV_data = data;
403:     data->old_isRestarting = d->isRestarting;
404:     d->isRestarting = dvd_isrestarting_fullV;
405:     d->updateV = dvd_updateV_extrapol;
406:     d->preTestConv = dvd_updateV_testConv;
407:     PetscCall(EPSDavidsonFLAdd(&d->startList,dvd_updateV_start));
408:     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_managementV_basic_d));
409:   }
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }