Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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"
      12             : 
      13             :    Step: calculate the best eigenpairs in the subspace V
      14             : 
      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             : */
      22             : 
      23             : #include "davidson.h"
      24             : #include <slepcblaslapack.h>
      25             : 
      26          95 : static PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
      27             : {
      28          95 :   PetscFunctionBegin;
      29          95 :   PetscCall(BVSetActiveColumns(d->eps->V,0,0));
      30          95 :   if (d->W) PetscCall(BVSetActiveColumns(d->W,0,0));
      31          95 :   PetscCall(BVSetActiveColumns(d->AX,0,0));
      32          95 :   if (d->BX) PetscCall(BVSetActiveColumns(d->BX,0,0));
      33          95 :   PetscFunctionReturn(PETSC_SUCCESS);
      34             : }
      35             : 
      36          95 : static PetscErrorCode dvd_calcpairs_qz_d(dvdDashboard *d)
      37             : {
      38          95 :   PetscFunctionBegin;
      39          95 :   PetscCall(BVDestroy(&d->W));
      40          95 :   PetscCall(BVDestroy(&d->AX));
      41          95 :   PetscCall(BVDestroy(&d->BX));
      42          95 :   PetscCall(BVDestroy(&d->auxBV));
      43          95 :   PetscCall(MatDestroy(&d->H));
      44          95 :   PetscCall(MatDestroy(&d->G));
      45          95 :   PetscCall(MatDestroy(&d->auxM));
      46          95 :   PetscCall(SlepcVecPoolDestroy(&d->auxV));
      47          95 :   PetscCall(PetscFree(d->nBds));
      48          95 :   PetscFunctionReturn(PETSC_SUCCESS);
      49             : }
      50             : 
      51             : /* in complex, d->size_H real auxiliary values are needed */
      52        7946 : static PetscErrorCode dvd_calcpairs_projeig_solve(dvdDashboard *d)
      53             : {
      54        7946 :   Vec               v;
      55        7946 :   Mat               A,B,H0,G0;
      56        7946 :   PetscScalar       *pA;
      57        7946 :   const PetscScalar *pv;
      58        7946 :   PetscInt          i,lV,kV,n,ld;
      59             : 
      60        7946 :   PetscFunctionBegin;
      61        7946 :   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
      62        7946 :   n = kV-lV;
      63        7946 :   PetscCall(DSSetDimensions(d->eps->ds,n,0,0));
      64        7946 :   PetscCall(DSGetMat(d->eps->ds,DS_MAT_A,&A));
      65        7946 :   PetscCall(MatDenseGetSubMatrix(d->H,lV,lV+n,lV,lV+n,&H0));
      66        7946 :   PetscCall(MatCopy(H0,A,SAME_NONZERO_PATTERN));
      67        7946 :   PetscCall(MatDenseRestoreSubMatrix(d->H,&H0));
      68        7946 :   PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_A,&A));
      69        7946 :   if (d->G) {
      70        1135 :     PetscCall(DSGetMat(d->eps->ds,DS_MAT_B,&B));
      71        1135 :     PetscCall(MatDenseGetSubMatrix(d->G,lV,lV+n,lV,lV+n,&G0));
      72        1135 :     PetscCall(MatCopy(G0,B,SAME_NONZERO_PATTERN));
      73        1135 :     PetscCall(MatDenseRestoreSubMatrix(d->G,&G0));
      74        1135 :     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_B,&B));
      75             :   }
      76             :   /* Set the signature on projected matrix B */
      77        7946 :   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
      78           0 :     PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
      79           0 :     PetscCall(DSGetArray(d->eps->ds,DS_MAT_B,&pA));
      80           0 :     PetscCall(PetscArrayzero(pA,n*ld));
      81           0 :     PetscCall(VecCreateSeq(PETSC_COMM_SELF,kV,&v));
      82           0 :     PetscCall(BVGetSignature(d->eps->V,v));
      83           0 :     PetscCall(VecGetArrayRead(v,&pv));
      84           0 :     for (i=0;i<n;i++) {
      85           0 :       pA[i+ld*i] = d->nBds[i] = PetscRealPart(pv[lV+i]);
      86             :     }
      87           0 :     PetscCall(VecRestoreArrayRead(v,&pv));
      88           0 :     PetscCall(VecDestroy(&v));
      89           0 :     PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_B,&pA));
      90             :   }
      91        7946 :   PetscCall(DSSetState(d->eps->ds,DS_STATE_RAW));
      92        7946 :   PetscCall(DSSolve(d->eps->ds,d->eigr,d->eigi));
      93        7946 :   PetscFunctionReturn(PETSC_SUCCESS);
      94             : }
      95             : 
      96             : /*
      97             :    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
      98             :  */
      99        1631 : static PetscErrorCode EPSXDUpdateProj(Mat Q,Mat Z,PetscInt l,Mat A,PetscInt lA,PetscInt kA,Mat aux)
     100             : {
     101        1631 :   PetscScalar       one=1.0,zero=0.0;
     102        1631 :   PetscInt          i,j,dA_=kA-lA,m0,n0,ldA_,ldQ_,ldZ_,nQ_;
     103        1631 :   PetscBLASInt      dA,nQ,ldA,ldQ,ldZ;
     104        1631 :   PetscScalar       *pA,*pW;
     105        1631 :   const PetscScalar *pQ,*pZ;
     106        1631 :   PetscBool         symm=PETSC_FALSE,set,flg;
     107             : 
     108        1631 :   PetscFunctionBegin;
     109        1631 :   PetscCall(MatGetSize(A,&m0,&n0));
     110        1631 :   PetscCall(MatDenseGetLDA(A,&ldA_));
     111        1631 :   PetscAssert(m0==n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"A should be square");
     112        1631 :   PetscAssert(lA>=0 && lA<=m0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial row, column in A");
     113        1631 :   PetscAssert(kA>=0 && kA>=lA && kA<=m0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid final row, column in A");
     114        1631 :   PetscCall(MatIsHermitianKnown(A,&set,&flg));
     115        1631 :   symm = set? flg: PETSC_FALSE;
     116        1631 :   PetscCall(MatGetSize(Q,&m0,&n0)); nQ_=m0;
     117        1631 :   PetscCall(MatDenseGetLDA(Q,&ldQ_));
     118        1631 :   PetscAssert(l>=0 && l<=n0 && l+dA_<=n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Q");
     119        1631 :   PetscCall(MatGetSize(Z,&m0,&n0));
     120        1631 :   PetscCall(MatDenseGetLDA(Z,&ldZ_));
     121        1631 :   PetscAssert(l>=0 && l<=n0 && l+dA_<=n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Z");
     122        1631 :   PetscCall(MatGetSize(aux,&m0,&n0));
     123        1631 :   PetscAssert(m0*n0>=nQ_*dA_,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"aux should be larger");
     124        1631 :   PetscCall(PetscBLASIntCast(dA_,&dA));
     125        1631 :   PetscCall(PetscBLASIntCast(nQ_,&nQ));
     126        1631 :   PetscCall(PetscBLASIntCast(ldA_,&ldA));
     127        1631 :   PetscCall(PetscBLASIntCast(ldQ_,&ldQ));
     128        1631 :   PetscCall(PetscBLASIntCast(ldZ_,&ldZ));
     129        1631 :   PetscCall(MatDenseGetArray(A,&pA));
     130        1631 :   PetscCall(MatDenseGetArrayRead(Q,&pQ));
     131        1631 :   if (Q!=Z) PetscCall(MatDenseGetArrayRead(Z,&pZ));
     132        1243 :   else pZ = pQ;
     133        1631 :   PetscCall(MatDenseGetArrayWrite(aux,&pW));
     134             :   /* W = A*Q */
     135        1631 :   if (symm) {
     136             :     /* symmetrize before multiplying */
     137       15299 :     for (i=lA+1;i<lA+nQ;i++) {
     138      117117 :       for (j=lA;j<i;j++) pA[i+j*ldA] = PetscConj(pA[j+i*ldA]);
     139             :     }
     140             :   }
     141        1631 :   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nQ,&dA,&nQ,&one,&pA[ldA*lA+lA],&ldA,&pQ[ldQ*l+l],&ldQ,&zero,pW,&nQ));
     142             :   /* A = Q'*W */
     143        1631 :   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&dA,&dA,&nQ,&one,&pZ[ldZ*l+l],&ldZ,pW,&nQ,&zero,&pA[ldA*lA+lA],&ldA));
     144        1631 :   PetscCall(MatDenseRestoreArray(A,&pA));
     145        1631 :   PetscCall(MatDenseRestoreArrayRead(Q,&pQ));
     146        1631 :   if (Q!=Z) PetscCall(MatDenseRestoreArrayRead(Z,&pZ));
     147        1631 :   PetscCall(MatDenseRestoreArrayWrite(aux,&pW));
     148        1631 :   PetscFunctionReturn(PETSC_SUCCESS);
     149             : }
     150             : 
     151        1437 : static PetscErrorCode dvd_calcpairs_updateproj(dvdDashboard *d)
     152             : {
     153        1437 :   Mat            Q,Z;
     154        1437 :   PetscInt       lV,kV;
     155        1437 :   PetscBool      symm;
     156             : 
     157        1437 :   PetscFunctionBegin;
     158        1437 :   PetscCall(DSGetMat(d->eps->ds,DS_MAT_Q,&Q));
     159        1437 :   if (d->W) PetscCall(DSGetMat(d->eps->ds,DS_MAT_Z,&Z));
     160        1243 :   else Z = Q;
     161        1437 :   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
     162        1437 :   PetscCall(EPSXDUpdateProj(Q,Z,0,d->H,lV,lV+d->V_tra_e,d->auxM));
     163        1437 :   if (d->G) PetscCall(EPSXDUpdateProj(Q,Z,0,d->G,lV,lV+d->V_tra_e,d->auxM));
     164        1437 :   PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q));
     165        1437 :   if (d->W) PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z));
     166             : 
     167        1437 :   PetscCall(PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSHEP,DSGHIEP,DSGHEP,""));
     168        1437 :   if (d->V_tra_s==0 || symm) PetscFunctionReturn(PETSC_SUCCESS);
     169             :   /* Compute upper part of H (and G): H(0:l-1,l:k-1) <- W(0:l-1)' * AV(l:k-1), where
     170             :      k=l+d->V_tra_s */
     171          67 :   PetscCall(BVSetActiveColumns(d->W?d->W:d->eps->V,0,lV));
     172          67 :   PetscCall(BVSetActiveColumns(d->AX,lV,lV+d->V_tra_s));
     173          67 :   PetscCall(BVDot(d->AX,d->W?d->W:d->eps->V,d->H));
     174          67 :   if (d->G) {
     175          55 :     PetscCall(BVSetActiveColumns(d->BX?d->BX:d->eps->V,lV,lV+d->V_tra_s));
     176          55 :     PetscCall(BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G));
     177             :   }
     178          67 :   PetscCall(PetscObjectTypeCompare((PetscObject)d->eps->ds,DSGHEP,&symm));
     179          67 :   if (!symm) {
     180          67 :     PetscCall(BVSetActiveColumns(d->W?d->W:d->eps->V,lV,lV+d->V_tra_s));
     181          67 :     PetscCall(BVSetActiveColumns(d->AX,0,lV));
     182          67 :     PetscCall(BVDot(d->AX,d->W?d->W:d->eps->V,d->H));
     183          67 :     if (d->G) {
     184          55 :       PetscCall(BVSetActiveColumns(d->BX?d->BX:d->eps->V,0,lV));
     185          55 :       PetscCall(BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G));
     186             :     }
     187             :   }
     188          67 :   PetscCall(BVSetActiveColumns(d->eps->V,lV,kV));
     189          67 :   PetscCall(BVSetActiveColumns(d->AX,lV,kV));
     190          67 :   if (d->BX) PetscCall(BVSetActiveColumns(d->BX,lV,kV));
     191          67 :   if (d->W) PetscCall(BVSetActiveColumns(d->W,lV,kV));
     192          67 :   if (d->W) PetscCall(dvd_harm_updateproj(d));
     193          67 :   PetscFunctionReturn(PETSC_SUCCESS);
     194             : }
     195             : 
     196             : /*
     197             :    BV <- BV*MT
     198             :  */
     199        3281 : static inline PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,BV bv,DSMatType mat)
     200             : {
     201        3281 :   PetscInt       l,k,n;
     202        3281 :   Mat            M,M0,auxM,auxM0;
     203             : 
     204        3281 :   PetscFunctionBegin;
     205        3281 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
     206        3281 :   PetscCall(DSGetDimensions(d->eps->ds,&n,NULL,NULL,NULL));
     207        3281 :   PetscAssert(k-l==n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
     208        3281 :   PetscCall(DSGetMat(d->eps->ds,mat,&M));
     209        3281 :   PetscCall(MatDenseGetSubMatrix(M,0,n,0,d->V_tra_e,&M0));
     210        3281 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&auxM));
     211        3281 :   PetscCall(MatDenseGetSubMatrix(auxM,l,l+n,l,l+d->V_tra_e,&auxM0));
     212        3281 :   PetscCall(MatCopy(M0,auxM0,SAME_NONZERO_PATTERN));
     213        3281 :   PetscCall(MatDenseRestoreSubMatrix(auxM,&auxM0));
     214        3281 :   PetscCall(MatDenseRestoreSubMatrix(M,&M0));
     215        3281 :   PetscCall(DSRestoreMat(d->eps->ds,mat,&M));
     216        3281 :   PetscCall(BVMultInPlace(bv,auxM,l,l+d->V_tra_e));
     217        3281 :   PetscCall(MatDestroy(&auxM));
     218        3281 :   PetscFunctionReturn(PETSC_SUCCESS);
     219             : }
     220             : 
     221        7946 : static PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
     222             : {
     223        7946 :   PetscInt       i,l,k;
     224        7946 :   Vec            v1,v2;
     225        7946 :   PetscScalar    *pv;
     226             : 
     227        7946 :   PetscFunctionBegin;
     228        7946 :   PetscCall(BVGetActiveColumns(d->eps->V,&l,&k));
     229             :   /* Update AV, BV, W and the projected matrices */
     230             :   /* 1. S <- S*MT */
     231        7946 :   if (d->V_tra_s != d->V_tra_e || d->V_tra_e > 0) {
     232        1437 :     PetscCall(dvd_calcpairs_updateBV0_gen(d,d->eps->V,DS_MAT_Q));
     233        1437 :     if (d->W) PetscCall(dvd_calcpairs_updateBV0_gen(d,d->W,DS_MAT_Z));
     234        1437 :     PetscCall(dvd_calcpairs_updateBV0_gen(d,d->AX,DS_MAT_Q));
     235        1437 :     if (d->BX) PetscCall(dvd_calcpairs_updateBV0_gen(d,d->BX,DS_MAT_Q));
     236        1437 :     PetscCall(dvd_calcpairs_updateproj(d));
     237             :     /* Update signature */
     238        1437 :     if (d->nBds) {
     239           0 :       PetscCall(VecCreateSeq(PETSC_COMM_SELF,l+d->V_tra_e,&v1));
     240           0 :       PetscCall(BVSetActiveColumns(d->eps->V,0,l+d->V_tra_e));
     241           0 :       PetscCall(BVGetSignature(d->eps->V,v1));
     242           0 :       PetscCall(VecGetArray(v1,&pv));
     243           0 :       for (i=0;i<d->V_tra_e;i++) pv[l+i] = d->nBds[i];
     244           0 :       PetscCall(VecRestoreArray(v1,&pv));
     245           0 :       PetscCall(BVSetSignature(d->eps->V,v1));
     246           0 :       PetscCall(BVSetActiveColumns(d->eps->V,l,k));
     247           0 :       PetscCall(VecDestroy(&v1));
     248             :     }
     249        1437 :     k = l+d->V_tra_e;
     250        1437 :     l+= d->V_tra_s;
     251             :   } else {
     252             :     /* 2. V <- orth(V, V_new) */
     253        6509 :     PetscCall(dvd_orthV(d->eps->V,l+d->V_new_s,l+d->V_new_e));
     254             :     /* 3. AV <- [AV A * V(V_new_s:V_new_e-1)] */
     255             :     /* Check consistency */
     256        6509 :     PetscAssert(k-l==d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
     257       16145 :     for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
     258        9636 :       PetscCall(BVGetColumn(d->eps->V,i,&v1));
     259        9636 :       PetscCall(BVGetColumn(d->AX,i,&v2));
     260        9636 :       PetscCall(MatMult(d->A,v1,v2));
     261        9636 :       PetscCall(BVRestoreColumn(d->eps->V,i,&v1));
     262        9636 :       PetscCall(BVRestoreColumn(d->AX,i,&v2));
     263             :     }
     264             :     /* 4. BV <- [BV B * V(V_new_s:V_new_e-1)] */
     265        6509 :     if (d->BX) {
     266             :       /* Check consistency */
     267        1236 :       PetscAssert(k-l==d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
     268        2728 :       for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
     269        1492 :         PetscCall(BVGetColumn(d->eps->V,i,&v1));
     270        1492 :         PetscCall(BVGetColumn(d->BX,i,&v2));
     271        1492 :         PetscCall(MatMult(d->B,v1,v2));
     272        1492 :         PetscCall(BVRestoreColumn(d->eps->V,i,&v1));
     273        1492 :         PetscCall(BVRestoreColumn(d->BX,i,&v2));
     274             :       }
     275             :     }
     276             :     /* 5. W <- [W f(AV,BV)] */
     277        6509 :     if (d->W) {
     278         941 :       PetscCall(d->calcpairs_W(d));
     279         941 :       PetscCall(dvd_orthV(d->W,l+d->V_new_s,l+d->V_new_e));
     280             :     }
     281             :     /* 6. H <- W' * AX; G <- W' * BX */
     282        6509 :     PetscCall(BVSetActiveColumns(d->eps->V,l+d->V_new_s,l+d->V_new_e));
     283        6509 :     PetscCall(BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e));
     284        6509 :     if (d->BX) PetscCall(BVSetActiveColumns(d->BX,l+d->V_new_s,l+d->V_new_e));
     285        6509 :     if (d->W) PetscCall(BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e));
     286        6509 :     PetscCall(BVMatProject(d->AX,NULL,d->W?d->W:d->eps->V,d->H));
     287        6509 :     if (d->G) PetscCall(BVMatProject(d->BX?d->BX:d->eps->V,NULL,d->W?d->W:d->eps->V,d->G));
     288        6509 :     PetscCall(BVSetActiveColumns(d->eps->V,l,k));
     289        6509 :     PetscCall(BVSetActiveColumns(d->AX,l,k));
     290        6509 :     if (d->BX) PetscCall(BVSetActiveColumns(d->BX,l,k));
     291        6509 :     if (d->W) PetscCall(BVSetActiveColumns(d->W,l,k));
     292             : 
     293             :     /* Perform the transformation on the projected problem */
     294        6509 :     if (d->W) PetscCall(d->calcpairs_proj_trans(d));
     295        6509 :     k = l+d->V_new_e;
     296             :   }
     297        7946 :   PetscCall(BVSetActiveColumns(d->eps->V,l,k));
     298        7946 :   PetscCall(BVSetActiveColumns(d->AX,l,k));
     299        7946 :   if (d->BX) PetscCall(BVSetActiveColumns(d->BX,l,k));
     300        7946 :   if (d->W) PetscCall(BVSetActiveColumns(d->W,l,k));
     301             : 
     302             :   /* Solve the projected problem */
     303        7946 :   PetscCall(dvd_calcpairs_projeig_solve(d));
     304             : 
     305        7946 :   d->V_tra_s = d->V_tra_e = 0;
     306        7946 :   d->V_new_s = d->V_new_e;
     307        7946 :   PetscFunctionReturn(PETSC_SUCCESS);
     308             : }
     309             : 
     310        9128 : static PetscErrorCode dvd_calcpairs_apply_arbitrary(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscScalar *rr,PetscScalar *ri)
     311             : {
     312        9128 :   PetscInt       i,k,ld;
     313        9128 :   PetscScalar    *pX;
     314        9128 :   Vec            *X,xr,xi;
     315             : #if defined(PETSC_USE_COMPLEX)
     316             :   PetscInt       N=1;
     317             : #else
     318        9128 :   PetscInt       N=2,j;
     319             : #endif
     320             : 
     321        9128 :   PetscFunctionBegin;
     322             :   /* Quick exit without neither arbitrary selection nor harmonic extraction */
     323        9128 :   if (!d->eps->arbitrary && !d->calcpairs_eig_backtrans) PetscFunctionReturn(PETSC_SUCCESS);
     324             : 
     325             :   /* Quick exit without arbitrary selection, but with harmonic extraction */
     326        2554 :   if (d->calcpairs_eig_backtrans) {
     327       26794 :     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]));
     328             :   }
     329        2554 :   if (!d->eps->arbitrary) PetscFunctionReturn(PETSC_SUCCESS);
     330             : 
     331        1192 :   PetscCall(SlepcVecPoolGetVecs(d->auxV,N,&X));
     332        1192 :   PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
     333       15116 :   for (i=r_s;i<r_e;i++) {
     334       13924 :     k = i;
     335       13924 :     PetscCall(DSVectors(d->eps->ds,DS_MAT_X,&k,NULL));
     336       13924 :     PetscCall(DSGetArray(d->eps->ds,DS_MAT_X,&pX));
     337       13924 :     PetscCall(dvd_improvex_compute_X(d,i,k+1,X,pX,ld));
     338       13924 :     PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_X,&pX));
     339             : #if !defined(PETSC_USE_COMPLEX)
     340       13924 :     if (d->nX[i] != 1.0) {
     341           0 :       for (j=i;j<k+1;j++) PetscCall(VecScale(X[j-i],1.0/d->nX[i]));
     342             :     }
     343       13924 :     xr = X[0];
     344       13924 :     xi = X[1];
     345       13924 :     if (i == k) PetscCall(VecSet(xi,0.0));
     346             : #else
     347             :     xr = X[0];
     348             :     xi = NULL;
     349             :     if (d->nX[i] != 1.0) PetscCall(VecScale(xr,1.0/d->nX[i]));
     350             : #endif
     351       13924 :     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));
     352             : #if !defined(PETSC_USE_COMPLEX)
     353       13924 :     if (i != k) {
     354           0 :       rr[i+1-r_s] = rr[i-r_s];
     355           0 :       ri[i+1-r_s] = ri[i-r_s];
     356           0 :       i++;
     357             :     }
     358             : #endif
     359             :   }
     360        1192 :   PetscCall(SlepcVecPoolRestoreVecs(d->auxV,N,&X));
     361        1192 :   PetscFunctionReturn(PETSC_SUCCESS);
     362             : }
     363             : 
     364        7851 : static PetscErrorCode dvd_calcpairs_selectPairs(dvdDashboard *d,PetscInt n)
     365             : {
     366        7851 :   PetscInt       k,lV,kV,nV;
     367        7851 :   PetscScalar    *rr,*ri;
     368             : 
     369        7851 :   PetscFunctionBegin;
     370        7851 :   PetscCall(BVGetActiveColumns(d->eps->V,&lV,&kV));
     371        7851 :   nV = kV - lV;
     372        7851 :   n = PetscMin(n,nV);
     373        7851 :   if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS);
     374             :   /* Put the best n pairs at the beginning. Useful for restarting */
     375        7851 :   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
     376        1277 :     PetscCall(PetscMalloc1(nV,&rr));
     377        1277 :     PetscCall(PetscMalloc1(nV,&ri));
     378        1277 :     PetscCall(dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri));
     379             :   } else {
     380        6574 :     rr = d->eigr;
     381        6574 :     ri = d->eigi;
     382             :   }
     383        7851 :   k = n;
     384        7851 :   PetscCall(DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k));
     385             :   /* Put the best pair at the beginning. Useful to check its residual */
     386             : #if !defined(PETSC_USE_COMPLEX)
     387        7851 :   if (n != 1 && (n != 2 || d->eigi[0] == 0.0))
     388             : #else
     389             :   if (n != 1)
     390             : #endif
     391             :   {
     392        7851 :     PetscCall(dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri));
     393        7851 :     k = 1;
     394        7851 :     PetscCall(DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k));
     395             :   }
     396        7851 :   PetscCall(DSSynchronize(d->eps->ds,d->eigr,d->eigi));
     397             : 
     398        7851 :   if (d->calcpairs_eigs_trans) PetscCall(d->calcpairs_eigs_trans(d));
     399        7851 :   if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
     400        1277 :     PetscCall(PetscFree(rr));
     401        1277 :     PetscCall(PetscFree(ri));
     402             :   }
     403        7851 :   PetscFunctionReturn(PETSC_SUCCESS);
     404             : }
     405             : 
     406          95 : static PetscErrorCode EPSXDComputeDSConv(dvdDashboard *d)
     407             : {
     408          95 :   PetscInt          i,ld;
     409          95 :   Vec               v;
     410          95 :   Mat               A,B,H0,G0;
     411          95 :   PetscScalar       *pA;
     412          95 :   const PetscScalar *pv;
     413          95 :   PetscBool         symm;
     414             : 
     415          95 :   PetscFunctionBegin;
     416          95 :   PetscCall(BVSetActiveColumns(d->eps->V,0,d->eps->nconv));
     417          95 :   PetscCall(PetscObjectTypeCompare((PetscObject)d->eps->ds,DSHEP,&symm));
     418          95 :   if (symm || !d->eps->nconv) PetscFunctionReturn(PETSC_SUCCESS);
     419          26 :   PetscCall(DSSetDimensions(d->eps->ds,d->eps->nconv,0,0));
     420          26 :   PetscCall(DSGetMat(d->eps->ds,DS_MAT_A,&A));
     421          26 :   PetscCall(MatDenseGetSubMatrix(d->H,0,d->eps->nconv,0,d->eps->nconv,&H0));
     422          26 :   PetscCall(MatCopy(H0,A,SAME_NONZERO_PATTERN));
     423          26 :   PetscCall(MatDenseRestoreSubMatrix(d->H,&H0));
     424          26 :   PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_A,&A));
     425          26 :   if (d->G) {
     426          23 :     PetscCall(DSGetMat(d->eps->ds,DS_MAT_B,&B));
     427          23 :     PetscCall(MatDenseGetSubMatrix(d->G,0,d->eps->nconv,0,d->eps->nconv,&G0));
     428          23 :     PetscCall(MatCopy(G0,B,SAME_NONZERO_PATTERN));
     429          23 :     PetscCall(MatDenseRestoreSubMatrix(d->G,&G0));
     430          23 :     PetscCall(DSRestoreMat(d->eps->ds,DS_MAT_B,&B));
     431             :   }
     432             :   /* Set the signature on projected matrix B */
     433          26 :   if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
     434           0 :     PetscCall(DSGetLeadingDimension(d->eps->ds,&ld));
     435           0 :     PetscCall(DSGetArray(d->eps->ds,DS_MAT_B,&pA));
     436           0 :     PetscCall(PetscArrayzero(pA,d->eps->nconv*ld));
     437           0 :     PetscCall(VecCreateSeq(PETSC_COMM_SELF,d->eps->nconv,&v));
     438           0 :     PetscCall(BVGetSignature(d->eps->V,v));
     439           0 :     PetscCall(VecGetArrayRead(v,&pv));
     440           0 :     for (i=0;i<d->eps->nconv;i++) pA[i+ld*i] = pv[i];
     441           0 :     PetscCall(VecRestoreArrayRead(v,&pv));
     442           0 :     PetscCall(VecDestroy(&v));
     443           0 :     PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_B,&pA));
     444             :   }
     445          26 :   PetscCall(DSSetState(d->eps->ds,DS_STATE_RAW));
     446          26 :   PetscCall(DSSolve(d->eps->ds,d->eps->eigr,d->eps->eigi));
     447          26 :   PetscCall(DSSynchronize(d->eps->ds,d->eps->eigr,d->eps->eigi));
     448          26 :   if (d->W) {
     449          80 :     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]));
     450             :   }
     451          26 :   PetscFunctionReturn(PETSC_SUCCESS);
     452             : }
     453             : 
     454             : /*
     455             :    Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
     456             :    the norm associated to the Schur pair, where i = r_s..r_e
     457             : */
     458        2453 : static PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d,PetscInt r_s,PetscInt r_e)
     459             : {
     460        2453 :   PetscInt       i,ldpX;
     461        2453 :   PetscScalar    *pX;
     462        2453 :   BV             BX = d->BX?d->BX:d->eps->V;
     463        2453 :   Vec            *R;
     464             : 
     465        2453 :   PetscFunctionBegin;
     466        2453 :   PetscCall(DSGetLeadingDimension(d->eps->ds,&ldpX));
     467        2453 :   PetscCall(DSGetArray(d->eps->ds,DS_MAT_Q,&pX));
     468             :   /* nX(i) <- ||X(i)|| */
     469        2453 :   PetscCall(dvd_improvex_compute_X(d,r_s,r_e,NULL,pX,ldpX));
     470        2453 :   PetscCall(SlepcVecPoolGetVecs(d->auxV,r_e-r_s,&R));
     471        4926 :   for (i=r_s;i<r_e;i++) {
     472             :     /* R(i-r_s) <- AV*pX(i) */
     473        2473 :     PetscCall(BVMultVec(d->AX,1.0,0.0,R[i-r_s],&pX[ldpX*i]));
     474             :     /* R(i-r_s) <- R(i-r_s) - eigr(i)*BX*pX(i) */
     475        2473 :     PetscCall(BVMultVec(BX,-d->eigr[i],1.0,R[i-r_s],&pX[ldpX*i]));
     476             :   }
     477        2453 :   PetscCall(DSRestoreArray(d->eps->ds,DS_MAT_Q,&pX));
     478        2453 :   PetscCall(d->calcpairs_proj_res(d,r_s,r_e,R));
     479        2453 :   PetscCall(SlepcVecPoolRestoreVecs(d->auxV,r_e-r_s,&R));
     480        2453 :   PetscFunctionReturn(PETSC_SUCCESS);
     481             : }
     482             : 
     483        9555 : static PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R)
     484             : {
     485        9555 :   PetscInt       i,l,k;
     486        9555 :   PetscBool      lindep=PETSC_FALSE;
     487        9555 :   BV             cX;
     488             : 
     489        9555 :   PetscFunctionBegin;
     490        9555 :   if (d->W) cX = d->W; /* If left subspace exists, R <- orth(cY, R), nR[i] <- ||R[i]|| */
     491        8229 :   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]|| */
     492             :   else cX = NULL; /* Otherwise, nR[i] <- ||R[i]|| */
     493             : 
     494        2040 :   if (cX) {
     495        2040 :     PetscCall(BVGetActiveColumns(cX,&l,&k));
     496        2040 :     PetscCall(BVSetActiveColumns(cX,0,l));
     497        4151 :     for (i=0;i<r_e-r_s;i++) PetscCall(BVOrthogonalizeVec(cX,R[i],NULL,&d->nR[r_s+i],&lindep));
     498        2040 :     PetscCall(BVSetActiveColumns(cX,l,k));
     499        2040 :     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]));
     500             :   } else {
     501       15030 :     for (i=0;i<r_e-r_s;i++) PetscCall(VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]));
     502       15030 :     for (i=0;i<r_e-r_s;i++) PetscCall(VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]));
     503             :   }
     504        9555 :   PetscFunctionReturn(PETSC_SUCCESS);
     505             : }
     506             : 
     507         285 : PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d,dvdBlackboard *b,PetscBool borth,PetscBool harm)
     508             : {
     509         285 :   PetscBool      std_probl,her_probl,ind_probl;
     510         285 :   DSType         dstype;
     511         285 :   Vec            v1;
     512             : 
     513         285 :   PetscFunctionBegin;
     514         285 :   std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
     515         285 :   her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
     516         285 :   ind_probl = DVD_IS(d->sEP,DVD_EP_INDEFINITE)? PETSC_TRUE: PETSC_FALSE;
     517             : 
     518             :   /* Setting configuration constrains */
     519         285 :   b->max_size_proj = PetscMax(b->max_size_proj,b->max_size_V);
     520         285 :   d->W_shift = d->B? PETSC_TRUE: PETSC_FALSE;
     521             : 
     522             :   /* Setup the step */
     523         285 :   if (b->state >= DVD_STATE_CONF) {
     524          95 :     d->max_size_P = b->max_size_P;
     525          95 :     d->max_size_proj = b->max_size_proj;
     526             :     /* Create a DS if the method works with Schur decompositions */
     527          95 :     d->calcPairs = dvd_calcpairs_proj;
     528          95 :     d->calcpairs_residual = dvd_calcpairs_res_0;
     529          95 :     d->calcpairs_proj_res = dvd_calcpairs_proj_res;
     530          95 :     d->calcpairs_selectPairs = dvd_calcpairs_selectPairs;
     531             :     /* Create and configure a DS for solving the projected problems */
     532          95 :     if (d->W) dstype = DSGNHEP;    /* If we use harmonics */
     533             :     else {
     534          95 :       if (ind_probl) dstype = DSGHIEP;
     535          95 :       else if (std_probl) dstype = her_probl? DSHEP : DSNHEP;
     536          23 :       else dstype = her_probl? DSGHEP : DSGNHEP;
     537             :     }
     538          95 :     PetscCall(DSSetType(d->eps->ds,dstype));
     539          95 :     PetscCall(DSAllocate(d->eps->ds,d->eps->ncv));
     540             :     /* Create various vector basis */
     541          95 :     if (harm) {
     542          23 :       PetscCall(BVDuplicateResize(d->eps->V,d->eps->ncv,&d->W));
     543          23 :       PetscCall(BVSetMatrix(d->W,NULL,PETSC_FALSE));
     544          72 :     } else d->W = NULL;
     545          95 :     PetscCall(BVDuplicateResize(d->eps->V,d->eps->ncv,&d->AX));
     546          95 :     PetscCall(BVSetMatrix(d->AX,NULL,PETSC_FALSE));
     547          95 :     PetscCall(BVDuplicateResize(d->eps->V,d->eps->ncv,&d->auxBV));
     548          95 :     PetscCall(BVSetMatrix(d->auxBV,NULL,PETSC_FALSE));
     549          95 :     if (d->B) {
     550          14 :       PetscCall(BVDuplicateResize(d->eps->V,d->eps->ncv,&d->BX));
     551          14 :       PetscCall(BVSetMatrix(d->BX,NULL,PETSC_FALSE));
     552          81 :     } else d->BX = NULL;
     553          95 :     PetscCall(MatCreateVecsEmpty(d->A,&v1,NULL));
     554          95 :     PetscCall(SlepcVecPoolCreate(v1,0,&d->auxV));
     555          95 :     PetscCall(VecDestroy(&v1));
     556             :     /* Create projected problem matrices */
     557          95 :     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->H));
     558          95 :     if (!std_probl) PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->G));
     559          72 :     else d->G = NULL;
     560          95 :     if (her_probl) {
     561          69 :       PetscCall(MatSetOption(d->H,MAT_HERMITIAN,PETSC_TRUE));
     562          69 :       if (d->G) PetscCall(MatSetOption(d->G,MAT_HERMITIAN,PETSC_TRUE));
     563             :     }
     564             : 
     565          95 :     if (ind_probl) PetscCall(PetscMalloc1(d->eps->ncv,&d->nBds));
     566          95 :     else d->nBds = NULL;
     567          95 :     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->auxM));
     568             : 
     569          95 :     PetscCall(EPSDavidsonFLAdd(&d->startList,dvd_calcpairs_qz_start));
     570          95 :     PetscCall(EPSDavidsonFLAdd(&d->endList,EPSXDComputeDSConv));
     571          95 :     PetscCall(EPSDavidsonFLAdd(&d->destroyList,dvd_calcpairs_qz_d));
     572             :   }
     573         285 :   PetscFunctionReturn(PETSC_SUCCESS);
     574             : }