Actual source code: dvdcalcpairs.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:    SLEPc eigensolver: "davidson"

 13:    Step: calculate the best eigenpairs in the subspace V

 15:    For that, performs these steps:
 16:      1) Update W <- A * V
 17:      2) Update H <- V' * W
 18:      3) Obtain eigenpairs of H
 19:      4) Select some eigenpairs
 20:      5) Compute the Ritz pairs of the selected ones
 21: */

 23: #include "davidson.h"
 24: #include <slepcblaslapack.h>

 26: static PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
 27: {

 31:   BVSetActiveColumns(d->eps->V,0,0);
 32:   if (d->W) { BVSetActiveColumns(d->W,0,0); }
 33:   BVSetActiveColumns(d->AX,0,0);
 34:   if (d->BX) { BVSetActiveColumns(d->BX,0,0); }
 35:   return(0);
 36: }

 38: static PetscErrorCode dvd_calcpairs_qz_d(dvdDashboard *d)
 39: {
 40:   PetscErrorCode  ierr;

 43:   BVDestroy(&d->W);
 44:   BVDestroy(&d->AX);
 45:   BVDestroy(&d->BX);
 46:   BVDestroy(&d->auxBV);
 47:   MatDestroy(&d->H);
 48:   MatDestroy(&d->G);
 49:   MatDestroy(&d->auxM);
 50:   SlepcVecPoolDestroy(&d->auxV);
 51:   PetscFree(d->nBds);
 52:   return(0);
 53: }

 55: /* in complex, d->size_H real auxiliary values are needed */
 56: static PetscErrorCode dvd_calcpairs_projeig_solve(dvdDashboard *d)
 57: {
 58:   PetscErrorCode    ierr;
 59:   Vec               v;
 60:   PetscScalar       *pA;
 61:   const PetscScalar *pv;
 62:   PetscInt          i,lV,kV,n,ld;

 65:   BVGetActiveColumns(d->eps->V,&lV,&kV);
 66:   n = kV-lV;
 67:   DSSetDimensions(d->eps->ds,n,0,0,0);
 68:   DSCopyMat(d->eps->ds,DS_MAT_A,0,0,d->H,lV,lV,n,n,PETSC_FALSE);
 69:   if (d->G) {
 70:     DSCopyMat(d->eps->ds,DS_MAT_B,0,0,d->G,lV,lV,n,n,PETSC_FALSE);
 71:   }
 72:   /* Set the signature on projected matrix B */
 73:   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
 74:     DSGetLeadingDimension(d->eps->ds,&ld);
 75:     DSGetArray(d->eps->ds,DS_MAT_B,&pA);
 76:     PetscMemzero(pA,sizeof(PetscScalar)*n*ld);
 77:     VecCreateSeq(PETSC_COMM_SELF,kV,&v);
 78:     BVGetSignature(d->eps->V,v);
 79:     VecGetArrayRead(v,&pv);
 80:     for (i=0;i<n;i++) {
 81:       pA[i+ld*i] = d->nBds[i] = PetscRealPart(pv[lV+i]);
 82:     }
 83:     VecRestoreArrayRead(v,&pv);
 84:     VecDestroy(&v);
 85:     DSRestoreArray(d->eps->ds,DS_MAT_B,&pA);
 86:   }
 87:   DSSetState(d->eps->ds,DS_STATE_RAW);
 88:   DSSolve(d->eps->ds,d->eigr,d->eigi);
 89:   return(0);
 90: }

 92: /*
 93:    A(lA:kA-1,lA:kA-1) <- Z(l:k-1)'*A(l:k-1,l:k-1)*Q(l,k-1), where k=l+kA-lA
 94:  */
 95: static PetscErrorCode EPSXDUpdateProj(Mat Q,Mat Z,PetscInt l,Mat A,PetscInt lA,PetscInt kA,Mat aux)
 96: {
 98:   PetscScalar    one=1.0,zero=0.0;
 99:   PetscInt       i,j,dA_=kA-lA,m0,n0,ldA_,ldQ_,ldZ_,nQ_;
100:   PetscBLASInt   dA,nQ,ldA,ldQ,ldZ;
101:   PetscScalar    *pA,*pQ,*pZ,*pW;
102:   PetscBool      symm=PETSC_FALSE,set,flg;

105:   MatGetSize(A,&m0,&n0); ldA_=m0;
106:   if (m0!=n0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"A should be square");
107:   if (lA<0 || lA>m0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial row, column in A");
108:   if (kA<0 || kA<lA || kA>m0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid final row, column in A");
109:   MatIsHermitianKnown(A,&set,&flg);
110:   symm = set? flg: PETSC_FALSE;
111:   MatGetSize(Q,&m0,&n0); ldQ_=nQ_=m0;
112:   if (l<0 || l>n0 || l+dA_>n0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Q");
113:   MatGetSize(Z,&m0,&n0); ldZ_=m0;
114:   if (l<0 || l>n0 || l+dA_>n0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Z");
115:   MatGetSize(aux,&m0,&n0);
116:   if (m0*n0<nQ_*dA_) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"aux should be larger");
117:   PetscBLASIntCast(dA_,&dA);
118:   PetscBLASIntCast(nQ_,&nQ);
119:   PetscBLASIntCast(ldA_,&ldA);
120:   PetscBLASIntCast(ldQ_,&ldQ);
121:   PetscBLASIntCast(ldZ_,&ldZ);
122:   MatDenseGetArray(A,&pA);
123:   MatDenseGetArray(Q,&pQ);
124:   if (Q!=Z) { MatDenseGetArray(Z,&pZ); }
125:   else pZ = pQ;
126: #if PETSC_USE_DEBUG
127:   /* Avoid valgrind warning in xgemm and xsymm */
128:   MatZeroEntries(aux);
129: #endif
130:   MatDenseGetArray(aux,&pW);
131:   /* W = A*Q */
132:   if (symm) {
133:     /* symmetrize before multiplying */
134:     for (i=lA+1;i<lA+nQ;i++) {
135:       for (j=lA;j<i;j++) pA[i+j*ldA] = PetscConj(pA[j+i*ldA]);
136:     }
137:   }
138:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&nQ,&dA,&nQ,&one,&pA[ldA*lA+lA],&ldA,&pQ[ldQ*l+l],&ldQ,&zero,pW,&nQ));
139:   /* A = Q'*W */
140:   PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&dA,&dA,&nQ,&one,&pZ[ldZ*l+l],&ldZ,pW,&nQ,&zero,&pA[ldA*lA+lA],&ldA));
141:   MatDenseRestoreArray(A,&pA);
142:   MatDenseRestoreArray(Q,&pQ);
143:   if (Q!=Z) { MatDenseGetArray(Z,&pZ); }
144:   else pZ = pQ;
145:   MatDenseRestoreArray(aux,&pW);
146:   return(0);
147: }

149: static PetscErrorCode dvd_calcpairs_updateproj(dvdDashboard *d)
150: {
152:   Mat            Q,Z;
153:   PetscInt       lV,kV;
154:   PetscBool      symm;

157:   DSGetMat(d->eps->ds,DS_MAT_Q,&Q);
158:   if (d->W) { DSGetMat(d->eps->ds,DS_MAT_Z,&Z); }
159:   else Z = Q;
160:   BVGetActiveColumns(d->eps->V,&lV,&kV);
161:   EPSXDUpdateProj(Q,Z,0,d->H,lV,lV+d->V_tra_e,d->auxM);
162:   if (d->G) { EPSXDUpdateProj(Q,Z,0,d->G,lV,lV+d->V_tra_e,d->auxM); }
163:   MatDestroy(&Q);
164:   if (d->W) { MatDestroy(&Z); }

166:   PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSHEP,DSGHIEP,DSGHEP,"");
167:   if (d->V_tra_s==0 || symm) return(0);
168:   /* Compute upper part of H (and G): H(0:l-1,l:k-1) <- W(0:l-1)' * AV(l:k-1), where
169:      k=l+d->V_tra_s */
170:   BVSetActiveColumns(d->W?d->W:d->eps->V,0,lV);
171:   BVSetActiveColumns(d->AX,lV,lV+d->V_tra_s);
172:   BVDot(d->AX,d->W?d->W:d->eps->V,d->H);
173:   if (d->G) {
174:     BVSetActiveColumns(d->BX?d->BX:d->eps->V,lV,lV+d->V_tra_s);
175:     BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G);
176:   }
177:   PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSGHEP,"");
178:   if (!symm) {
179:     BVSetActiveColumns(d->W?d->W:d->eps->V,lV,lV+d->V_tra_s);
180:     BVSetActiveColumns(d->AX,0,lV);
181:     BVDot(d->AX,d->W?d->W:d->eps->V,d->H);
182:     if (d->G) {
183:       BVSetActiveColumns(d->BX?d->BX:d->eps->V,0,lV);
184:       BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G);
185:     }
186:   }
187:   BVSetActiveColumns(d->eps->V,lV,kV);
188:   BVSetActiveColumns(d->AX,lV,kV);
189:   if (d->BX) { BVSetActiveColumns(d->BX,lV,kV); }
190:   if (d->W) { BVSetActiveColumns(d->W,lV,kV); }
191:   if (d->W) { dvd_harm_updateproj(d); }
192:   return(0);
193: }

195: /*
196:    BV <- BV*MT
197:  */
198: PETSC_STATIC_INLINE PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,BV bv,DSMatType mat)
199: {
201:   PetscInt       l,k,n;
202:   Mat            auxM;

205:   BVGetActiveColumns(d->eps->V,&l,&k);
206:   MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&auxM);
207:   MatZeroEntries(auxM);
208:   DSGetDimensions(d->eps->ds,&n,NULL,NULL,NULL,NULL);
209:   if (k-l!=n) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
210:   DSCopyMat(d->eps->ds,mat,0,0,auxM,l,l,n,d->V_tra_e,PETSC_TRUE);
211:   BVMultInPlace(bv,auxM,l,l+d->V_tra_e);
212:   MatDestroy(&auxM);
213:   return(0);
214: }

216: static PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
217: {
219:   PetscInt       i,l,k;
220:   Vec            v1,v2;
221:   PetscScalar    *pv;

224:   BVGetActiveColumns(d->eps->V,&l,&k);
225:   /* Update AV, BV, W and the projected matrices */
226:   /* 1. S <- S*MT */
227:   if (d->V_tra_s != d->V_tra_e || d->V_tra_e > 0) {
228:     dvd_calcpairs_updateBV0_gen(d,d->eps->V,DS_MAT_Q);
229:     if (d->W) { dvd_calcpairs_updateBV0_gen(d,d->W,DS_MAT_Z); }
230:     dvd_calcpairs_updateBV0_gen(d,d->AX,DS_MAT_Q);
231:     if (d->BX) { dvd_calcpairs_updateBV0_gen(d,d->BX,DS_MAT_Q); }
232:     dvd_calcpairs_updateproj(d);
233:     /* Update signature */
234:     if (d->nBds) {
235:       VecCreateSeq(PETSC_COMM_SELF,l+d->V_tra_e,&v1);
236:       BVSetActiveColumns(d->eps->V,0,l+d->V_tra_e);
237:       BVGetSignature(d->eps->V,v1);
238:       VecGetArray(v1,&pv);
239:       for (i=0;i<d->V_tra_e;i++) pv[l+i] = d->nBds[i];
240:       VecRestoreArray(v1,&pv);
241:       BVSetSignature(d->eps->V,v1);
242:       BVSetActiveColumns(d->eps->V,l,k);
243:       VecDestroy(&v1);
244:     }
245:     k = l+d->V_tra_e;
246:     l+= d->V_tra_s;
247:   } else {
248:     /* 2. V <- orth(V, V_new) */
249:     dvd_orthV(d->eps->V,l+d->V_new_s,l+d->V_new_e);
250:     /* 3. AV <- [AV A * V(V_new_s:V_new_e-1)] */
251:     /* Check consistency */
252:     if (k-l != d->V_new_s) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
253:     for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
254:       BVGetColumn(d->eps->V,i,&v1);
255:       BVGetColumn(d->AX,i,&v2);
256:       MatMult(d->A,v1,v2);
257:       BVRestoreColumn(d->eps->V,i,&v1);
258:       BVRestoreColumn(d->AX,i,&v2);
259:     }
260:     /* 4. BV <- [BV B * V(V_new_s:V_new_e-1)] */
261:     if (d->BX) {
262:       /* Check consistency */
263:       if (k-l != d->V_new_s) SETERRQ(PETSC_COMM_SELF,1,"Consistency broken");
264:       for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
265:         BVGetColumn(d->eps->V,i,&v1);
266:         BVGetColumn(d->BX,i,&v2);
267:         MatMult(d->B,v1,v2);
268:         BVRestoreColumn(d->eps->V,i,&v1);
269:         BVRestoreColumn(d->BX,i,&v2);
270:       }
271:     }
272:     /* 5. W <- [W f(AV,BV)] */
273:     if (d->W) {
274:       d->calcpairs_W(d);
275:       dvd_orthV(d->W,l+d->V_new_s,l+d->V_new_e);
276:     }
277:     /* 6. H <- W' * AX; G <- W' * BX */
278:     BVSetActiveColumns(d->eps->V,l+d->V_new_s,l+d->V_new_e);
279:     BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e);
280:     if (d->BX) { BVSetActiveColumns(d->BX,l+d->V_new_s,l+d->V_new_e); }
281:     if (d->W) { BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e); }
282:     BVMatProject(d->AX,NULL,d->W?d->W:d->eps->V,d->H);
283:     if (d->G) { BVMatProject(d->BX?d->BX:d->eps->V,NULL,d->W?d->W:d->eps->V,d->G); }
284:     BVSetActiveColumns(d->eps->V,l,k);
285:     BVSetActiveColumns(d->AX,l,k);
286:     if (d->BX) { BVSetActiveColumns(d->BX,l,k); }
287:     if (d->W) { BVSetActiveColumns(d->W,l,k); }

289:     /* Perform the transformation on the projected problem */
290:     if (d->W) {
291:       d->calcpairs_proj_trans(d);
292:     }
293:     k = l+d->V_new_e;
294:   }
295:   BVSetActiveColumns(d->eps->V,l,k);
296:   BVSetActiveColumns(d->AX,l,k);
297:   if (d->BX) { BVSetActiveColumns(d->BX,l,k); }
298:   if (d->W) { BVSetActiveColumns(d->W,l,k); }

300:   /* Solve the projected problem */
301:   dvd_calcpairs_projeig_solve(d);

303:   d->V_tra_s = d->V_tra_e = 0;
304:   d->V_new_s = d->V_new_e;
305:   return(0);
306: }

308: static PetscErrorCode dvd_calcpairs_apply_arbitrary(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscScalar *rr,PetscScalar *ri)
309: {
310:   PetscInt       i,k,ld;
311:   PetscScalar    *pX;
312:   Vec            *X,xr,xi;
314: #if defined(PETSC_USE_COMPLEX)
315:   PetscInt       N=1;
316: #else
317:   PetscInt       N=2,j;
318: #endif

321:   /* Quick exit without neither arbitrary selection nor harmonic extraction */
322:   if (!d->eps->arbitrary && !d->calcpairs_eig_backtrans) return(0);

324:   /* Quick exit without arbitrary selection, but with harmonic extraction */
325:   if (d->calcpairs_eig_backtrans) {
326:     for (i=r_s; i<r_e; i++) {
327:       d->calcpairs_eig_backtrans(d,d->eigr[i],d->eigi[i],&rr[i-r_s],&ri[i-r_s]);
328:     }
329:   }
330:   if (!d->eps->arbitrary) return(0);

332:   SlepcVecPoolGetVecs(d->auxV,N,&X);
333:   DSGetLeadingDimension(d->eps->ds,&ld);
334:   for (i=r_s;i<r_e;i++) {
335:     k = i;
336:     DSVectors(d->eps->ds,DS_MAT_X,&k,NULL);
337:     DSGetArray(d->eps->ds,DS_MAT_X,&pX);
338:     dvd_improvex_compute_X(d,i,k+1,X,pX,ld);
339:     DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);
340: #if !defined(PETSC_USE_COMPLEX)
341:     if (d->nX[i] != 1.0) {
342:       for (j=i;j<k+1;j++) {
343:         VecScale(X[j-i],1.0/d->nX[i]);
344:       }
345:     }
346:     xr = X[0];
347:     xi = X[1];
348:     if (i == k) {
349:       VecSet(xi,0.0);
350:     }
351: #else
352:     xr = X[0];
353:     xi = NULL;
354:     if (d->nX[i] != 1.0) {
355:       VecScale(xr,1.0/d->nX[i]);
356:     }
357: #endif
358:     (d->eps->arbitrary)(rr[i-r_s],ri[i-r_s],xr,xi,&rr[i-r_s],&ri[i-r_s],d->eps->arbitraryctx);
359: #if !defined(PETSC_USE_COMPLEX)
360:     if (i != k) {
361:       rr[i+1-r_s] = rr[i-r_s];
362:       ri[i+1-r_s] = ri[i-r_s];
363:       i++;
364:     }
365: #endif
366:   }
367:   SlepcVecPoolRestoreVecs(d->auxV,N,&X);
368:   return(0);
369: }

371: static PetscErrorCode dvd_calcpairs_selectPairs(dvdDashboard *d,PetscInt n)
372: {
373:   PetscInt       k,lV,kV,nV;
374:   PetscScalar    *rr,*ri;

378:   BVGetActiveColumns(d->eps->V,&lV,&kV);
379:   nV = kV - lV;
380:   n = PetscMin(n,nV);
381:   if (n <= 0) return(0);
382:   /* Put the best n pairs at the beginning. Useful for restarting */
383:   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
384:     PetscMalloc1(nV,&rr);
385:     PetscMalloc1(nV,&ri);
386:     dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri);
387:   } else {
388:     rr = d->eigr;
389:     ri = d->eigi;
390:   }
391:   k = n;
392:   DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k);
393:   /* Put the best pair at the beginning. Useful to check its residual */
394: #if !defined(PETSC_USE_COMPLEX)
395:   if (n != 1 && (n != 2 || d->eigi[0] == 0.0))
396: #else
397:   if (n != 1)
398: #endif
399:   {
400:     dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri);
401:     k = 1;
402:     DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k);
403:   }
404:   DSSynchronize(d->eps->ds,d->eigr,d->eigi);

406:   if (d->calcpairs_eigs_trans) {
407:     d->calcpairs_eigs_trans(d);
408:   }
409:   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
410:     PetscFree(rr);
411:     PetscFree(ri);
412:   }
413:   return(0);
414: }

416: static PetscErrorCode EPSXDComputeDSConv(dvdDashboard *d)
417: {
418:   PetscErrorCode    ierr;
419:   PetscInt          i,ld;
420:   Vec               v;
421:   PetscScalar       *pA;
422:   const PetscScalar *pv;
423:   PetscBool         symm;

426:   BVSetActiveColumns(d->eps->V,0,d->eps->nconv);
427:   PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSHEP,"");
428:   if (symm) return(0);
429:   DSSetDimensions(d->eps->ds,d->eps->nconv,0,0,0);
430:   DSCopyMat(d->eps->ds,DS_MAT_A,0,0,d->H,0,0,d->eps->nconv,d->eps->nconv,PETSC_FALSE);
431:   if (d->G) {
432:     DSCopyMat(d->eps->ds,DS_MAT_B,0,0,d->G,0,0,d->eps->nconv,d->eps->nconv,PETSC_FALSE);
433:   }
434:   /* Set the signature on projected matrix B */
435:   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
436:     DSGetLeadingDimension(d->eps->ds,&ld);
437:     DSGetArray(d->eps->ds,DS_MAT_B,&pA);
438:     PetscMemzero(pA,sizeof(PetscScalar)*d->eps->nconv*ld);
439:     VecCreateSeq(PETSC_COMM_SELF,d->eps->nconv,&v);
440:     BVGetSignature(d->eps->V,v);
441:     VecGetArrayRead(v,&pv);
442:     for (i=0;i<d->eps->nconv;i++) pA[i+ld*i] = pv[i];
443:     VecRestoreArrayRead(v,&pv);
444:     VecDestroy(&v);
445:     DSRestoreArray(d->eps->ds,DS_MAT_B,&pA);
446:   }
447:   DSSetState(d->eps->ds,DS_STATE_RAW);
448:   DSSolve(d->eps->ds,d->eps->eigr,d->eps->eigi);
449:   DSSynchronize(d->eps->ds,d->eps->eigr,d->eps->eigi);
450:   if (d->W) {
451:     for (i=0; i<d->eps->nconv; i++) {
452:       d->calcpairs_eig_backtrans(d,d->eps->eigr[i],d->eps->eigi[i],&d->eps->eigr[i],&d->eps->eigi[i]);
453:     }
454:   }
455:   return(0);
456: }

458: /*
459:    Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
460:    the norm associated to the Schur pair, where i = r_s..r_e
461: */
462: static PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d,PetscInt r_s,PetscInt r_e)
463: {
464:   PetscInt       i,ldpX;
465:   PetscScalar    *pX;
467:   BV             BX = d->BX?d->BX:d->eps->V;
468:   Vec            *R;

471:   DSGetLeadingDimension(d->eps->ds,&ldpX);
472:   DSGetArray(d->eps->ds,DS_MAT_Q,&pX);
473:   /* nX(i) <- ||X(i)|| */
474:   dvd_improvex_compute_X(d,r_s,r_e,NULL,pX,ldpX);
475:   SlepcVecPoolGetVecs(d->auxV,r_e-r_s,&R);
476:   for (i=r_s;i<r_e;i++) {
477:     /* R(i-r_s) <- AV*pX(i) */
478:     BVMultVec(d->AX,1.0,0.0,R[i-r_s],&pX[ldpX*i]);
479:     /* R(i-r_s) <- R(i-r_s) - eigr(i)*BX*pX(i) */
480:     BVMultVec(BX,-d->eigr[i],1.0,R[i-r_s],&pX[ldpX*i]);
481:   }
482:   DSRestoreArray(d->eps->ds,DS_MAT_Q,&pX);
483:   d->calcpairs_proj_res(d,r_s,r_e,R);
484:   SlepcVecPoolRestoreVecs(d->auxV,r_e-r_s,&R);
485:   return(0);
486: }

488: static PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R)
489: {
490:   PetscInt       i,l,k;
492:   PetscBool      lindep=PETSC_FALSE;
493:   BV             cX;

496:   if (d->W) cX = d->W; /* If left subspace exists, R <- orth(cY, R), nR[i] <- ||R[i]|| */
497:   else if (!(DVD_IS(d->sEP, DVD_EP_STD) && DVD_IS(d->sEP, DVD_EP_HERMITIAN))) cX = d->eps->V; /* If not HEP, R <- orth(cX, R), nR[i] <- ||R[i]|| */
498:   else cX = NULL; /* Otherwise, nR[i] <- ||R[i]|| */

500:   if (cX) {
501:     BVGetActiveColumns(cX,&l,&k);
502:     BVSetActiveColumns(cX,0,l);
503:     for (i=0;i<r_e-r_s;i++) {
504:       BVOrthogonalizeVec(cX,R[i],NULL,&d->nR[r_s+i],&lindep);
505:     }
506:     BVSetActiveColumns(cX,l,k);
507:     if (lindep || (PetscAbs(d->nR[r_s+i]) < PETSC_MACHINE_EPSILON)) {
508:       PetscInfo2(d->eps,"The computed eigenvector residual %D is too low, %g!\n",r_s+i,(double)(d->nR[r_s+i]));
509:     }
510:   } else {
511:     for (i=0;i<r_e-r_s;i++) {
512:       VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]);
513:     }
514:     for (i=0;i<r_e-r_s;i++) {
515:       VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]);
516:     }
517:   }
518:   return(0);
519: }

521: PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d,dvdBlackboard *b,PetscBool borth,PetscBool harm)
522: {
524:   PetscBool      std_probl,her_probl,ind_probl;
525:   DSType         dstype;
526:   Vec            v1;

529:   std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
530:   her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
531:   ind_probl = DVD_IS(d->sEP,DVD_EP_INDEFINITE)? PETSC_TRUE: PETSC_FALSE;

533:   /* Setting configuration constrains */
534:   b->max_size_proj = PetscMax(b->max_size_proj,b->max_size_V);
535:   d->W_shift = d->B? PETSC_TRUE: PETSC_FALSE;

537:   /* Setup the step */
538:   if (b->state >= DVD_STATE_CONF) {
539:     d->max_size_P = b->max_size_P;
540:     d->max_size_proj = b->max_size_proj;
541:     /* Create a DS if the method works with Schur decompositions */
542:     d->calcPairs = dvd_calcpairs_proj;
543:     d->calcpairs_residual = dvd_calcpairs_res_0;
544:     d->calcpairs_proj_res = dvd_calcpairs_proj_res;
545:     d->calcpairs_selectPairs = dvd_calcpairs_selectPairs;
546:     /* Create and configure a DS for solving the projected problems */
547:     if (d->W) dstype = DSGNHEP;    /* If we use harmonics */
548:     else {
549:       if (ind_probl) dstype = DSGHIEP;
550:       else if (std_probl) dstype = her_probl? DSHEP : DSNHEP;
551:       else dstype = her_probl? DSGHEP : DSGNHEP;
552:     }
553:     DSSetType(d->eps->ds,dstype);
554:     DSAllocate(d->eps->ds,d->eps->ncv);
555:     /* Create various vector basis */
556:     if (harm) {
557:       BVDuplicateResize(d->eps->V,d->eps->ncv,&d->W);
558:       BVSetMatrix(d->W,NULL,PETSC_FALSE);
559:     } else d->W = NULL;
560:     BVDuplicateResize(d->eps->V,d->eps->ncv,&d->AX);
561:     BVSetMatrix(d->AX,NULL,PETSC_FALSE);
562:     BVDuplicateResize(d->eps->V,d->eps->ncv,&d->auxBV);
563:     BVSetMatrix(d->auxBV,NULL,PETSC_FALSE);
564:     if (d->B) {
565:       BVDuplicateResize(d->eps->V,d->eps->ncv,&d->BX);
566:       BVSetMatrix(d->BX,NULL,PETSC_FALSE);
567:     } else d->BX = NULL;
568:     MatCreateVecsEmpty(d->A,&v1,NULL);
569:     SlepcVecPoolCreate(v1,0,&d->auxV);
570:     VecDestroy(&v1);
571:     /* Create projected problem matrices */
572:     MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->H);
573:     if (!std_probl) {
574:       MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->G);
575:     } else d->G = NULL;
576:     if (her_probl) {
577:       MatSetOption(d->H,MAT_HERMITIAN,PETSC_TRUE);
578:       if (d->G) { MatSetOption(d->G,MAT_HERMITIAN,PETSC_TRUE); }
579:     }

581:     if (ind_probl) {
582:       PetscMalloc1(d->eps->ncv,&d->nBds);
583:     } else d->nBds = NULL;
584:     MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->auxM);

586:     EPSDavidsonFLAdd(&d->startList,dvd_calcpairs_qz_start);
587:     EPSDavidsonFLAdd(&d->endList,EPSXDComputeDSConv);
588:     EPSDavidsonFLAdd(&d->destroyList,dvd_calcpairs_qz_d);
589:   }
590:   return(0);
591: }