Actual source code: dvdupdatev.c
slepc-3.18.3 2023-03-24
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: for (i=0;i<d->eps->ncv;i++) d->eigi[i] = 0.0;
36: d->nR = d->real_nR;
37: for (i=0;i<d->eps->ncv;i++) d->nR[i] = 1.0;
38: d->nX = d->real_nX;
39: for (i=0;i<d->eps->ncv;i++) d->errest[i] = 1.0;
40: data->size_oldU = 0;
41: d->nconv = 0;
42: d->npreconv = 0;
43: d->V_tra_s = d->V_tra_e = d->V_new_s = d->V_new_e = 0;
44: d->size_D = 0;
45: return 0;
46: }
48: static PetscErrorCode dvd_isrestarting_fullV(dvdDashboard *d,PetscBool *r)
49: {
50: PetscInt l,k;
51: PetscBool restart;
52: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
54: BVGetActiveColumns(d->eps->V,&l,&k);
55: restart = (k+2 > d->eps->ncv)? PETSC_TRUE: PETSC_FALSE;
57: /* Check old isRestarting function */
58: if (PetscUnlikely(!restart && data->old_isRestarting)) data->old_isRestarting(d,&restart);
59: *r = restart;
60: return 0;
61: }
63: static PetscErrorCode dvd_managementV_basic_d(dvdDashboard *d)
64: {
65: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
67: /* Restore changes in dvdDashboard */
68: d->updateV_data = data->old_updateV_data;
70: /* Free local data */
71: MatDestroy(&data->oldU);
72: MatDestroy(&data->oldV);
73: PetscFree(d->real_nR);
74: PetscFree(d->real_nX);
75: PetscFree(data);
76: return 0;
77: }
79: static PetscErrorCode dvd_updateV_conv_gen(dvdDashboard *d)
80: {
81: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
82: PetscInt npreconv,cMT,cMTX,lV,kV,nV;
83: Mat Z,Z0,Q,Q0;
84: PetscBool t;
85: #if !defined(PETSC_USE_COMPLEX)
86: PetscInt i;
87: #endif
89: npreconv = d->npreconv;
90: /* Constrains the converged pairs to nev */
91: #if !defined(PETSC_USE_COMPLEX)
92: /* Tries to maintain together conjugate eigenpairs */
93: 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));
94: npreconv = i;
95: #else
96: npreconv = PetscMax(PetscMin(d->nev-d->nconv,npreconv),0);
97: #endif
98: /* For GHEP without B-ortho, converge all of the requested pairs at once */
99: PetscObjectTypeCompare((PetscObject)d->eps->ds,DSGHEP,&t);
100: if (t && d->nconv+npreconv<d->nev) npreconv = 0;
101: /* Quick exit */
102: if (npreconv == 0) return 0;
104: BVGetActiveColumns(d->eps->V,&lV,&kV);
105: nV = kV - lV;
106: cMT = nV - npreconv;
107: /* Harmonics restarts with right eigenvectors, and other with the left ones.
108: If the problem is standard or hermitian, left and right vectors are the same */
109: if (!(d->W||DVD_IS(d->sEP,DVD_EP_STD)||DVD_IS(d->sEP,DVD_EP_HERMITIAN))) {
110: /* ps.Q <- [ps.Q(0:npreconv-1) ps.Z(npreconv:size_H-1)] */
111: DSGetMat(d->eps->ds,DS_MAT_Q,&Q);
112: DSGetMat(d->eps->ds,DS_MAT_Z,&Z);
113: MatDenseGetSubMatrix(Q,0,npreconv,nV,npreconv+cMT,&Q0);
114: MatDenseGetSubMatrix(Z,0,npreconv,nV,npreconv+cMT,&Z0);
115: MatCopy(Z0,Q0,SAME_NONZERO_PATTERN);
116: MatDenseRestoreSubMatrix(Q,&Q0);
117: MatDenseRestoreSubMatrix(Z,&Z0);
118: DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q);
119: DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z);
120: }
121: if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) DSPseudoOrthogonalize(d->eps->ds,DS_MAT_Q,nV,d->nBds,&cMTX,d->nBds);
122: else DSOrthogonalize(d->eps->ds,DS_MAT_Q,nV,&cMTX);
123: cMT = cMTX - npreconv;
125: if (d->W) {
126: DSOrthogonalize(d->eps->ds,DS_MAT_Z,nV,&cMTX);
127: cMT = PetscMin(cMT,cMTX - npreconv);
128: }
130: /* Lock the converged pairs */
131: d->eigr+= npreconv;
132: #if !defined(PETSC_USE_COMPLEX)
133: if (d->eigi) d->eigi+= npreconv;
134: #endif
135: d->nconv+= npreconv;
136: d->errest+= npreconv;
137: /* Notify the changes in V and update the other subspaces */
138: d->V_tra_s = npreconv; d->V_tra_e = nV;
139: d->V_new_s = cMT; d->V_new_e = d->V_new_s;
140: /* Remove oldU */
141: data->size_oldU = 0;
143: d->npreconv-= npreconv;
144: return 0;
145: }
147: static PetscErrorCode dvd_updateV_restart_gen(dvdDashboard *d)
148: {
149: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
150: PetscInt lV,kV,nV,size_plusk,size_X,cMTX,cMTY,max_restart_size;
151: Mat Q,Q0,Z,Z0,U,V;
153: /* Select size_X desired pairs from V */
154: /* The restarted basis should:
155: - have at least one spot to add a new direction;
156: - keep converged vectors, npreconv;
157: - keep at least 1 oldU direction if possible.
158: */
159: BVGetActiveColumns(d->eps->V,&lV,&kV);
160: nV = kV - lV;
161: max_restart_size = PetscMax(0,PetscMin(d->eps->mpd - 1,d->eps->ncv - lV - 2));
162: 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);
164: /* Add plusk eigenvectors from the previous iteration */
165: size_plusk = PetscMax(0,PetscMin(PetscMin(PetscMin(data->plusk,data->size_oldU),max_restart_size - size_X),nV - size_X));
167: d->size_MT = nV;
168: /* ps.Q <- orth([pX(0:size_X-1) [oldU(0:size_plusk-1); 0] ]) */
169: /* Harmonics restarts with right eigenvectors, and other with the left ones.
170: If the problem is standard or hermitian, left and right vectors are the same */
171: if (!(d->W||DVD_IS(d->sEP,DVD_EP_STD)||DVD_IS(d->sEP,DVD_EP_HERMITIAN))) {
172: DSGetMat(d->eps->ds,DS_MAT_Q,&Q);
173: DSGetMat(d->eps->ds,DS_MAT_Z,&Z);
174: MatDenseGetSubMatrix(Q,0,nV,0,size_X,&Q0);
175: MatDenseGetSubMatrix(Z,0,nV,0,size_X,&Z0);
176: MatCopy(Z0,Q0,SAME_NONZERO_PATTERN);
177: MatDenseRestoreSubMatrix(Q,&Q0);
178: MatDenseRestoreSubMatrix(Z,&Z0);
179: DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q);
180: DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z);
181: }
183: if (size_plusk > 0) {
184: DSGetMat(d->eps->ds,DS_MAT_Q,&Q);
185: MatDenseGetSubMatrix(Q,0,nV,size_X,size_X+size_plusk,&Q0);
186: MatDenseGetSubMatrix(data->oldU,0,nV,0,size_plusk,&U);
187: MatCopy(U,Q0,SAME_NONZERO_PATTERN);
188: MatDenseRestoreSubMatrix(Q,&Q0);
189: MatDenseRestoreSubMatrix(data->oldU,&U);
190: DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q);
191: }
192: if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) DSPseudoOrthogonalize(d->eps->ds,DS_MAT_Q,size_X,d->nBds,&cMTX,d->nBds);
193: else DSOrthogonalize(d->eps->ds,DS_MAT_Q,size_X+size_plusk,&cMTX);
195: if (d->W && size_plusk > 0) {
196: /* ps.Z <- orth([ps.Z(0:size_X-1) [oldV(0:size_plusk-1); 0] ]) */
197: DSGetMat(d->eps->ds,DS_MAT_Z,&Z);
198: MatDenseGetSubMatrix(Z,0,nV,size_X,size_X+size_plusk,&Z0);
199: MatDenseGetSubMatrix(data->oldV,0,nV,0,size_plusk,&V);
200: MatCopy(V,Z0,SAME_NONZERO_PATTERN);
201: MatDenseRestoreSubMatrix(Z,&Z0);
202: MatDenseRestoreSubMatrix(data->oldV,&V);
203: DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z);
204: DSOrthogonalize(d->eps->ds,DS_MAT_Z,size_X+size_plusk,&cMTY);
205: cMTX = PetscMin(cMTX, cMTY);
206: }
207: PetscAssert(cMTX<=size_X+size_plusk,PETSC_COMM_SELF,PETSC_ERR_SUP,"Invalid number of columns to restart");
209: /* Notify the changes in V and update the other subspaces */
210: d->V_tra_s = 0; d->V_tra_e = cMTX;
211: d->V_new_s = d->V_tra_e; d->V_new_e = d->V_new_s;
213: /* Remove oldU */
214: data->size_oldU = 0;
216: /* Remove npreconv */
217: d->npreconv = 0;
218: return 0;
219: }
221: static PetscErrorCode dvd_updateV_testConv(dvdDashboard *d,PetscInt s,PetscInt pre,PetscInt e,PetscInt *nConv)
222: {
223: PetscInt i,j,b;
224: PetscReal norm;
225: PetscBool conv, c;
226: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
228: if (nConv) *nConv = s;
229: for (i=s,conv=PETSC_TRUE;(conv || data->allResiduals) && (i < e);i+=b) {
230: #if !defined(PETSC_USE_COMPLEX)
231: b = d->eigi[i]!=0.0?2:1;
232: #else
233: b = 1;
234: #endif
235: if (i+b-1 >= pre) d->calcpairs_residual(d,i,i+b);
236: /* Test the Schur vector */
237: for (j=0,c=PETSC_TRUE;j<b && c;j++) {
238: norm = d->nR[i+j]/d->nX[i+j];
239: c = d->testConv(d,d->eigr[i+j],d->eigi[i+j],norm,&d->errest[i+j]);
240: }
241: if (conv && c) { if (nConv) *nConv = i+b; }
242: else conv = PETSC_FALSE;
243: }
244: pre = PetscMax(pre,i);
246: #if !defined(PETSC_USE_COMPLEX)
247: /* Enforce converged conjugate complex eigenpairs */
248: if (nConv) {
249: for (j=0;j<*nConv;j++) if (d->eigi[j] != 0.0) j++;
250: if (j>*nConv) (*nConv)--;
251: }
252: #endif
253: for (i=pre;i<e;i++) d->errest[i] = d->nR[i] = 1.0;
254: return 0;
255: }
257: static PetscErrorCode dvd_updateV_update_gen(dvdDashboard *d)
258: {
259: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
260: PetscInt size_D,s,lV,kV,nV;
261: Mat Q,Q0,Z,Z0,U,V;
263: /* Select the desired pairs */
264: BVGetActiveColumns(d->eps->V,&lV,&kV);
265: nV = kV - lV;
266: size_D = PetscMin(PetscMin(PetscMin(d->bs,nV),d->eps->ncv-nV),nV);
267: if (size_D == 0) return 0;
269: /* Fill V with D */
270: d->improveX(d,d->npreconv,d->npreconv+size_D,&size_D);
272: /* If D is empty, exit */
273: d->size_D = size_D;
274: if (size_D == 0) return 0;
276: /* Get the residual of all pairs */
277: #if !defined(PETSC_USE_COMPLEX)
278: s = (d->eigi[0]!=0.0)? 2: 1;
279: #else
280: s = 1;
281: #endif
282: BVGetActiveColumns(d->eps->V,&lV,&kV);
283: nV = kV - lV;
284: dvd_updateV_testConv(d,s,s,data->allResiduals?nV:size_D,NULL);
286: /* Notify the changes in V */
287: d->V_tra_s = 0; d->V_tra_e = 0;
288: d->V_new_s = nV; d->V_new_e = nV+size_D;
290: /* Save the projected eigenvectors */
291: if (data->plusk > 0) {
292: MatZeroEntries(data->oldU);
293: data->size_oldU = nV;
294: DSGetMat(d->eps->ds,DS_MAT_Q,&Q);
295: MatDenseGetSubMatrix(Q,0,nV,0,nV,&Q0);
296: MatDenseGetSubMatrix(data->oldU,0,nV,0,nV,&U);
297: MatCopy(Q0,U,SAME_NONZERO_PATTERN);
298: MatDenseRestoreSubMatrix(Q,&Q0);
299: MatDenseRestoreSubMatrix(data->oldU,&U);
300: DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q);
301: if (d->W) {
302: MatZeroEntries(data->oldV);
303: DSGetMat(d->eps->ds,DS_MAT_Z,&Z);
304: MatDenseGetSubMatrix(Z,0,nV,0,nV,&Z0);
305: MatDenseGetSubMatrix(data->oldV,0,nV,0,nV,&V);
306: MatCopy(Z0,V,SAME_NONZERO_PATTERN);
307: MatDenseRestoreSubMatrix(Z,&Z0);
308: MatDenseRestoreSubMatrix(data->oldV,&V);
309: DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z);
310: }
311: }
312: return 0;
313: }
315: static PetscErrorCode dvd_updateV_extrapol(dvdDashboard *d)
316: {
317: dvdManagV_basic *data = (dvdManagV_basic*)d->updateV_data;
318: PetscInt i;
319: PetscBool restart,t;
321: /* TODO: restrict select pairs to each case */
322: d->calcpairs_selectPairs(d, data->min_size_V+d->npreconv);
324: /* If the subspaces doesn't need restart, add new vector */
325: d->isRestarting(d,&restart);
326: if (!restart) {
327: d->size_D = 0;
328: dvd_updateV_update_gen(d);
330: /* If no vector were converged, exit */
331: /* For GHEP without B-ortho, converge all of the requested pairs at once */
332: PetscObjectTypeCompare((PetscObject)d->eps->ds,DSGHEP,&t);
333: if (d->nconv+d->npreconv < d->nev && (t || d->npreconv == 0)) return 0;
334: }
336: /* If some eigenpairs were converged, lock them */
337: if (d->npreconv > 0) {
338: i = d->npreconv;
339: dvd_updateV_conv_gen(d);
341: /* If some eigenpair was locked, exit */
342: if (i > d->npreconv) return 0;
343: }
345: /* Else, a restarting is performed */
346: dvd_updateV_restart_gen(d);
347: return 0;
348: }
350: PetscErrorCode dvd_managementV_basic(dvdDashboard *d,dvdBlackboard *b,PetscInt bs,PetscInt mpd,PetscInt min_size_V,PetscInt plusk,PetscBool harm,PetscBool allResiduals)
351: {
352: dvdManagV_basic *data;
353: #if !defined(PETSC_USE_COMPLEX)
354: PetscBool her_probl,std_probl;
355: #endif
357: /* Setting configuration constrains */
358: #if !defined(PETSC_USE_COMPLEX)
359: /* if the last converged eigenvalue is complex its conjugate pair is also
360: converged */
361: her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
362: std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
363: b->max_size_X = PetscMax(b->max_size_X,bs+((her_probl && std_probl)?0:1));
364: #else
365: b->max_size_X = PetscMax(b->max_size_X,bs);
366: #endif
368: b->max_size_V = PetscMax(b->max_size_V,mpd);
369: min_size_V = PetscMin(min_size_V,mpd-bs);
370: b->size_V = PetscMax(b->size_V,b->max_size_V+b->max_size_P+b->max_nev);
371: b->max_size_oldX = plusk;
373: /* Setup the step */
374: if (b->state >= DVD_STATE_CONF) {
375: PetscNew(&data);
376: data->mpd = b->max_size_V;
377: data->min_size_V = min_size_V;
378: d->bs = bs;
379: data->plusk = plusk;
380: data->allResiduals = allResiduals;
382: d->eigr = d->eps->eigr;
383: d->eigi = d->eps->eigi;
384: d->errest = d->eps->errest;
385: PetscMalloc1(d->eps->ncv,&d->real_nR);
386: PetscMalloc1(d->eps->ncv,&d->real_nX);
387: if (plusk > 0) MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&data->oldU);
388: else data->oldU = NULL;
389: if (harm && plusk>0) MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&data->oldV);
390: else data->oldV = NULL;
392: data->old_updateV_data = d->updateV_data;
393: d->updateV_data = data;
394: data->old_isRestarting = d->isRestarting;
395: d->isRestarting = dvd_isrestarting_fullV;
396: d->updateV = dvd_updateV_extrapol;
397: d->preTestConv = dvd_updateV_testConv;
398: EPSDavidsonFLAdd(&d->startList,dvd_updateV_start);
399: EPSDavidsonFLAdd(&d->destroyList,dvd_managementV_basic_d);
400: }
401: return 0;
402: }