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