Actual source code: dvdcalcpairs.c

slepc-3.23.1 2025-05-01
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: 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: {
 28:   PetscFunctionBegin;
 29:   PetscCall(BVSetActiveColumns(d->eps->V,0,0));
 30:   if (d->W) PetscCall(BVSetActiveColumns(d->W,0,0));
 31:   PetscCall(BVSetActiveColumns(d->AX,0,0));
 32:   if (d->BX) PetscCall(BVSetActiveColumns(d->BX,0,0));
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: static PetscErrorCode dvd_calcpairs_qz_d(dvdDashboard *d)
 37: {
 38:   PetscFunctionBegin;
 39:   PetscCall(BVDestroy(&d->W));
 40:   PetscCall(BVDestroy(&d->AX));
 41:   PetscCall(BVDestroy(&d->BX));
 42:   PetscCall(BVDestroy(&d->auxBV));
 43:   PetscCall(MatDestroy(&d->H));
 44:   PetscCall(MatDestroy(&d->G));
 45:   PetscCall(MatDestroy(&d->auxM));
 46:   PetscCall(SlepcVecPoolDestroy(&d->auxV));
 47:   PetscCall(PetscFree(d->nBds));
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

 51: /* in complex, d->size_H real auxiliary values are needed */
 52: static PetscErrorCode dvd_calcpairs_projeig_solve(dvdDashboard *d)
 53: {
 54:   Vec               v;
 55:   Mat               A,B,H0,G0;
 56:   PetscScalar       *pA;
 57:   const PetscScalar *pv;
 58:   PetscInt          i,lV,kV,n,ld;

 60:   PetscFunctionBegin;
 61:   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
 62:   n = kV-lV;
 63:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
 64:   PetscCall(DSSetDimensions(d->eps->ds,n,0,0));
 65:   PetscCall(DSGetMat(d->eps->ds,DS_MAT_A,&A));
 66:   PetscCall(MatDenseGetSubMatrix(d->H,lV,lV+n,lV,lV+n,&H0));
 67:   PetscCall(MatCopy(H0,A,SAME_NONZERO_PATTERN));
 68:   PetscCall(MatDenseRestoreSubMatrix(d->H,&H0));
 69:   PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_A,&A));
 70:   if (d->G) {
 71:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_B,&B));
 72:     PetscCall(MatDenseGetSubMatrix(d->G,lV,lV+n,lV,lV+n,&G0));
 73:     PetscCall(MatCopy(G0,B,SAME_NONZERO_PATTERN));
 74:     PetscCall(MatDenseRestoreSubMatrix(d->G,&G0));
 75:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_B,&B));
 76:   }
 77:   /* Set the signature on projected matrix B */
 78:   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
 79:     PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
 80:     PetscCall(DSGetArray(d->eps->ds,DS_MAT_B,&pA));
 81:     PetscCall(PetscArrayzero(pA,n*ld));
 82:     PetscCall(VecCreateSeq(PETSC_COMM_SELF,kV,&v));
 83:     PetscCall(BVGetSignature(d->eps->V,v));
 84:     PetscCall(VecGetArrayRead(v,&pv));
 85:     for (i=0;i<n;i++) {
 86:       pA[i+ld*i] = d->nBds[i] = PetscRealPart(pv[lV+i]);
 87:     }
 88:     PetscCall(VecRestoreArrayRead(v,&pv));
 89:     PetscCall(VecDestroy(&v));
 90:     PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_B,&pA));
 91:   }
 92:   PetscCall(DSSetState(d->eps->ds,DS_STATE_RAW));
 93:   PetscCall(DSSolve(d->eps->ds,d->eigr,d->eigi));
 94:   PetscFunctionReturn(PETSC_SUCCESS);
 95: }

 97: /*
 98:    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
 99:  */
100: static PetscErrorCode EPSXDUpdateProj(Mat Q,Mat Z,PetscInt l,Mat A,PetscInt lA,PetscInt kA,Mat aux)
101: {
102:   PetscScalar       one=1.0,zero=0.0;
103:   PetscInt          i,j,dA_=kA-lA,m0,n0,ldA_,ldQ_,ldZ_,nQ_;
104:   PetscBLASInt      dA,nQ,ldA,ldQ,ldZ;
105:   PetscScalar       *pA,*pW;
106:   const PetscScalar *pQ,*pZ;
107:   PetscBool         symm=PETSC_FALSE,set,flg;

109:   PetscFunctionBegin;
110:   PetscCall(MatGetSize(A,&m0,&n0));
111:   PetscCall(MatDenseGetLDA(A,&ldA_));
112:   PetscAssert(m0==n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"A should be square");
113:   PetscAssert(lA>=0 && lA<=m0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial row, column in A");
114:   PetscAssert(kA>=0 && kA>=lA && kA<=m0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid final row, column in A");
115:   PetscCall(MatIsHermitianKnown(A,&set,&flg));
116:   symm = set? flg: PETSC_FALSE;
117:   PetscCall(MatGetSize(Q,&m0,&n0)); nQ_=m0;
118:   PetscCall(MatDenseGetLDA(Q,&ldQ_));
119:   PetscAssert(l>=0 && l<=n0 && l+dA_<=n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Q");
120:   PetscCall(MatGetSize(Z,&m0,&n0));
121:   PetscCall(MatDenseGetLDA(Z,&ldZ_));
122:   PetscAssert(l>=0 && l<=n0 && l+dA_<=n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Z");
123:   PetscCall(MatGetSize(aux,&m0,&n0));
124:   PetscAssert(m0*n0>=nQ_*dA_,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"aux should be larger");
125:   PetscCall(PetscBLASIntCast(dA_,&dA));
126:   PetscCall(PetscBLASIntCast(nQ_,&nQ));
127:   PetscCall(PetscBLASIntCast(ldA_,&ldA));
128:   PetscCall(PetscBLASIntCast(ldQ_,&ldQ));
129:   PetscCall(PetscBLASIntCast(ldZ_,&ldZ));
130:   PetscCall(MatDenseGetArray(A,&pA));
131:   PetscCall(MatDenseGetArrayRead(Q,&pQ));
132:   if (Q!=Z) PetscCall(MatDenseGetArrayRead(Z,&pZ));
133:   else pZ = pQ;
134:   PetscCall(MatDenseGetArrayWrite(aux,&pW));
135:   /* W = A*Q */
136:   if (symm) {
137:     /* symmetrize before multiplying */
138:     for (i=lA+1;i<lA+nQ;i++) {
139:       for (j=lA;j<i;j++) pA[i+j*ldA] = PetscConj(pA[j+i*ldA]);
140:     }
141:   }
142:   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nQ,&dA,&nQ,&one,&pA[ldA*lA+lA],&ldA,&pQ[ldQ*l+l],&ldQ,&zero,pW,&nQ));
143:   /* A = Q'*W */
144:   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&dA,&dA,&nQ,&one,&pZ[ldZ*l+l],&ldZ,pW,&nQ,&zero,&pA[ldA*lA+lA],&ldA));
145:   PetscCall(MatDenseRestoreArray(A,&pA));
146:   PetscCall(MatDenseRestoreArrayRead(Q,&pQ));
147:   if (Q!=Z) PetscCall(MatDenseRestoreArrayRead(Z,&pZ));
148:   PetscCall(MatDenseRestoreArrayWrite(aux,&pW));
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

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

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

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

197: /*
198:    BV <- BV*MT
199:  */
200: static inline PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,BV bv,DSMatType mat)
201: {
202:   PetscInt       l,k,n;
203:   Mat            M,M0,auxM,auxM0;

205:   PetscFunctionBegin;
206:   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
207:   PetscCall(DSGetDimensions(d->eps->ds,&n,NULL,NULL,NULL));
208:   PetscAssert(k-l==n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
209:   PetscCall(DSGetMat(d->eps->ds,mat,&M));
210:   PetscCall(MatDenseGetSubMatrix(M,0,n,0,d->V_tra_e,&M0));
211:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&auxM));
212:   PetscCall(MatDenseGetSubMatrix(auxM,l,l+n,l,l+d->V_tra_e,&auxM0));
213:   PetscCall(MatCopy(M0,auxM0,SAME_NONZERO_PATTERN));
214:   PetscCall(MatDenseRestoreSubMatrix(auxM,&auxM0));
215:   PetscCall(MatDenseRestoreSubMatrix(M,&M0));
216:   PetscCall(DSRestoreMat(d->eps->ds,mat,&M));
217:   PetscCall(BVMultInPlace(bv,auxM,l,l+d->V_tra_e));
218:   PetscCall(MatDestroy(&auxM));
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: static PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
223: {
224:   PetscInt       i,l,k;
225:   Vec            v1,v2;
226:   PetscScalar    *pv;

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

294:     /* Perform the transformation on the projected problem */
295:     if (d->W) PetscCall(d->calcpairs_proj_trans(d));
296:     k = l+d->V_new_e;
297:   }
298:   PetscCall(BVSetActiveColumns(d->eps->V,l,k));
299:   PetscCall(BVSetActiveColumns(d->AX,l,k));
300:   if (d->BX) PetscCall(BVSetActiveColumns(d->BX,l,k));
301:   if (d->W) PetscCall(BVSetActiveColumns(d->W,l,k));

303:   /* Solve the projected problem */
304:   PetscCall(dvd_calcpairs_projeig_solve(d));

306:   d->V_tra_s = d->V_tra_e = 0;
307:   d->V_new_s = d->V_new_e;
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }

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

322:   PetscFunctionBegin;
323:   /* Quick exit without neither arbitrary selection nor harmonic extraction */
324:   if (!d->eps->arbitrary && !d->calcpairs_eig_backtrans) PetscFunctionReturn(PETSC_SUCCESS);

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

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

365: static PetscErrorCode dvd_calcpairs_selectPairs(dvdDashboard *d,PetscInt n)
366: {
367:   PetscInt       k,lV,kV,nV;
368:   PetscScalar    *rr,*ri;

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

399:   if (d->calcpairs_eigs_trans) PetscCall(d->calcpairs_eigs_trans(d));
400:   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
401:     PetscCall(PetscFree(rr));
402:     PetscCall(PetscFree(ri));
403:   }
404:   PetscFunctionReturn(PETSC_SUCCESS);
405: }

407: static PetscErrorCode EPSXDComputeDSConv(dvdDashboard *d)
408: {
409:   PetscInt          i,ld;
410:   Vec               v;
411:   Mat               A,B,H0,G0;
412:   PetscScalar       *pA;
413:   const PetscScalar *pv;
414:   PetscBool         symm;

416:   PetscFunctionBegin;
417:   PetscCall(BVSetActiveColumns(d->eps->V,0,d->eps->nconv));
418:   PetscCall(PetscObjectTypeCompare((PetscObject)d->eps->ds,DSHEP,&symm));
419:   if (symm || !d->eps->nconv) PetscFunctionReturn(PETSC_SUCCESS);
420:   PetscCall(DSSetDimensions(d->eps->ds,d->eps->nconv,0,0));
421:   PetscCall(DSGetMat(d->eps->ds,DS_MAT_A,&A));
422:   PetscCall(MatDenseGetSubMatrix(d->H,0,d->eps->nconv,0,d->eps->nconv,&H0));
423:   PetscCall(MatCopy(H0,A,SAME_NONZERO_PATTERN));
424:   PetscCall(MatDenseRestoreSubMatrix(d->H,&H0));
425:   PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_A,&A));
426:   if (d->G) {
427:     PetscCall(DSGetMat(d->eps->ds,DS_MAT_B,&B));
428:     PetscCall(MatDenseGetSubMatrix(d->G,0,d->eps->nconv,0,d->eps->nconv,&G0));
429:     PetscCall(MatCopy(G0,B,SAME_NONZERO_PATTERN));
430:     PetscCall(MatDenseRestoreSubMatrix(d->G,&G0));
431:     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_B,&B));
432:   }
433:   /* Set the signature on projected matrix B */
434:   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
435:     PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
436:     PetscCall(DSGetArray(d->eps->ds,DS_MAT_B,&pA));
437:     PetscCall(PetscArrayzero(pA,d->eps->nconv*ld));
438:     PetscCall(VecCreateSeq(PETSC_COMM_SELF,d->eps->nconv,&v));
439:     PetscCall(BVGetSignature(d->eps->V,v));
440:     PetscCall(VecGetArrayRead(v,&pv));
441:     for (i=0;i<d->eps->nconv;i++) pA[i+ld*i] = pv[i];
442:     PetscCall(VecRestoreArrayRead(v,&pv));
443:     PetscCall(VecDestroy(&v));
444:     PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_B,&pA));
445:   }
446:   PetscCall(DSSetState(d->eps->ds,DS_STATE_RAW));
447:   PetscCall(DSSolve(d->eps->ds,d->eps->eigr,d->eps->eigi));
448:   PetscCall(DSSynchronize(d->eps->ds,d->eps->eigr,d->eps->eigi));
449:   if (d->W) {
450:     for (i=0;i<d->eps->nconv;i++) PetscCall(d->calcpairs_eig_backtrans(d,d->eps->eigr[i],d->eps->eigi[i],&d->eps->eigr[i],&d->eps->eigi[i]));
451:   }
452:   PetscFunctionReturn(PETSC_SUCCESS);
453: }

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

466:   PetscFunctionBegin;
467:   PetscCall(DSGetLeadingDimension(d->eps->ds,&ldpX));
468:   PetscCall(DSGetArray(d->eps->ds,DS_MAT_Q,&pX));
469:   /* nX(i) <- ||X(i)|| */
470:   PetscCall(dvd_improvex_compute_X(d,r_s,r_e,NULL,pX,ldpX));
471:   PetscCall(SlepcVecPoolGetVecs(d->auxV,r_e-r_s,&R));
472:   for (i=r_s;i<r_e;i++) {
473:     /* R(i-r_s) <- AV*pX(i) */
474:     PetscCall(BVMultVec(d->AX,1.0,0.0,R[i-r_s],&pX[ldpX*i]));
475:     /* R(i-r_s) <- R(i-r_s) - eigr(i)*BX*pX(i) */
476:     PetscCall(BVMultVec(BX,-d->eigr[i],1.0,R[i-r_s],&pX[ldpX*i]));
477:   }
478:   PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_Q,&pX));
479:   PetscCall(d->calcpairs_proj_res(d,r_s,r_e,R));
480:   PetscCall(SlepcVecPoolRestoreVecs(d->auxV,r_e-r_s,&R));
481:   PetscFunctionReturn(PETSC_SUCCESS);
482: }

484: static PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R)
485: {
486:   PetscInt       i,l,k;
487:   PetscBool      lindep=PETSC_FALSE;
488:   BV             cX;

490:   PetscFunctionBegin;
491:   if (d->W) cX = d->W; /* If left subspace exists, R <- orth(cY, R), nR[i] <- ||R[i]|| */
492:   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]|| */
493:   else cX = NULL; /* Otherwise, nR[i] <- ||R[i]|| */

495:   if (cX) {
496:     PetscCall(BVGetActiveColumns(cX,&l,&k));
497:     PetscCall(BVSetActiveColumns(cX,0,l));
498:     for (i=0;i<r_e-r_s;i++) PetscCall(BVOrthogonalizeVec(cX,R[i],NULL,&d->nR[r_s+i],&lindep));
499:     PetscCall(BVSetActiveColumns(cX,l,k));
500:     if (lindep || (PetscAbs(d->nR[r_s+i]) < PETSC_MACHINE_EPSILON)) PetscCall(PetscInfo(d->eps,"The computed eigenvector residual %" PetscInt_FMT " is too low, %g!\n",r_s+i,(double)d->nR[r_s+i]));
501:   } else {
502:     for (i=0;i<r_e-r_s;i++) PetscCall(VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]));
503:     for (i=0;i<r_e-r_s;i++) PetscCall(VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]));
504:   }
505:   PetscFunctionReturn(PETSC_SUCCESS);
506: }

508: PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d,dvdBlackboard *b,PetscBool borth,PetscBool harm)
509: {
510:   PetscBool      std_probl,her_probl,ind_probl;
511:   DSType         dstype;
512:   Vec            v1;

514:   PetscFunctionBegin;
515:   std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
516:   her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
517:   ind_probl = DVD_IS(d->sEP,DVD_EP_INDEFINITE)? PETSC_TRUE: PETSC_FALSE;

519:   /* Setting configuration constrains */
520:   b->max_size_proj = PetscMax(b->max_size_proj,b->max_size_V);
521:   d->W_shift = d->B? PETSC_TRUE: PETSC_FALSE;

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

566:     if (ind_probl) PetscCall(PetscMalloc1(d->eps->ncv,&d->nBds));
567:     else d->nBds = NULL;
568:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->auxM));

570:     PetscCall(EPSDavidsonFLAdd(&d->startList,dvd_calcpairs_qz_start));
571:     PetscCall(EPSDavidsonFLAdd(&d->endList,EPSXDComputeDSConv));
572:     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_calcpairs_qz_d));
573:   }
574:   PetscFunctionReturn(PETSC_SUCCESS);
575: }