LCOV - code coverage report
Current view: top level - pep/impls/krylov/toar - nrefine.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 1045 1079 96.8 %
Date: 2024-03-28 00:52:16 Functions: 17 17 100.0 %
Legend: Lines: hit not hit

          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             :    Newton refinement for polynomial eigenproblems.
      12             : 
      13             :    References:
      14             : 
      15             :        [1] T. Betcke and D. Kressner, "Perturbation, extraction and refinement
      16             :            of invariant pairs for matrix polynomials", Linear Algebra Appl.
      17             :            435(3):514-536, 2011.
      18             : 
      19             :        [2] C. Campos and J.E. Roman, "Parallel iterative refinement in
      20             :            polynomial eigenvalue problems", Numer. Linear Algebra Appl. 23(4):
      21             :            730-745, 2016.
      22             : */
      23             : 
      24             : #include <slepc/private/pepimpl.h>
      25             : #include <slepcblaslapack.h>
      26             : 
      27             : typedef struct {
      28             :   Mat          *A,M1;
      29             :   BV           V,M2,M3,W;
      30             :   PetscInt     k,nmat;
      31             :   PetscScalar  *fih,*work,*M4;
      32             :   PetscBLASInt *pM4;
      33             :   PetscBool    compM1;
      34             :   Vec          t;
      35             : } PEP_REFINE_MATSHELL;
      36             : 
      37             : typedef struct {
      38             :   Mat          E[2],M1;
      39             :   Vec          tN,ttN,t1,vseq;
      40             :   VecScatter   scatterctx;
      41             :   PetscBool    compM1;
      42             :   PetscInt     *map0,*map1,*idxg,*idxp;
      43             :   PetscSubcomm subc;
      44             :   VecScatter   scatter_sub;
      45             :   VecScatter   *scatter_id,*scatterp_id;
      46             :   Mat          *A;
      47             :   BV           V,W,M2,M3,Wt;
      48             :   PetscScalar  *M4,*w,*wt,*d,*dt;
      49             :   Vec          t,tg,Rv,Vi,tp,tpg;
      50             :   PetscInt     idx,*cols;
      51             : } PEP_REFINE_EXPLICIT;
      52             : 
      53        1616 : static PetscErrorCode MatMult_FS(Mat M ,Vec x,Vec y)
      54             : {
      55        1616 :   PEP_REFINE_MATSHELL *ctx;
      56        1616 :   PetscInt            k,i;
      57        1616 :   PetscScalar         *c;
      58        1616 :   PetscBLASInt        k_,one=1,info;
      59             : 
      60        1616 :   PetscFunctionBegin;
      61        1616 :   PetscCall(MatShellGetContext(M,&ctx));
      62        1616 :   PetscCall(VecCopy(x,ctx->t));
      63        1616 :   k    = ctx->k;
      64        1616 :   c    = ctx->work;
      65        1616 :   PetscCall(PetscBLASIntCast(k,&k_));
      66        1616 :   PetscCall(MatMult(ctx->M1,x,y));
      67        1616 :   PetscCall(VecConjugate(ctx->t));
      68        1616 :   PetscCall(BVDotVec(ctx->M3,ctx->t,c));
      69       21786 :   for (i=0;i<k;i++) c[i] = PetscConj(c[i]);
      70        1616 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
      71        1616 :   PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&k_,&one,ctx->M4,&k_,ctx->pM4,c,&k_,&info));
      72        1616 :   PetscCall(PetscFPTrapPop());
      73        1616 :   SlepcCheckLapackInfo("getrs",info);
      74        1616 :   PetscCall(BVMultVec(ctx->M2,-1.0,1.0,y,c));
      75        1616 :   PetscFunctionReturn(PETSC_SUCCESS);
      76             : }
      77             : 
      78             : /*
      79             :   Evaluates the first d elements of the polynomial basis
      80             :   on a given matrix H which is considered to be triangular
      81             : */
      82          42 : static PetscErrorCode PEPEvaluateBasisforMatrix(PEP pep,PetscInt nm,PetscInt k,PetscScalar *H,PetscInt ldh,PetscScalar *fH)
      83             : {
      84          42 :   PetscInt       i,j,ldfh=nm*k,off,nmat=pep->nmat;
      85          42 :   PetscReal      *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat,t;
      86          42 :   PetscScalar    corr=0.0,alpha,beta;
      87          42 :   PetscBLASInt   k_,ldh_,ldfh_;
      88             : 
      89          42 :   PetscFunctionBegin;
      90          42 :   PetscCall(PetscBLASIntCast(ldh,&ldh_));
      91          42 :   PetscCall(PetscBLASIntCast(k,&k_));
      92          42 :   PetscCall(PetscBLASIntCast(ldfh,&ldfh_));
      93          42 :   PetscCall(PetscArrayzero(fH,nm*k*k));
      94         381 :   if (nm>0) for (j=0;j<k;j++) fH[j+j*ldfh] = 1.0;
      95          42 :   if (nm>1) {
      96          42 :     t = b[0]/a[0];
      97          42 :     off = k;
      98         381 :     for (j=0;j<k;j++) {
      99        3708 :       for (i=0;i<k;i++) fH[off+i+j*ldfh] = H[i+j*ldh]/a[0];
     100         339 :       fH[j+j*ldfh] -= t;
     101             :     }
     102             :   }
     103          84 :   for (i=2;i<nm;i++) {
     104          42 :     off = i*k;
     105          42 :     if (i==2) {
     106         381 :       for (j=0;j<k;j++) {
     107         339 :         fH[off+j+j*ldfh] = 1.0;
     108         339 :         H[j+j*ldh] -= b[1];
     109             :       }
     110             :     } else {
     111           0 :       for (j=0;j<k;j++) {
     112           0 :         PetscCall(PetscArraycpy(fH+off+j*ldfh,fH+(i-2)*k+j*ldfh,k));
     113           0 :         H[j+j*ldh] += corr-b[i-1];
     114             :       }
     115             :     }
     116          42 :     corr  = b[i-1];
     117          42 :     beta  = -g[i-1]/a[i-1];
     118          42 :     alpha = 1/a[i-1];
     119          42 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&alpha,H,&ldh_,fH+(i-1)*k,&ldfh_,&beta,fH+off,&ldfh_));
     120             :   }
     121         381 :   for (j=0;j<k;j++) H[j+j*ldh] += corr;
     122          42 :   PetscFunctionReturn(PETSC_SUCCESS);
     123             : }
     124             : 
     125          67 : static PetscErrorCode NRefSysSetup_shell(PEP pep,PetscInt k,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar h,PEP_REFINE_MATSHELL *ctx)
     126             : {
     127          67 :   PetscScalar       *DHii,*T12,*Tr,*Ts,*array,s,ss,sone=1.0,zero=0.0,*M4=ctx->M4,t,*v,*T;
     128          67 :   const PetscScalar *m3,*m2;
     129          67 :   PetscInt          i,d,j,nmat=pep->nmat,lda=nmat*k,deg=nmat-1,nloc,ld2,ld3;
     130          67 :   PetscReal         *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
     131          67 :   PetscBLASInt      k_,lda_,lds_,nloc_,ld2_,one=1,info;
     132          67 :   Mat               *A=ctx->A,Mk,M1=ctx->M1,P;
     133          67 :   BV                V=ctx->V,M2=ctx->M2,M3=ctx->M3,W=ctx->W;
     134          67 :   MatStructure      str;
     135          67 :   Vec               vc;
     136             : 
     137          67 :   PetscFunctionBegin;
     138          67 :   PetscCall(STGetMatStructure(pep->st,&str));
     139          67 :   PetscCall(PetscMalloc3(nmat*k*k,&T12,k*k,&Tr,PetscMax(k*k,nmat),&Ts));
     140          67 :   DHii = T12;
     141          67 :   PetscCall(PetscArrayzero(DHii,k*k*nmat));
     142         916 :   for (i=0;i<k;i++) DHii[k+i+i*lda] = 1.0/a[0];
     143         134 :   for (d=2;d<nmat;d++) {
     144         916 :     for (j=0;j<k;j++) {
     145       12232 :       for (i=0;i<k;i++) {
     146       11383 :         DHii[d*k+i+j*lda] = ((h-b[d-1])*DHii[(d-1)*k+i+j*lda]+fH[(d-1)*k+i+j*lda]-g[d-1]*DHii[(d-2)*k+i+j*lda])/(a[d-1]);
     147             :       }
     148             :     }
     149             :   }
     150             :   /* T11 */
     151          67 :   if (!ctx->compM1) {
     152          61 :     PetscCall(MatCopy(A[0],M1,DIFFERENT_NONZERO_PATTERN));
     153          61 :     PetscCall(PEPEvaluateBasis(pep,h,0,Ts,NULL));
     154         183 :     for (j=1;j<nmat;j++) PetscCall(MatAXPY(M1,Ts[j],A[j],str));
     155             :   }
     156             : 
     157             :   /* T22 */
     158          67 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     159          67 :   PetscCall(PetscBLASIntCast(k,&k_));
     160          67 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     161          67 :   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,S,&lds_,S,&lds_,&zero,Tr,&k_));
     162         134 :   for (i=1;i<deg;i++) {
     163          67 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,Tr,&k_,DHii+i*k,&lda_,&zero,Ts,&k_));
     164          67 :     s = (i==1)?0.0:1.0;
     165          67 :     PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,fH+i*k,&lda_,Ts,&k_,&s,M4,&k_));
     166             :   }
     167        6183 :   for (i=0;i<k;i++) for (j=0;j<i;j++) { t=M4[i+j*k];M4[i+j*k]=M4[j+i*k];M4[j+i*k]=t; }
     168             : 
     169             :   /* T12 */
     170          67 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&Mk));
     171         201 :   for (i=1;i<nmat;i++) {
     172         134 :     PetscCall(MatDenseGetArrayWrite(Mk,&array));
     173         134 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,DHii+i*k,&lda_,&zero,array,&k_));
     174         134 :     PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
     175         134 :     PetscCall(BVSetActiveColumns(W,0,k));
     176         134 :     PetscCall(BVMult(W,1.0,0.0,V,Mk));
     177         134 :     if (i==1) PetscCall(BVMatMult(W,A[i],M2));
     178             :     else {
     179          67 :       PetscCall(BVMatMult(W,A[i],M3)); /* using M3 as work space */
     180         134 :       PetscCall(BVMult(M2,1.0,1.0,M3,NULL));
     181             :     }
     182             :   }
     183             : 
     184             :   /* T21 */
     185          67 :   PetscCall(MatDenseGetArrayWrite(Mk,&array));
     186         134 :   for (i=1;i<deg;i++) {
     187          67 :     s = (i==1)?0.0:1.0;
     188          67 :     ss = PetscConj(fh[i]);
     189          67 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&ss,S,&lds_,fH+i*k,&lda_,&s,array,&k_));
     190             :   }
     191          67 :   PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
     192          67 :   PetscCall(BVSetActiveColumns(M3,0,k));
     193          67 :   PetscCall(BVMult(M3,1.0,0.0,V,Mk));
     194         916 :   for (i=0;i<k;i++) {
     195         849 :     PetscCall(BVGetColumn(M3,i,&vc));
     196         849 :     PetscCall(VecConjugate(vc));
     197         849 :     PetscCall(BVRestoreColumn(M3,i,&vc));
     198             :   }
     199          67 :   PetscCall(MatDestroy(&Mk));
     200          67 :   PetscCall(PetscFree3(T12,Tr,Ts));
     201             : 
     202          67 :   PetscCall(VecGetLocalSize(ctx->t,&nloc));
     203          67 :   PetscCall(PetscBLASIntCast(nloc,&nloc_));
     204          67 :   PetscCall(PetscMalloc1(nloc*k,&T));
     205          67 :   PetscCall(KSPGetOperators(pep->refineksp,NULL,&P));
     206          67 :   if (!ctx->compM1) PetscCall(MatCopy(ctx->M1,P,SAME_NONZERO_PATTERN));
     207          67 :   PetscCall(BVGetArrayRead(ctx->M2,&m2));
     208          67 :   PetscCall(BVGetLeadingDimension(ctx->M2,&ld2));
     209          67 :   PetscCall(PetscBLASIntCast(ld2,&ld2_));
     210          67 :   PetscCall(BVGetArrayRead(ctx->M3,&m3));
     211          67 :   PetscCall(BVGetLeadingDimension(ctx->M3,&ld3));
     212          67 :   PetscCall(VecGetArray(ctx->t,&v));
     213       27547 :   for (i=0;i<nloc;i++) for (j=0;j<k;j++) T[j+i*k] = m3[i+j*ld3];
     214          67 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     215          67 :   PetscCallBLAS("LAPACKgesv",LAPACKgesv_(&k_,&nloc_,ctx->M4,&k_,ctx->pM4,T,&k_,&info));
     216          67 :   PetscCall(PetscFPTrapPop());
     217          67 :   SlepcCheckLapackInfo("gesv",info);
     218        2077 :   for (i=0;i<nloc;i++) v[i] = BLASdot_(&k_,m2+i,&ld2_,T+i*k,&one);
     219          67 :   PetscCall(VecRestoreArray(ctx->t,&v));
     220          67 :   PetscCall(BVRestoreArrayRead(ctx->M2,&m2));
     221          67 :   PetscCall(BVRestoreArrayRead(ctx->M3,&m3));
     222          67 :   PetscCall(MatDiagonalSet(P,ctx->t,ADD_VALUES));
     223          67 :   PetscCall(PetscFree(T));
     224          67 :   PetscCall(KSPSetUp(pep->refineksp));
     225          67 :   PetscFunctionReturn(PETSC_SUCCESS);
     226             : }
     227             : 
     228          67 : static PetscErrorCode NRefSysSolve_shell(KSP ksp,PetscInt nmat,Vec Rv,PetscScalar *Rh,PetscInt k,Vec dVi,PetscScalar *dHi)
     229             : {
     230          67 :   PetscScalar         *t0;
     231          67 :   PetscBLASInt        k_,one=1,info,lda_;
     232          67 :   PetscInt            i,lda=nmat*k;
     233          67 :   Mat                 M;
     234          67 :   PEP_REFINE_MATSHELL *ctx;
     235             : 
     236          67 :   PetscFunctionBegin;
     237          67 :   PetscCall(KSPGetOperators(ksp,&M,NULL));
     238          67 :   PetscCall(MatShellGetContext(M,&ctx));
     239          67 :   PetscCall(PetscCalloc1(k,&t0));
     240          67 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     241          67 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     242          67 :   PetscCall(PetscBLASIntCast(k,&k_));
     243         916 :   for (i=0;i<k;i++) t0[i] = Rh[i];
     244          67 :   PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&k_,&one,ctx->M4,&k_,ctx->pM4,t0,&k_,&info));
     245          67 :   SlepcCheckLapackInfo("getrs",info);
     246          67 :   PetscCall(BVMultVec(ctx->M2,-1.0,1.0,Rv,t0));
     247          67 :   PetscCall(KSPSolve(ksp,Rv,dVi));
     248          67 :   PetscCall(VecConjugate(dVi));
     249          67 :   PetscCall(BVDotVec(ctx->M3,dVi,dHi));
     250          67 :   PetscCall(VecConjugate(dVi));
     251         916 :   for (i=0;i<k;i++) dHi[i] = Rh[i]-PetscConj(dHi[i]);
     252          67 :   PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&k_,&one,ctx->M4,&k_,ctx->pM4,dHi,&k_,&info));
     253          67 :   SlepcCheckLapackInfo("getrs",info);
     254          67 :   PetscCall(PetscFPTrapPop());
     255          67 :   PetscCall(PetscFree(t0));
     256          67 :   PetscFunctionReturn(PETSC_SUCCESS);
     257             : }
     258             : 
     259             : /*
     260             :    Computes the residual P(H,V*S)*e_j for the polynomial
     261             : */
     262         113 : static PetscErrorCode NRefRightSide(PetscInt nmat,PetscReal *pcf,Mat *A,PetscInt k,BV V,PetscScalar *S,PetscInt lds,PetscInt j,PetscScalar *H,PetscInt ldh,PetscScalar *fH,PetscScalar *DfH,PetscScalar *dH,BV dV,PetscScalar *dVS,PetscInt rds,Vec Rv,PetscScalar *Rh,BV W,Vec t)
     263             : {
     264         113 :   PetscScalar    *DS0,*DS1,*F,beta=0.0,sone=1.0,none=-1.0,tt=0.0,*h,zero=0.0,*Z,*c0;
     265         113 :   PetscReal      *a=pcf,*b=pcf+nmat,*g=b+nmat;
     266         113 :   PetscInt       i,ii,jj,lda;
     267         113 :   PetscBLASInt   lda_,k_,ldh_,lds_,nmat_,k2_,krds_,j_,one=1;
     268         113 :   Mat            M0;
     269         113 :   Vec            w;
     270             : 
     271         113 :   PetscFunctionBegin;
     272         113 :   PetscCall(PetscMalloc4(k*nmat,&h,k*k,&DS0,k*k,&DS1,k*k,&Z));
     273         113 :   lda = k*nmat;
     274         113 :   PetscCall(PetscBLASIntCast(k,&k_));
     275         113 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     276         113 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     277         113 :   PetscCall(PetscBLASIntCast(nmat,&nmat_));
     278         113 :   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&nmat_,&k_,&sone,S,&lds_,fH+j*lda,&k_,&zero,h,&k_));
     279         113 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,nmat,h,&M0));
     280         113 :   PetscCall(BVSetActiveColumns(W,0,nmat));
     281         113 :   PetscCall(BVMult(W,1.0,0.0,V,M0));
     282         113 :   PetscCall(MatDestroy(&M0));
     283             : 
     284         113 :   PetscCall(BVGetColumn(W,0,&w));
     285         113 :   PetscCall(MatMult(A[0],w,Rv));
     286         113 :   PetscCall(BVRestoreColumn(W,0,&w));
     287         339 :   for (i=1;i<nmat;i++) {
     288         226 :     PetscCall(BVGetColumn(W,i,&w));
     289         226 :     PetscCall(MatMult(A[i],w,t));
     290         226 :     PetscCall(BVRestoreColumn(W,i,&w));
     291         226 :     PetscCall(VecAXPY(Rv,1.0,t));
     292             :   }
     293             :   /* Update right-hand side */
     294         113 :   if (j) {
     295          99 :     PetscCall(PetscBLASIntCast(ldh,&ldh_));
     296          99 :     PetscCall(PetscArrayzero(Z,k*k));
     297          99 :     PetscCall(PetscArrayzero(DS0,k*k));
     298          99 :     PetscCall(PetscArraycpy(Z+(j-1)*k,dH+(j-1)*k,k));
     299             :     /* Update DfH */
     300         297 :     for (i=1;i<nmat;i++) {
     301         198 :       if (i>1) {
     302          99 :         beta = -g[i-1];
     303          99 :         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,fH+(i-1)*k,&lda_,Z,&k_,&beta,DS0,&k_));
     304          99 :         tt += -b[i-1];
     305        1109 :         for (ii=0;ii<k;ii++) H[ii+ii*ldh] += tt;
     306          99 :         tt = b[i-1];
     307          99 :         beta = 1.0/a[i-1];
     308          99 :         PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&beta,DS1,&k_,H,&ldh_,&beta,DS0,&k_));
     309          99 :         F = DS0; DS0 = DS1; DS1 = F;
     310             :       } else {
     311          99 :         PetscCall(PetscArrayzero(DS1,k*k));
     312        1109 :         for (ii=0;ii<k;ii++) DS1[ii+(j-1)*k] = Z[ii+(j-1)*k]/a[0];
     313             :       }
     314        1208 :       for (jj=j;jj<k;jj++) {
     315       12948 :         for (ii=0;ii<k;ii++) DfH[k*i+ii+jj*lda] += DS1[ii+jj*k];
     316             :       }
     317             :     }
     318        1109 :     for (ii=0;ii<k;ii++) H[ii+ii*ldh] += tt;
     319             :     /* Update right-hand side */
     320          99 :     PetscCall(PetscBLASIntCast(2*k,&k2_));
     321          99 :     PetscCall(PetscBLASIntCast(j,&j_));
     322          99 :     PetscCall(PetscBLASIntCast(k+rds,&krds_));
     323          99 :     c0 = DS0;
     324          99 :     PetscCall(PetscArrayzero(Rh,k));
     325         396 :     for (i=0;i<nmat;i++) {
     326         297 :       PetscCallBLAS("BLASgemv",BLASgemv_("N",&krds_,&j_,&sone,dVS,&k2_,fH+j*lda+i*k,&one,&zero,h,&one));
     327         297 :       PetscCallBLAS("BLASgemv",BLASgemv_("N",&k_,&k_,&sone,S,&lds_,DfH+i*k+j*lda,&one,&sone,h,&one));
     328         297 :       PetscCall(BVMultVec(V,1.0,0.0,t,h));
     329         297 :       PetscCall(BVSetActiveColumns(dV,0,rds));
     330         297 :       PetscCall(BVMultVec(dV,1.0,1.0,t,h+k));
     331         297 :       PetscCall(BVGetColumn(W,i,&w));
     332         297 :       PetscCall(MatMult(A[i],t,w));
     333         297 :       PetscCall(BVRestoreColumn(W,i,&w));
     334         297 :       if (i>0 && i<nmat-1) {
     335          99 :         PetscCallBLAS("BLASgemv",BLASgemv_("C",&k_,&k_,&sone,S,&lds_,h,&one,&zero,c0,&one));
     336         297 :         PetscCallBLAS("BLASgemv",BLASgemv_("C",&k_,&k_,&none,fH+i*k,&lda_,c0,&one,&sone,Rh,&one));
     337             :       }
     338             :     }
     339             : 
     340         396 :     for (i=0;i<nmat;i++) h[i] = -1.0;
     341          99 :     PetscCall(BVMultVec(W,1.0,1.0,Rv,h));
     342             :   }
     343         113 :   PetscCall(PetscFree4(h,DS0,DS1,Z));
     344         113 :   PetscFunctionReturn(PETSC_SUCCESS);
     345             : }
     346             : 
     347         219 : static PetscErrorCode NRefSysSolve_mbe(PetscInt k,PetscInt sz,BV W,PetscScalar *w,BV Wt,PetscScalar *wt,PetscScalar *d,PetscScalar *dt,KSP ksp,BV T2,BV T3 ,PetscScalar *T4,PetscBool trans,Vec x1,PetscScalar *x2,Vec sol1,PetscScalar *sol2,Vec vw)
     348             : {
     349         219 :   PetscInt       i,j,incf,incc;
     350         219 :   PetscScalar    *y,*g,*xx2,*ww,y2,*dd;
     351         219 :   Vec            v,t,xx1;
     352         219 :   BV             WW,T;
     353             : 
     354         219 :   PetscFunctionBegin;
     355         219 :   PetscCall(PetscMalloc3(sz,&y,sz,&g,k,&xx2));
     356         219 :   if (trans) {
     357             :     WW = W; ww = w; dd = d; T = T3; incf = 0; incc = 1;
     358             :   } else {
     359         118 :     WW = Wt; ww = wt; dd = dt; T = T2; incf = 1; incc = 0;
     360             :   }
     361         219 :   xx1 = vw;
     362         219 :   PetscCall(VecCopy(x1,xx1));
     363         219 :   PetscCall(PetscArraycpy(xx2,x2,sz));
     364         219 :   PetscCall(PetscArrayzero(sol2,k));
     365         842 :   for (i=sz-1;i>=0;i--) {
     366         623 :     PetscCall(BVGetColumn(WW,i,&v));
     367         623 :     PetscCall(VecConjugate(v));
     368         623 :     PetscCall(VecDot(xx1,v,y+i));
     369         623 :     PetscCall(VecConjugate(v));
     370         623 :     PetscCall(BVRestoreColumn(WW,i,&v));
     371        1646 :     for (j=0;j<i;j++) y[i] += ww[j+i*k]*xx2[j];
     372         623 :     y[i] = -(y[i]-xx2[i])/dd[i];
     373         623 :     PetscCall(BVGetColumn(T,i,&t));
     374         623 :     PetscCall(VecAXPY(xx1,-y[i],t));
     375         623 :     PetscCall(BVRestoreColumn(T,i,&t));
     376        2269 :     for (j=0;j<=i;j++) xx2[j] -= y[i]*T4[j*incf+incc*i+(i*incf+incc*j)*k];
     377         623 :     g[i] = xx2[i];
     378             :   }
     379         219 :   if (trans) PetscCall(KSPSolveTranspose(ksp,xx1,sol1));
     380         118 :   else PetscCall(KSPSolve(ksp,xx1,sol1));
     381         219 :   if (trans) {
     382             :     WW = Wt; ww = wt; dd = dt; T = T2; incf = 1; incc = 0;
     383             :   } else {
     384         118 :     WW = W; ww = w; dd = d; T = T3; incf = 0; incc = 1;
     385             :   }
     386         842 :   for (i=0;i<sz;i++) {
     387         623 :     PetscCall(BVGetColumn(T,i,&t));
     388         623 :     PetscCall(VecConjugate(t));
     389         623 :     PetscCall(VecDot(sol1,t,&y2));
     390         623 :     PetscCall(VecConjugate(t));
     391         623 :     PetscCall(BVRestoreColumn(T,i,&t));
     392        1646 :     for (j=0;j<i;j++) y2 += sol2[j]*T4[j*incf+incc*i+(i*incf+incc*j)*k];
     393         623 :     y2 = (g[i]-y2)/dd[i];
     394         623 :     PetscCall(BVGetColumn(WW,i,&v));
     395         623 :     PetscCall(VecAXPY(sol1,-y2,v));
     396        1646 :     for (j=0;j<i;j++) sol2[j] -= ww[j+i*k]*y2;
     397         623 :     sol2[i] = y[i]+y2;
     398         623 :     PetscCall(BVRestoreColumn(WW,i,&v));
     399             :   }
     400         219 :   PetscCall(PetscFree3(y,g,xx2));
     401         219 :   PetscFunctionReturn(PETSC_SUCCESS);
     402             : }
     403             : 
     404          17 : static PetscErrorCode NRefSysSetup_mbe(PEP pep,PetscInt k,KSP ksp,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar h,BV V,PEP_REFINE_EXPLICIT *matctx)
     405             : {
     406          17 :   PetscInt       i,j,l,nmat=pep->nmat,lda=nmat*k,deg=nmat-1;
     407          17 :   Mat            M1=matctx->M1,*A,*At,Mk;
     408          17 :   PetscReal      *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
     409          17 :   PetscScalar    s,ss,*DHii,*T12,*array,*Ts,*Tr,*M4=matctx->M4,sone=1.0,zero=0.0;
     410          17 :   PetscScalar    *w=matctx->w,*wt=matctx->wt,*d=matctx->d,*dt=matctx->dt;
     411          17 :   PetscBLASInt   lds_,lda_,k_;
     412          17 :   MatStructure   str;
     413          17 :   PetscBool      flg;
     414          17 :   BV             M2=matctx->M2,M3=matctx->M3,W=matctx->W,Wt=matctx->Wt;
     415          17 :   Vec            vc,vc2;
     416             : 
     417          17 :   PetscFunctionBegin;
     418          17 :   PetscCall(PetscMalloc3(nmat*k*k,&T12,k*k,&Tr,PetscMax(k*k,nmat),&Ts));
     419          17 :   PetscCall(STGetMatStructure(pep->st,&str));
     420          17 :   PetscCall(STGetTransform(pep->st,&flg));
     421          17 :   if (flg) {
     422           7 :     PetscCall(PetscMalloc1(pep->nmat,&At));
     423          28 :     for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
     424          10 :   } else At = pep->A;
     425          17 :   if (matctx->subc) A = matctx->A;
     426          11 :   else A = At;
     427             :   /* Form the explicit system matrix */
     428          17 :   DHii = T12;
     429          17 :   PetscCall(PetscArrayzero(DHii,k*k*nmat));
     430         118 :   for (i=0;i<k;i++) DHii[k+i+i*lda] = 1.0/a[0];
     431          34 :   for (l=2;l<nmat;l++) {
     432         118 :     for (j=0;j<k;j++) {
     433         724 :       for (i=0;i<k;i++) {
     434         623 :         DHii[l*k+i+j*lda] = ((h-b[l-1])*DHii[(l-1)*k+i+j*lda]+fH[(l-1)*k+i+j*lda]-g[l-1]*DHii[(l-2)*k+i+j*lda])/a[l-1];
     435             :       }
     436             :     }
     437             :   }
     438             : 
     439             :   /* T11 */
     440          17 :   if (!matctx->compM1) {
     441          13 :     PetscCall(MatCopy(A[0],M1,DIFFERENT_NONZERO_PATTERN));
     442          13 :     PetscCall(PEPEvaluateBasis(pep,h,0,Ts,NULL));
     443          39 :     for (j=1;j<nmat;j++) PetscCall(MatAXPY(M1,Ts[j],A[j],str));
     444             :   }
     445             : 
     446             :   /* T22 */
     447          17 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     448          17 :   PetscCall(PetscBLASIntCast(k,&k_));
     449          17 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     450          17 :   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,S,&lds_,S,&lds_,&zero,Tr,&k_));
     451          34 :   for (i=1;i<deg;i++) {
     452          17 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,Tr,&k_,DHii+i*k,&lda_,&zero,Ts,&k_));
     453          17 :     s = (i==1)?0.0:1.0;
     454          17 :     PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,fH+i*k,&lda_,Ts,&k_,&s,M4,&k_));
     455             :   }
     456             : 
     457             :   /* T12 */
     458          17 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&Mk));
     459          51 :   for (i=1;i<nmat;i++) {
     460          34 :     PetscCall(MatDenseGetArrayWrite(Mk,&array));
     461          34 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,DHii+i*k,&lda_,&zero,array,&k_));
     462          34 :     PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
     463          34 :     PetscCall(BVSetActiveColumns(W,0,k));
     464          34 :     PetscCall(BVMult(W,1.0,0.0,V,Mk));
     465          34 :     if (i==1) PetscCall(BVMatMult(W,A[i],M2));
     466             :     else {
     467          17 :       PetscCall(BVMatMult(W,A[i],M3)); /* using M3 as work space */
     468          34 :       PetscCall(BVMult(M2,1.0,1.0,M3,NULL));
     469             :     }
     470             :   }
     471             : 
     472             :   /* T21 */
     473          17 :   PetscCall(MatDenseGetArrayWrite(Mk,&array));
     474          34 :   for (i=1;i<deg;i++) {
     475          17 :     s = (i==1)?0.0:1.0;
     476          17 :     ss = PetscConj(fh[i]);
     477          17 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&ss,S,&lds_,fH+i*k,&lda_,&s,array,&k_));
     478             :   }
     479          17 :   PetscCall(MatDenseRestoreArrayWrite(Mk,&array));
     480          17 :   PetscCall(BVSetActiveColumns(M3,0,k));
     481          17 :   PetscCall(BVMult(M3,1.0,0.0,V,Mk));
     482         118 :   for (i=0;i<k;i++) {
     483         101 :     PetscCall(BVGetColumn(M3,i,&vc));
     484         101 :     PetscCall(VecConjugate(vc));
     485         101 :     PetscCall(BVRestoreColumn(M3,i,&vc));
     486             :   }
     487             : 
     488          17 :   PetscCall(PEP_KSPSetOperators(ksp,M1,M1));
     489          17 :   PetscCall(KSPSetUp(ksp));
     490          17 :   PetscCall(MatDestroy(&Mk));
     491             : 
     492             :   /* Set up for BEMW */
     493         118 :   for (i=0;i<k;i++) {
     494         101 :     PetscCall(BVGetColumn(M2,i,&vc));
     495         101 :     PetscCall(BVGetColumn(W,i,&vc2));
     496         101 :     PetscCall(NRefSysSolve_mbe(k,i,W,w,Wt,wt,d,dt,ksp,M2,M3,M4,PETSC_FALSE,vc,M4+i*k,vc2,w+i*k,matctx->t));
     497         101 :     PetscCall(BVRestoreColumn(M2,i,&vc));
     498         101 :     PetscCall(BVGetColumn(M3,i,&vc));
     499         101 :     PetscCall(VecConjugate(vc));
     500         101 :     PetscCall(VecDot(vc2,vc,&d[i]));
     501         101 :     PetscCall(VecConjugate(vc));
     502         101 :     PetscCall(BVRestoreColumn(M3,i,&vc));
     503         362 :     for (j=0;j<i;j++) d[i] += M4[i+j*k]*w[j+i*k];
     504         101 :     d[i] = M4[i+i*k]-d[i];
     505         101 :     PetscCall(BVRestoreColumn(W,i,&vc2));
     506             : 
     507         101 :     PetscCall(BVGetColumn(M3,i,&vc));
     508         101 :     PetscCall(BVGetColumn(Wt,i,&vc2));
     509         463 :     for (j=0;j<=i;j++) Ts[j] = M4[i+j*k];
     510         101 :     PetscCall(NRefSysSolve_mbe(k,i,W,w,Wt,wt,d,dt,ksp,M2,M3,M4,PETSC_TRUE,vc,Ts,vc2,wt+i*k,matctx->t));
     511         101 :     PetscCall(BVRestoreColumn(M3,i,&vc));
     512         101 :     PetscCall(BVGetColumn(M2,i,&vc));
     513         101 :     PetscCall(VecConjugate(vc2));
     514         101 :     PetscCall(VecDot(vc,vc2,&dt[i]));
     515         101 :     PetscCall(VecConjugate(vc2));
     516         101 :     PetscCall(BVRestoreColumn(M2,i,&vc));
     517         362 :     for (j=0;j<i;j++) dt[i] += M4[j+i*k]*wt[j+i*k];
     518         101 :     dt[i] = M4[i+i*k]-dt[i];
     519         101 :     PetscCall(BVRestoreColumn(Wt,i,&vc2));
     520             :   }
     521             : 
     522          17 :   if (flg) PetscCall(PetscFree(At));
     523          17 :   PetscCall(PetscFree3(T12,Tr,Ts));
     524          17 :   PetscFunctionReturn(PETSC_SUCCESS);
     525             : }
     526             : 
     527          17 : static PetscErrorCode NRefSysSetup_explicit(PEP pep,PetscInt k,KSP ksp,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar h,BV V,PEP_REFINE_EXPLICIT *matctx,BV W)
     528             : {
     529          17 :   PetscInt          i,j,d,n,n0,m0,n1,m1,nmat=pep->nmat,lda=nmat*k,deg=nmat-1;
     530          17 :   PetscInt          *idxg=matctx->idxg,*idxp=matctx->idxp,idx,ncols;
     531          17 :   Mat               M,*E=matctx->E,*A,*At,Mk,Md;
     532          17 :   PetscReal         *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;
     533          17 :   PetscScalar       s,ss,*DHii,*T22,*T21,*T12,*Ts,*Tr,*array,*ts,sone=1.0,zero=0.0;
     534          17 :   PetscBLASInt      lds_,lda_,k_;
     535          17 :   const PetscInt    *idxmc;
     536          17 :   const PetscScalar *valsc,*carray;
     537          17 :   MatStructure      str;
     538          17 :   Vec               vc,vc0;
     539          17 :   PetscBool         flg;
     540             : 
     541          17 :   PetscFunctionBegin;
     542          17 :   PetscCall(PetscMalloc5(k*k,&T22,k*k,&T21,nmat*k*k,&T12,k*k,&Tr,k*k,&Ts));
     543          17 :   PetscCall(STGetMatStructure(pep->st,&str));
     544          17 :   PetscCall(KSPGetOperators(ksp,&M,NULL));
     545          17 :   PetscCall(MatGetOwnershipRange(E[1],&n1,&m1));
     546          17 :   PetscCall(MatGetOwnershipRange(E[0],&n0,&m0));
     547          17 :   PetscCall(MatGetOwnershipRange(M,&n,NULL));
     548          17 :   PetscCall(PetscMalloc1(nmat,&ts));
     549          17 :   PetscCall(STGetTransform(pep->st,&flg));
     550          17 :   if (flg) {
     551           7 :     PetscCall(PetscMalloc1(pep->nmat,&At));
     552          28 :     for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
     553          10 :   } else At = pep->A;
     554          17 :   if (matctx->subc) A = matctx->A;
     555          11 :   else A = At;
     556             :   /* Form the explicit system matrix */
     557          17 :   DHii = T12;
     558          17 :   PetscCall(PetscArrayzero(DHii,k*k*nmat));
     559         118 :   for (i=0;i<k;i++) DHii[k+i+i*lda] = 1.0/a[0];
     560          34 :   for (d=2;d<nmat;d++) {
     561         118 :     for (j=0;j<k;j++) {
     562         724 :       for (i=0;i<k;i++) {
     563         623 :         DHii[d*k+i+j*lda] = ((h-b[d-1])*DHii[(d-1)*k+i+j*lda]+fH[(d-1)*k+i+j*lda]-g[d-1]*DHii[(d-2)*k+i+j*lda])/a[d-1];
     564             :       }
     565             :     }
     566             :   }
     567             : 
     568             :   /* T11 */
     569          17 :   if (!matctx->compM1) {
     570          13 :     PetscCall(MatCopy(A[0],E[0],DIFFERENT_NONZERO_PATTERN));
     571          13 :     PetscCall(PEPEvaluateBasis(pep,h,0,Ts,NULL));
     572          39 :     for (j=1;j<nmat;j++) PetscCall(MatAXPY(E[0],Ts[j],A[j],str));
     573             :   }
     574         527 :   for (i=n0;i<m0;i++) {
     575         510 :     PetscCall(MatGetRow(E[0],i,&ncols,&idxmc,&valsc));
     576         510 :     idx = n+i-n0;
     577        2006 :     for (j=0;j<ncols;j++) {
     578        1496 :       idxg[j] = matctx->map0[idxmc[j]];
     579             :     }
     580         510 :     PetscCall(MatSetValues(M,1,&idx,ncols,idxg,valsc,INSERT_VALUES));
     581         510 :     PetscCall(MatRestoreRow(E[0],i,&ncols,&idxmc,&valsc));
     582             :   }
     583             : 
     584             :   /* T22 */
     585          17 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     586          17 :   PetscCall(PetscBLASIntCast(k,&k_));
     587          17 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     588          17 :   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,S,&lds_,S,&lds_,&zero,Tr,&k_));
     589          34 :   for (i=1;i<deg;i++) {
     590          17 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,Tr,&k_,DHii+i*k,&lda_,&zero,Ts,&k_));
     591          17 :     s = (i==1)?0.0:1.0;
     592          17 :     PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&k_,&k_,&k_,&sone,fH+i*k,&lda_,Ts,&k_,&s,T22,&k_));
     593             :   }
     594         118 :   for (j=0;j<k;j++) idxp[j] = matctx->map1[j];
     595         118 :   for (i=0;i<m1-n1;i++) {
     596         101 :     idx = n+m0-n0+i;
     597         724 :     for (j=0;j<k;j++) {
     598         623 :       Tr[j] = T22[n1+i+j*k];
     599             :     }
     600         101 :     PetscCall(MatSetValues(M,1,&idx,k,idxp,Tr,INSERT_VALUES));
     601             :   }
     602             : 
     603             :   /* T21 */
     604          34 :   for (i=1;i<deg;i++) {
     605          17 :     s = (i==1)?0.0:1.0;
     606          17 :     ss = PetscConj(fh[i]);
     607          17 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&ss,S,&lds_,fH+i*k,&lda_,&s,T21,&k_));
     608             :   }
     609          17 :   PetscCall(BVSetActiveColumns(W,0,k));
     610          17 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,T21,&Mk));
     611          17 :   PetscCall(BVMult(W,1.0,0.0,V,Mk));
     612         118 :   for (i=0;i<k;i++) {
     613         101 :     PetscCall(BVGetColumn(W,i,&vc));
     614         101 :     PetscCall(VecConjugate(vc));
     615         101 :     PetscCall(VecGetArrayRead(vc,&carray));
     616         101 :     idx = matctx->map1[i];
     617         101 :     PetscCall(MatSetValues(M,1,&idx,m0-n0,matctx->map0+n0,carray,INSERT_VALUES));
     618         101 :     PetscCall(VecRestoreArrayRead(vc,&carray));
     619         101 :     PetscCall(BVRestoreColumn(W,i,&vc));
     620             :   }
     621             : 
     622             :   /* T12 */
     623          51 :   for (i=1;i<nmat;i++) {
     624          34 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,DHii+i*k,&lda_,&zero,Ts,&k_));
     625         236 :     for (j=0;j<k;j++) PetscCall(PetscArraycpy(T12+i*k+j*lda,Ts+j*k,k));
     626             :   }
     627          17 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,nmat-1,NULL,&Md));
     628          68 :   for (i=0;i<nmat;i++) ts[i] = 1.0;
     629         118 :   for (j=0;j<k;j++) {
     630         101 :     PetscCall(MatDenseGetArrayWrite(Md,&array));
     631         101 :     PetscCall(PetscArraycpy(array,T12+k+j*lda,(nmat-1)*k));
     632         101 :     PetscCall(MatDenseRestoreArrayWrite(Md,&array));
     633         101 :     PetscCall(BVSetActiveColumns(W,0,nmat-1));
     634         101 :     PetscCall(BVMult(W,1.0,0.0,V,Md));
     635         303 :     for (i=nmat-1;i>0;i--) {
     636         202 :       PetscCall(BVGetColumn(W,i-1,&vc0));
     637         202 :       PetscCall(BVGetColumn(W,i,&vc));
     638         202 :       PetscCall(MatMult(A[i],vc0,vc));
     639         202 :       PetscCall(BVRestoreColumn(W,i-1,&vc0));
     640         202 :       PetscCall(BVRestoreColumn(W,i,&vc));
     641             :     }
     642         101 :     PetscCall(BVSetActiveColumns(W,1,nmat));
     643         101 :     PetscCall(BVGetColumn(W,0,&vc0));
     644         101 :     PetscCall(BVMultVec(W,1.0,0.0,vc0,ts));
     645         101 :     PetscCall(VecGetArrayRead(vc0,&carray));
     646         101 :     idx = matctx->map1[j];
     647         101 :     PetscCall(MatSetValues(M,m0-n0,matctx->map0+n0,1,&idx,carray,INSERT_VALUES));
     648         101 :     PetscCall(VecRestoreArrayRead(vc0,&carray));
     649         101 :     PetscCall(BVRestoreColumn(W,0,&vc0));
     650             :   }
     651          17 :   PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
     652          17 :   PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
     653          17 :   PetscCall(PEP_KSPSetOperators(ksp,M,M));
     654          17 :   PetscCall(KSPSetUp(ksp));
     655          17 :   PetscCall(PetscFree(ts));
     656          17 :   PetscCall(MatDestroy(&Mk));
     657          17 :   PetscCall(MatDestroy(&Md));
     658          17 :   if (flg) PetscCall(PetscFree(At));
     659          17 :   PetscCall(PetscFree5(T22,T21,T12,Tr,Ts));
     660          17 :   PetscFunctionReturn(PETSC_SUCCESS);
     661             : }
     662             : 
     663          17 : static PetscErrorCode NRefSysSolve_explicit(PetscInt k,KSP ksp,Vec Rv,PetscScalar *Rh,Vec dVi,PetscScalar *dHi,PEP_REFINE_EXPLICIT *matctx)
     664             : {
     665          17 :   PetscInt          n0,m0,n1,m1,i;
     666          17 :   PetscScalar       *arrayV;
     667          17 :   const PetscScalar *array;
     668             : 
     669          17 :   PetscFunctionBegin;
     670          17 :   PetscCall(MatGetOwnershipRange(matctx->E[1],&n1,&m1));
     671          17 :   PetscCall(MatGetOwnershipRange(matctx->E[0],&n0,&m0));
     672             : 
     673             :   /* Right side */
     674          17 :   PetscCall(VecGetArrayRead(Rv,&array));
     675          17 :   PetscCall(VecSetValues(matctx->tN,m0-n0,matctx->map0+n0,array,INSERT_VALUES));
     676          17 :   PetscCall(VecRestoreArrayRead(Rv,&array));
     677          17 :   PetscCall(VecSetValues(matctx->tN,m1-n1,matctx->map1+n1,Rh+n1,INSERT_VALUES));
     678          17 :   PetscCall(VecAssemblyBegin(matctx->tN));
     679          17 :   PetscCall(VecAssemblyEnd(matctx->tN));
     680             : 
     681             :   /* Solve */
     682          17 :   PetscCall(KSPSolve(ksp,matctx->tN,matctx->ttN));
     683             : 
     684             :   /* Retrieve solution */
     685          17 :   PetscCall(VecGetArray(dVi,&arrayV));
     686          17 :   PetscCall(VecGetArrayRead(matctx->ttN,&array));
     687          17 :   PetscCall(PetscArraycpy(arrayV,array,m0-n0));
     688          17 :   PetscCall(VecRestoreArray(dVi,&arrayV));
     689          17 :   if (!matctx->subc) {
     690          11 :     PetscCall(VecGetArray(matctx->t1,&arrayV));
     691          76 :     for (i=0;i<m1-n1;i++) arrayV[i] =  array[m0-n0+i];
     692          11 :     PetscCall(VecRestoreArray(matctx->t1,&arrayV));
     693          11 :     PetscCall(VecRestoreArrayRead(matctx->ttN,&array));
     694          11 :     PetscCall(VecScatterBegin(matctx->scatterctx,matctx->t1,matctx->vseq,INSERT_VALUES,SCATTER_FORWARD));
     695          11 :     PetscCall(VecScatterEnd(matctx->scatterctx,matctx->t1,matctx->vseq,INSERT_VALUES,SCATTER_FORWARD));
     696          11 :     PetscCall(VecGetArrayRead(matctx->vseq,&array));
     697          76 :     for (i=0;i<k;i++) dHi[i] = array[i];
     698          11 :     PetscCall(VecRestoreArrayRead(matctx->vseq,&array));
     699             :   }
     700          17 :   PetscFunctionReturn(PETSC_SUCCESS);
     701             : }
     702             : 
     703         113 : static PetscErrorCode NRefSysIter(PetscInt i,PEP pep,PetscInt k,KSP ksp,PetscScalar *fH,PetscScalar *S,PetscInt lds,PetscScalar *fh,PetscScalar *H,PetscInt ldh,Vec Rv,PetscScalar *Rh,BV V,Vec dVi,PetscScalar *dHi,PEP_REFINE_EXPLICIT *matctx,BV W)
     704             : {
     705         113 :   PetscInt            j,m,lda=pep->nmat*k,n0,m0,idx;
     706         113 :   PetscMPIInt         root,len;
     707         113 :   PetscScalar         *array2,h;
     708         113 :   const PetscScalar   *array;
     709         113 :   Vec                 R,Vi;
     710         113 :   PEP_REFINE_MATSHELL *ctx;
     711         113 :   Mat                 M;
     712             : 
     713         113 :   PetscFunctionBegin;
     714         113 :   if (!matctx || !matctx->subc) {
     715         356 :     for (j=0;j<pep->nmat;j++) fh[j] = fH[j*k+i+i*lda];
     716          89 :     h   = H[i+i*ldh];
     717          89 :     idx = i;
     718          89 :     R   = Rv;
     719          89 :     Vi  = dVi;
     720          89 :     switch (pep->scheme) {
     721          11 :     case PEP_REFINE_SCHEME_EXPLICIT:
     722          11 :       PetscCall(NRefSysSetup_explicit(pep,k,ksp,fH,S,lds,fh,h,V,matctx,W));
     723          11 :       matctx->compM1 = PETSC_FALSE;
     724          11 :       break;
     725          11 :     case PEP_REFINE_SCHEME_MBE:
     726          11 :       PetscCall(NRefSysSetup_mbe(pep,k,ksp,fH,S,lds,fh,h,V,matctx));
     727          11 :       matctx->compM1 = PETSC_FALSE;
     728          11 :       break;
     729          67 :     case PEP_REFINE_SCHEME_SCHUR:
     730          67 :       PetscCall(KSPGetOperators(ksp,&M,NULL));
     731          67 :       PetscCall(MatShellGetContext(M,&ctx));
     732          67 :       PetscCall(NRefSysSetup_shell(pep,k,fH,S,lds,fh,h,ctx));
     733          67 :       ctx->compM1 = PETSC_FALSE;
     734          67 :       break;
     735             :     }
     736             :   } else {
     737          24 :     if (i%matctx->subc->n==0 && (idx=i+matctx->subc->color)<k) {
     738          48 :       for (j=0;j<pep->nmat;j++) fh[j] = fH[j*k+idx+idx*lda];
     739          12 :       h = H[idx+idx*ldh];
     740          12 :       matctx->idx = idx;
     741          12 :       switch (pep->scheme) {
     742           6 :       case PEP_REFINE_SCHEME_EXPLICIT:
     743           6 :         PetscCall(NRefSysSetup_explicit(pep,k,ksp,fH,S,lds,fh,h,matctx->V,matctx,matctx->W));
     744           6 :         matctx->compM1 = PETSC_FALSE;
     745           6 :         break;
     746           6 :       case PEP_REFINE_SCHEME_MBE:
     747           6 :         PetscCall(NRefSysSetup_mbe(pep,k,ksp,fH,S,lds,fh,h,matctx->V,matctx));
     748           6 :         matctx->compM1 = PETSC_FALSE;
     749           6 :         break;
     750             :       case PEP_REFINE_SCHEME_SCHUR:
     751             :         break;
     752             :       }
     753          12 :     } else idx = matctx->idx;
     754          24 :     PetscCall(VecScatterBegin(matctx->scatter_id[i%matctx->subc->n],Rv,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
     755          24 :     PetscCall(VecScatterEnd(matctx->scatter_id[i%matctx->subc->n],Rv,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
     756          24 :     PetscCall(VecGetArrayRead(matctx->tg,&array));
     757          24 :     PetscCall(VecPlaceArray(matctx->t,array));
     758          24 :     PetscCall(VecCopy(matctx->t,matctx->Rv));
     759          24 :     PetscCall(VecResetArray(matctx->t));
     760          24 :     PetscCall(VecRestoreArrayRead(matctx->tg,&array));
     761          24 :     R  = matctx->Rv;
     762          24 :     Vi = matctx->Vi;
     763             :   }
     764         113 :   if (idx==i && idx<k) {
     765         101 :     switch (pep->scheme) {
     766          17 :       case PEP_REFINE_SCHEME_EXPLICIT:
     767          17 :         PetscCall(NRefSysSolve_explicit(k,ksp,R,Rh,Vi,dHi,matctx));
     768             :         break;
     769          17 :       case PEP_REFINE_SCHEME_MBE:
     770          17 :         PetscCall(NRefSysSolve_mbe(k,k,matctx->W,matctx->w,matctx->Wt,matctx->wt,matctx->d,matctx->dt,ksp,matctx->M2,matctx->M3 ,matctx->M4,PETSC_FALSE,R,Rh,Vi,dHi,matctx->t));
     771             :         break;
     772          67 :       case PEP_REFINE_SCHEME_SCHUR:
     773          67 :         PetscCall(NRefSysSolve_shell(ksp,pep->nmat,R,Rh,k,Vi,dHi));
     774             :         break;
     775             :     }
     776          12 :   }
     777         113 :   if (matctx && matctx->subc) {
     778          24 :     PetscCall(VecGetLocalSize(Vi,&m));
     779          24 :     PetscCall(VecGetArrayRead(Vi,&array));
     780          24 :     PetscCall(VecGetArray(matctx->tg,&array2));
     781          24 :     PetscCall(PetscArraycpy(array2,array,m));
     782          24 :     PetscCall(VecRestoreArray(matctx->tg,&array2));
     783          24 :     PetscCall(VecRestoreArrayRead(Vi,&array));
     784          24 :     PetscCall(VecScatterBegin(matctx->scatter_id[i%matctx->subc->n],matctx->tg,dVi,INSERT_VALUES,SCATTER_REVERSE));
     785          24 :     PetscCall(VecScatterEnd(matctx->scatter_id[i%matctx->subc->n],matctx->tg,dVi,INSERT_VALUES,SCATTER_REVERSE));
     786          24 :     switch (pep->scheme) {
     787          12 :     case PEP_REFINE_SCHEME_EXPLICIT:
     788          12 :       PetscCall(MatGetOwnershipRange(matctx->E[0],&n0,&m0));
     789          12 :       PetscCall(VecGetArrayRead(matctx->ttN,&array));
     790          12 :       PetscCall(VecPlaceArray(matctx->tp,array+m0-n0));
     791          12 :       PetscCall(VecScatterBegin(matctx->scatterp_id[i%matctx->subc->n],matctx->tp,matctx->tpg,INSERT_VALUES,SCATTER_FORWARD));
     792          12 :       PetscCall(VecScatterEnd(matctx->scatterp_id[i%matctx->subc->n],matctx->tp,matctx->tpg,INSERT_VALUES,SCATTER_FORWARD));
     793          12 :       PetscCall(VecResetArray(matctx->tp));
     794          12 :       PetscCall(VecRestoreArrayRead(matctx->ttN,&array));
     795          12 :       PetscCall(VecGetArrayRead(matctx->tpg,&array));
     796          84 :       for (j=0;j<k;j++) dHi[j] = array[j];
     797          12 :       PetscCall(VecRestoreArrayRead(matctx->tpg,&array));
     798             :       break;
     799             :      case PEP_REFINE_SCHEME_MBE:
     800             :       root = 0;
     801          18 :       for (j=0;j<i%matctx->subc->n;j++) root += matctx->subc->subsize[j];
     802          12 :       PetscCall(PetscMPIIntCast(k,&len));
     803          24 :       PetscCallMPI(MPI_Bcast(dHi,len,MPIU_SCALAR,root,matctx->subc->dupparent));
     804             :       break;
     805             :     case PEP_REFINE_SCHEME_SCHUR:
     806             :       break;
     807             :     }
     808          89 :   }
     809         113 :   PetscFunctionReturn(PETSC_SUCCESS);
     810             : }
     811             : 
     812          14 : static PetscErrorCode PEPNRefForwardSubstitution(PEP pep,PetscInt k,PetscScalar *S,PetscInt lds,PetscScalar *H,PetscInt ldh,PetscScalar *fH,BV dV,PetscScalar *dVS,PetscInt *rds,PetscScalar *dH,PetscInt lddh,KSP ksp,PEP_REFINE_EXPLICIT *matctx)
     813             : {
     814          14 :   PetscInt            i,nmat=pep->nmat,lda=nmat*k;
     815          14 :   PetscScalar         *fh,*Rh,*DfH;
     816          14 :   PetscReal           norm;
     817          14 :   BV                  W;
     818          14 :   Vec                 Rv,t,dvi;
     819          14 :   PEP_REFINE_MATSHELL *ctx;
     820          14 :   Mat                 M,*At;
     821          14 :   PetscBool           flg,lindep;
     822             : 
     823          14 :   PetscFunctionBegin;
     824          14 :   PetscCall(PetscMalloc2(nmat*k*k,&DfH,k,&Rh));
     825          14 :   *rds = 0;
     826          14 :   PetscCall(BVCreateVec(pep->V,&Rv));
     827          14 :   switch (pep->scheme) {
     828           4 :   case PEP_REFINE_SCHEME_EXPLICIT:
     829           4 :     PetscCall(BVCreateVec(pep->V,&t));
     830           4 :     PetscCall(BVDuplicateResize(pep->V,PetscMax(k,nmat),&W));
     831           4 :     PetscCall(PetscMalloc1(nmat,&fh));
     832             :     break;
     833           4 :   case PEP_REFINE_SCHEME_MBE:
     834           4 :     if (matctx->subc) {
     835           2 :       PetscCall(BVCreateVec(pep->V,&t));
     836           2 :       PetscCall(BVDuplicateResize(pep->V,PetscMax(k,nmat),&W));
     837             :     } else {
     838           2 :       W = matctx->W;
     839           2 :       PetscCall(PetscObjectReference((PetscObject)W));
     840           2 :       t = matctx->t;
     841           2 :       PetscCall(PetscObjectReference((PetscObject)t));
     842             :     }
     843           4 :     PetscCall(BVScale(matctx->W,0.0));
     844           4 :     PetscCall(BVScale(matctx->Wt,0.0));
     845           4 :     PetscCall(BVScale(matctx->M2,0.0));
     846           4 :     PetscCall(BVScale(matctx->M3,0.0));
     847           4 :     PetscCall(PetscMalloc1(nmat,&fh));
     848             :     break;
     849           6 :   case PEP_REFINE_SCHEME_SCHUR:
     850           6 :     PetscCall(KSPGetOperators(ksp,&M,NULL));
     851           6 :     PetscCall(MatShellGetContext(M,&ctx));
     852           6 :     PetscCall(BVCreateVec(pep->V,&t));
     853           6 :     PetscCall(BVDuplicateResize(pep->V,PetscMax(k,nmat),&W));
     854           6 :     fh = ctx->fih;
     855           6 :     break;
     856             :   }
     857          14 :   PetscCall(PetscArrayzero(dVS,2*k*k));
     858          14 :   PetscCall(PetscArrayzero(DfH,lda*k));
     859          14 :   PetscCall(STGetTransform(pep->st,&flg));
     860          14 :   if (flg) {
     861           3 :     PetscCall(PetscMalloc1(pep->nmat,&At));
     862          12 :     for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
     863          11 :   } else At = pep->A;
     864             : 
     865             :   /* Main loop for computing the i-th columns of dX and dS */
     866         127 :   for (i=0;i<k;i++) {
     867             :     /* Compute and update i-th column of the right hand side */
     868         113 :     PetscCall(PetscArrayzero(Rh,k));
     869         113 :     PetscCall(NRefRightSide(nmat,pep->pbc,At,k,pep->V,S,lds,i,H,ldh,fH,DfH,dH,dV,dVS,*rds,Rv,Rh,W,t));
     870             : 
     871             :     /* Update and solve system */
     872         113 :     PetscCall(BVGetColumn(dV,i,&dvi));
     873         113 :     PetscCall(NRefSysIter(i,pep,k,ksp,fH,S,lds,fh,H,ldh,Rv,Rh,pep->V,dvi,dH+i*k,matctx,W));
     874             :     /* Orthogonalize computed solution */
     875         113 :     PetscCall(BVOrthogonalizeVec(pep->V,dvi,dVS+i*2*k,&norm,&lindep));
     876         113 :     PetscCall(BVRestoreColumn(dV,i,&dvi));
     877         113 :     if (!lindep) {
     878         113 :       PetscCall(BVOrthogonalizeColumn(dV,i,dVS+k+i*2*k,&norm,&lindep));
     879         113 :       if (!lindep) {
     880         113 :         dVS[k+i+i*2*k] = norm;
     881         113 :         PetscCall(BVScaleColumn(dV,i,1.0/norm));
     882         113 :         (*rds)++;
     883             :       }
     884             :     }
     885             :   }
     886          14 :   PetscCall(BVSetActiveColumns(dV,0,*rds));
     887          14 :   PetscCall(VecDestroy(&t));
     888          14 :   PetscCall(VecDestroy(&Rv));
     889          14 :   PetscCall(BVDestroy(&W));
     890          14 :   if (flg) PetscCall(PetscFree(At));
     891          14 :   PetscCall(PetscFree2(DfH,Rh));
     892          14 :   if (pep->scheme!=PEP_REFINE_SCHEME_SCHUR) PetscCall(PetscFree(fh));
     893          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     894             : }
     895             : 
     896          28 : static PetscErrorCode NRefOrthogStep(PEP pep,PetscInt k,PetscScalar *H,PetscInt ldh,PetscScalar *fH,PetscScalar *S,PetscInt lds)
     897             : {
     898          28 :   PetscInt       j,nmat=pep->nmat,deg=nmat-1,lda=nmat*k,ldg;
     899          28 :   PetscScalar    *G,*tau,sone=1.0,zero=0.0,*work;
     900          28 :   PetscBLASInt   lds_,k_,ldh_,info,ldg_,lda_;
     901             : 
     902          28 :   PetscFunctionBegin;
     903          28 :   PetscCall(PetscMalloc3(k,&tau,k,&work,deg*k*k,&G));
     904          28 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     905          28 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     906          28 :   PetscCall(PetscBLASIntCast(k,&k_));
     907             : 
     908             :   /* Form auxiliary matrix for the orthogonalization step */
     909          28 :   ldg = deg*k;
     910          28 :   PetscCall(PEPEvaluateBasisforMatrix(pep,nmat,k,H,ldh,fH));
     911          28 :   PetscCall(PetscBLASIntCast(ldg,&ldg_));
     912          28 :   PetscCall(PetscBLASIntCast(ldh,&ldh_));
     913          84 :   for (j=0;j<deg;j++) {
     914          56 :     PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&sone,S,&lds_,fH+j*k,&lda_,&zero,G+j*k,&ldg_));
     915             :   }
     916             :   /* Orthogonalize and update S */
     917          28 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     918          28 :   PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&ldg_,&k_,G,&ldg_,tau,work,&k_,&info));
     919          28 :   PetscCall(PetscFPTrapPop());
     920          28 :   SlepcCheckLapackInfo("geqrf",info);
     921          28 :   PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&k_,&k_,&sone,G,&ldg_,S,&lds_));
     922             : 
     923             :   /* Update H */
     924          28 :   PetscCallBLAS("BLAStrmm",BLAStrmm_("L","U","N","N",&k_,&k_,&sone,G,&ldg_,H,&ldh_));
     925          28 :   PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&k_,&k_,&sone,G,&ldg_,H,&ldh_));
     926          28 :   PetscCall(PetscFree3(tau,work,G));
     927          28 :   PetscFunctionReturn(PETSC_SUCCESS);
     928             : }
     929             : 
     930          14 : static PetscErrorCode PEPNRefUpdateInvPair(PEP pep,PetscInt k,PetscScalar *H,PetscInt ldh,PetscScalar *fH,PetscScalar *dH,PetscScalar *S,PetscInt lds,BV dV,PetscScalar *dVS,PetscInt rds)
     931             : {
     932          14 :   PetscInt       i,j,nmat=pep->nmat,lda=nmat*k;
     933          14 :   PetscScalar    *tau,*array,*work;
     934          14 :   PetscBLASInt   lds_,k_,lda_,ldh_,kdrs_,info,k2_;
     935          14 :   Mat            M0;
     936             : 
     937          14 :   PetscFunctionBegin;
     938          14 :   PetscCall(PetscMalloc2(k,&tau,k,&work));
     939          14 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     940          14 :   PetscCall(PetscBLASIntCast(lda,&lda_));
     941          14 :   PetscCall(PetscBLASIntCast(ldh,&ldh_));
     942          14 :   PetscCall(PetscBLASIntCast(k,&k_));
     943          14 :   PetscCall(PetscBLASIntCast(2*k,&k2_));
     944          14 :   PetscCall(PetscBLASIntCast((k+rds),&kdrs_));
     945             :   /* Update H */
     946         127 :   for (j=0;j<k;j++) {
     947        1236 :     for (i=0;i<k;i++) H[i+j*ldh] -= dH[i+j*k];
     948             :   }
     949             :   /* Update V */
     950         127 :   for (j=0;j<k;j++) {
     951        1236 :     for (i=0;i<k;i++) dVS[i+j*2*k] = -dVS[i+j*2*k]+S[i+j*lds];
     952        1236 :     for (i=k;i<2*k;i++) dVS[i+j*2*k] = -dVS[i+j*2*k];
     953             :   }
     954          14 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     955          14 :   PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&kdrs_,&k_,dVS,&k2_,tau,work,&k_,&info));
     956          14 :   SlepcCheckLapackInfo("geqrf",info);
     957             :   /* Copy triangular matrix in S */
     958         127 :   for (j=0;j<k;j++) {
     959         731 :     for (i=0;i<=j;i++) S[i+j*lds] = dVS[i+j*2*k];
     960         618 :     for (i=j+1;i<k;i++) S[i+j*lds] = 0.0;
     961             :   }
     962          14 :   PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&k2_,&k_,&k_,dVS,&k2_,tau,work,&k_,&info));
     963          14 :   SlepcCheckLapackInfo("orgqr",info);
     964          14 :   PetscCall(PetscFPTrapPop());
     965          14 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M0));
     966          14 :   PetscCall(MatDenseGetArrayWrite(M0,&array));
     967         127 :   for (j=0;j<k;j++) PetscCall(PetscArraycpy(array+j*k,dVS+j*2*k,k));
     968          14 :   PetscCall(MatDenseRestoreArrayWrite(M0,&array));
     969          14 :   PetscCall(BVMultInPlace(pep->V,M0,0,k));
     970          14 :   if (rds) {
     971          14 :     PetscCall(MatDenseGetArrayWrite(M0,&array));
     972         127 :     for (j=0;j<k;j++) PetscCall(PetscArraycpy(array+j*k,dVS+k+j*2*k,rds));
     973          14 :     PetscCall(MatDenseRestoreArrayWrite(M0,&array));
     974          14 :     PetscCall(BVMultInPlace(dV,M0,0,k));
     975          14 :     PetscCall(BVMult(pep->V,1.0,1.0,dV,NULL));
     976             :   }
     977          14 :   PetscCall(MatDestroy(&M0));
     978          14 :   PetscCall(NRefOrthogStep(pep,k,H,ldh,fH,S,lds));
     979          14 :   PetscCall(PetscFree2(tau,work));
     980          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     981             : }
     982             : 
     983          14 : static PetscErrorCode PEPNRefSetUp(PEP pep,PetscInt k,PetscScalar *H,PetscInt ldh,PEP_REFINE_EXPLICIT *matctx,PetscBool ini)
     984             : {
     985          14 :   PEP_REFINE_MATSHELL *ctx;
     986          14 :   PetscScalar         t,*coef;
     987          14 :   const PetscScalar   *array;
     988          14 :   MatStructure        str;
     989          14 :   PetscInt            j,nmat=pep->nmat,n0,m0,n1,m1,n0_,m0_,n1_,m1_,N0,N1,p,*idx1,*idx2,count,si,i,l0;
     990          14 :   MPI_Comm            comm;
     991          14 :   PetscMPIInt         np;
     992          14 :   const PetscInt      *rgs0,*rgs1;
     993          14 :   Mat                 B,C,*E,*A,*At;
     994          14 :   IS                  is1,is2;
     995          14 :   Vec                 v;
     996          14 :   PetscBool           flg;
     997          14 :   Mat                 M,P;
     998             : 
     999          14 :   PetscFunctionBegin;
    1000          14 :   PetscCall(PetscMalloc1(nmat,&coef));
    1001          14 :   PetscCall(STGetTransform(pep->st,&flg));
    1002          14 :   if (flg) {
    1003           3 :     PetscCall(PetscMalloc1(pep->nmat,&At));
    1004          12 :     for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&At[i]));
    1005          11 :   } else At = pep->A;
    1006          14 :   switch (pep->scheme) {
    1007           4 :   case PEP_REFINE_SCHEME_EXPLICIT:
    1008           4 :     if (ini) {
    1009           4 :       if (matctx->subc) {
    1010           2 :         A = matctx->A;
    1011           2 :         PetscCall(PetscSubcommGetChild(matctx->subc,&comm));
    1012             :       } else {
    1013           2 :         A = At;
    1014           2 :         PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
    1015             :       }
    1016           4 :       E = matctx->E;
    1017           4 :       PetscCall(STGetMatStructure(pep->st,&str));
    1018           4 :       PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&E[0]));
    1019           4 :       j = (matctx->subc)?matctx->subc->color:0;
    1020           4 :       PetscCall(PEPEvaluateBasis(pep,H[j+j*ldh],0,coef,NULL));
    1021          12 :       for (j=1;j<nmat;j++) PetscCall(MatAXPY(E[0],coef[j],A[j],str));
    1022           4 :       PetscCall(MatCreateDense(comm,PETSC_DECIDE,PETSC_DECIDE,k,k,NULL,&E[1]));
    1023           4 :       PetscCall(MatGetOwnershipRange(E[0],&n0,&m0));
    1024           4 :       PetscCall(MatGetOwnershipRange(E[1],&n1,&m1));
    1025           4 :       PetscCall(MatGetOwnershipRangeColumn(E[0],&n0_,&m0_));
    1026           4 :       PetscCall(MatGetOwnershipRangeColumn(E[1],&n1_,&m1_));
    1027             :       /* T12 and T21 are computed from V and V*, so,
    1028             :          they must have the same column and row ranges */
    1029           4 :       PetscCheck(m0_-n0_==m0-n0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent dimensions");
    1030           4 :       PetscCall(MatCreateDense(comm,m0-n0,m1_-n1_,PETSC_DECIDE,PETSC_DECIDE,NULL,&B));
    1031           4 :       PetscCall(MatCreateDense(comm,m1-n1,m0_-n0_,PETSC_DECIDE,PETSC_DECIDE,NULL,&C));
    1032           4 :       PetscCall(MatCreateTile(1.0,E[0],1.0,B,1.0,C,1.0,E[1],&M));
    1033           4 :       PetscCall(MatDestroy(&B));
    1034           4 :       PetscCall(MatDestroy(&C));
    1035           4 :       matctx->compM1 = PETSC_TRUE;
    1036           4 :       PetscCall(MatGetSize(E[0],NULL,&N0));
    1037           4 :       PetscCall(MatGetSize(E[1],NULL,&N1));
    1038           4 :       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)M),&np));
    1039           4 :       PetscCall(MatGetOwnershipRanges(E[0],&rgs0));
    1040           4 :       PetscCall(MatGetOwnershipRanges(E[1],&rgs1));
    1041           4 :       PetscCall(PetscMalloc4(PetscMax(k,N1),&matctx->idxp,N0,&matctx->idxg,N0,&matctx->map0,N1,&matctx->map1));
    1042             :       /* Create column (and row) mapping */
    1043           8 :       for (p=0;p<np;p++) {
    1044         124 :         for (j=rgs0[p];j<rgs0[p+1];j++) matctx->map0[j] = j+rgs1[p];
    1045          27 :         for (j=rgs1[p];j<rgs1[p+1];j++) matctx->map1[j] = j+rgs0[p+1];
    1046             :       }
    1047           4 :       PetscCall(MatCreateVecs(M,NULL,&matctx->tN));
    1048           4 :       PetscCall(MatCreateVecs(matctx->E[1],NULL,&matctx->t1));
    1049           4 :       PetscCall(VecDuplicate(matctx->tN,&matctx->ttN));
    1050           4 :       if (matctx->subc) {
    1051           2 :         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pep),&np));
    1052           2 :         count = np*k;
    1053           2 :         PetscCall(PetscMalloc2(count,&idx1,count,&idx2));
    1054           2 :         PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pep),m1-n1,PETSC_DECIDE,&matctx->tp));
    1055           2 :         PetscCall(VecGetOwnershipRange(matctx->tp,&l0,NULL));
    1056           2 :         PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pep),k,PETSC_DECIDE,&matctx->tpg));
    1057           6 :         for (si=0;si<matctx->subc->n;si++) {
    1058           4 :           if (matctx->subc->color==si) {
    1059             :             j=0;
    1060             :             if (matctx->subc->color==si) {
    1061           6 :               for (p=0;p<np;p++) {
    1062          28 :                 for (i=n1;i<m1;i++) {
    1063          24 :                   idx1[j] = l0+i-n1;
    1064          24 :                   idx2[j++] =p*k+i;
    1065             :                 }
    1066             :               }
    1067             :             }
    1068           2 :             count = np*(m1-n1);
    1069             :           } else count =0;
    1070           4 :           PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),count,idx1,PETSC_COPY_VALUES,&is1));
    1071           4 :           PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),count,idx2,PETSC_COPY_VALUES,&is2));
    1072           4 :           PetscCall(VecScatterCreate(matctx->tp,is1,matctx->tpg,is2,&matctx->scatterp_id[si]));
    1073           4 :           PetscCall(ISDestroy(&is1));
    1074           4 :           PetscCall(ISDestroy(&is2));
    1075             :         }
    1076           2 :         PetscCall(PetscFree2(idx1,idx2));
    1077           2 :       } else PetscCall(VecScatterCreateToAll(matctx->t1,&matctx->scatterctx,&matctx->vseq));
    1078           4 :       P = M;
    1079             :     } else {
    1080           0 :       if (matctx->subc) {
    1081             :         /* Scatter vectors pep->V */
    1082           0 :         for (i=0;i<k;i++) {
    1083           0 :           PetscCall(BVGetColumn(pep->V,i,&v));
    1084           0 :           PetscCall(VecScatterBegin(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
    1085           0 :           PetscCall(VecScatterEnd(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
    1086           0 :           PetscCall(BVRestoreColumn(pep->V,i,&v));
    1087           0 :           PetscCall(VecGetArrayRead(matctx->tg,&array));
    1088           0 :           PetscCall(VecPlaceArray(matctx->t,(const PetscScalar*)array));
    1089           0 :           PetscCall(BVInsertVec(matctx->V,i,matctx->t));
    1090           0 :           PetscCall(VecResetArray(matctx->t));
    1091           0 :           PetscCall(VecRestoreArrayRead(matctx->tg,&array));
    1092             :         }
    1093             :       }
    1094             :     }
    1095             :     break;
    1096           4 :   case PEP_REFINE_SCHEME_MBE:
    1097           4 :     if (ini) {
    1098           4 :       if (matctx->subc) {
    1099           2 :         A = matctx->A;
    1100           2 :         PetscCall(PetscSubcommGetChild(matctx->subc,&comm));
    1101             :       } else {
    1102           2 :         matctx->V = pep->V;
    1103           2 :         A = At;
    1104           2 :         PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
    1105           2 :         PetscCall(MatCreateVecs(pep->A[0],&matctx->t,NULL));
    1106             :       }
    1107           4 :       PetscCall(STGetMatStructure(pep->st,&str));
    1108           4 :       PetscCall(MatDuplicate(A[0],MAT_COPY_VALUES,&matctx->M1));
    1109           4 :       j = (matctx->subc)?matctx->subc->color:0;
    1110           4 :       PetscCall(PEPEvaluateBasis(pep,H[j+j*ldh],0,coef,NULL));
    1111          12 :       for (j=1;j<nmat;j++) PetscCall(MatAXPY(matctx->M1,coef[j],A[j],str));
    1112           4 :       PetscCall(BVDuplicateResize(matctx->V,PetscMax(k,pep->nmat),&matctx->W));
    1113           4 :       PetscCall(BVDuplicateResize(matctx->V,k,&matctx->M2));
    1114           4 :       PetscCall(BVDuplicate(matctx->M2,&matctx->M3));
    1115           4 :       PetscCall(BVDuplicate(matctx->M2,&matctx->Wt));
    1116           4 :       PetscCall(PetscMalloc5(k*k,&matctx->M4,k*k,&matctx->w,k*k,&matctx->wt,k,&matctx->d,k,&matctx->dt));
    1117           4 :       matctx->compM1 = PETSC_TRUE;
    1118           4 :       M = matctx->M1;
    1119           4 :       P = M;
    1120             :     }
    1121             :     break;
    1122           6 :   case PEP_REFINE_SCHEME_SCHUR:
    1123           6 :     if (ini) {
    1124           6 :       PetscCall(PetscObjectGetComm((PetscObject)pep,&comm));
    1125           6 :       PetscCall(MatGetSize(At[0],&m0,&n0));
    1126           6 :       PetscCall(PetscMalloc1(1,&ctx));
    1127           6 :       PetscCall(STGetMatStructure(pep->st,&str));
    1128             :       /* Create a shell matrix to solve the linear system */
    1129           6 :       ctx->V = pep->V;
    1130           6 :       ctx->k = k; ctx->nmat = nmat;
    1131           6 :       PetscCall(PetscMalloc5(nmat,&ctx->A,k*k,&ctx->M4,k,&ctx->pM4,2*k*k,&ctx->work,nmat,&ctx->fih));
    1132          24 :       for (i=0;i<nmat;i++) ctx->A[i] = At[i];
    1133           6 :       PetscCall(PetscArrayzero(ctx->M4,k*k));
    1134           6 :       PetscCall(MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,m0,n0,ctx,&M));
    1135           6 :       PetscCall(MatShellSetOperation(M,MATOP_MULT,(void(*)(void))MatMult_FS));
    1136           6 :       PetscCall(BVDuplicateResize(ctx->V,PetscMax(k,pep->nmat),&ctx->W));
    1137           6 :       PetscCall(BVDuplicateResize(ctx->V,k,&ctx->M2));
    1138           6 :       PetscCall(BVDuplicate(ctx->M2,&ctx->M3));
    1139           6 :       PetscCall(BVCreateVec(pep->V,&ctx->t));
    1140           6 :       PetscCall(MatDuplicate(At[0],MAT_COPY_VALUES,&ctx->M1));
    1141           6 :       PetscCall(PEPEvaluateBasis(pep,H[0],0,coef,NULL));
    1142          18 :       for (j=1;j<nmat;j++) PetscCall(MatAXPY(ctx->M1,coef[j],At[j],str));
    1143           6 :       PetscCall(MatDuplicate(At[0],MAT_COPY_VALUES,&P));
    1144             :       /* Compute a precond matrix for the system */
    1145           6 :       t = H[0];
    1146           6 :       PetscCall(PEPEvaluateBasis(pep,t,0,coef,NULL));
    1147          18 :       for (j=1;j<nmat;j++) PetscCall(MatAXPY(P,coef[j],At[j],str));
    1148           6 :       ctx->compM1 = PETSC_TRUE;
    1149             :     }
    1150             :     break;
    1151             :   }
    1152          14 :   if (ini) {
    1153          14 :     PetscCall(PEPRefineGetKSP(pep,&pep->refineksp));
    1154          14 :     PetscCall(KSPSetErrorIfNotConverged(pep->refineksp,PETSC_TRUE));
    1155          14 :     PetscCall(PEP_KSPSetOperators(pep->refineksp,M,P));
    1156          14 :     PetscCall(KSPSetFromOptions(pep->refineksp));
    1157             :   }
    1158             : 
    1159          14 :   if (!ini && matctx && matctx->subc) {
    1160             :      /* Scatter vectors pep->V */
    1161           0 :     for (i=0;i<k;i++) {
    1162           0 :       PetscCall(BVGetColumn(pep->V,i,&v));
    1163           0 :       PetscCall(VecScatterBegin(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
    1164           0 :       PetscCall(VecScatterEnd(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
    1165           0 :       PetscCall(BVRestoreColumn(pep->V,i,&v));
    1166           0 :       PetscCall(VecGetArrayRead(matctx->tg,&array));
    1167           0 :       PetscCall(VecPlaceArray(matctx->t,(const PetscScalar*)array));
    1168           0 :       PetscCall(BVInsertVec(matctx->V,i,matctx->t));
    1169           0 :       PetscCall(VecResetArray(matctx->t));
    1170           0 :       PetscCall(VecRestoreArrayRead(matctx->tg,&array));
    1171             :     }
    1172             :    }
    1173          14 :   PetscCall(PetscFree(coef));
    1174          14 :   if (flg) PetscCall(PetscFree(At));
    1175          14 :   PetscFunctionReturn(PETSC_SUCCESS);
    1176             : }
    1177             : 
    1178           4 : static PetscErrorCode NRefSubcommSetup(PEP pep,PetscInt k,PEP_REFINE_EXPLICIT *matctx,PetscInt nsubc)
    1179             : {
    1180           4 :   PetscInt          i,si,j,m0,n0,nloc0,nloc_sub,*idx1,*idx2;
    1181           4 :   IS                is1,is2;
    1182           4 :   BVType            type;
    1183           4 :   Vec               v;
    1184           4 :   const PetscScalar *array;
    1185           4 :   Mat               *A;
    1186           4 :   PetscBool         flg;
    1187           4 :   MPI_Comm          contpar,child;
    1188             : 
    1189           4 :   PetscFunctionBegin;
    1190           4 :   PetscCall(STGetTransform(pep->st,&flg));
    1191           4 :   if (flg) {
    1192           0 :     PetscCall(PetscMalloc1(pep->nmat,&A));
    1193           0 :     for (i=0;i<pep->nmat;i++) PetscCall(STGetMatrixTransformed(pep->st,i,&A[i]));
    1194           4 :   } else A = pep->A;
    1195           4 :   PetscCall(PetscSubcommGetChild(matctx->subc,&child));
    1196           4 :   PetscCall(PetscSubcommGetContiguousParent(matctx->subc,&contpar));
    1197             : 
    1198             :   /* Duplicate pep matrices */
    1199           4 :   PetscCall(PetscMalloc3(pep->nmat,&matctx->A,nsubc,&matctx->scatter_id,nsubc,&matctx->scatterp_id));
    1200          16 :   for (i=0;i<pep->nmat;i++) PetscCall(MatCreateRedundantMatrix(A[i],0,child,MAT_INITIAL_MATRIX,&matctx->A[i]));
    1201             : 
    1202             :   /* Create Scatter */
    1203           4 :   PetscCall(MatCreateVecs(matctx->A[0],&matctx->t,NULL));
    1204           4 :   PetscCall(MatGetLocalSize(matctx->A[0],&nloc_sub,NULL));
    1205           4 :   PetscCall(VecCreateMPI(contpar,nloc_sub,PETSC_DECIDE,&matctx->tg));
    1206           4 :   PetscCall(BVGetColumn(pep->V,0,&v));
    1207           4 :   PetscCall(VecGetOwnershipRange(v,&n0,&m0));
    1208           4 :   nloc0 = m0-n0;
    1209           4 :   PetscCall(PetscMalloc2(matctx->subc->n*nloc0,&idx1,matctx->subc->n*nloc0,&idx2));
    1210             :   j = 0;
    1211          12 :   for (si=0;si<matctx->subc->n;si++) {
    1212         128 :     for (i=n0;i<m0;i++) {
    1213         120 :       idx1[j]   = i;
    1214         120 :       idx2[j++] = i+pep->n*si;
    1215             :     }
    1216             :   }
    1217           4 :   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),matctx->subc->n*nloc0,idx1,PETSC_COPY_VALUES,&is1));
    1218           4 :   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),matctx->subc->n*nloc0,idx2,PETSC_COPY_VALUES,&is2));
    1219           4 :   PetscCall(VecScatterCreate(v,is1,matctx->tg,is2,&matctx->scatter_sub));
    1220           4 :   PetscCall(ISDestroy(&is1));
    1221           4 :   PetscCall(ISDestroy(&is2));
    1222          12 :   for (si=0;si<matctx->subc->n;si++) {
    1223           8 :     j=0;
    1224         128 :     for (i=n0;i<m0;i++) {
    1225         120 :       idx1[j] = i;
    1226         120 :       idx2[j++] = i+pep->n*si;
    1227             :     }
    1228           8 :     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),nloc0,idx1,PETSC_COPY_VALUES,&is1));
    1229           8 :     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pep),nloc0,idx2,PETSC_COPY_VALUES,&is2));
    1230           8 :     PetscCall(VecScatterCreate(v,is1,matctx->tg,is2,&matctx->scatter_id[si]));
    1231           8 :     PetscCall(ISDestroy(&is1));
    1232           8 :     PetscCall(ISDestroy(&is2));
    1233             :   }
    1234           4 :   PetscCall(BVRestoreColumn(pep->V,0,&v));
    1235           4 :   PetscCall(PetscFree2(idx1,idx2));
    1236             : 
    1237             :   /* Duplicate pep->V vecs */
    1238           4 :   PetscCall(BVGetType(pep->V,&type));
    1239           4 :   PetscCall(BVCreate(child,&matctx->V));
    1240           4 :   PetscCall(BVSetType(matctx->V,type));
    1241           4 :   PetscCall(BVSetSizesFromVec(matctx->V,matctx->t,k));
    1242           4 :   if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) PetscCall(BVDuplicateResize(matctx->V,PetscMax(k,pep->nmat),&matctx->W));
    1243          28 :   for (i=0;i<k;i++) {
    1244          24 :     PetscCall(BVGetColumn(pep->V,i,&v));
    1245          24 :     PetscCall(VecScatterBegin(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
    1246          24 :     PetscCall(VecScatterEnd(matctx->scatter_sub,v,matctx->tg,INSERT_VALUES,SCATTER_FORWARD));
    1247          24 :     PetscCall(BVRestoreColumn(pep->V,i,&v));
    1248          24 :     PetscCall(VecGetArrayRead(matctx->tg,&array));
    1249          24 :     PetscCall(VecPlaceArray(matctx->t,(const PetscScalar*)array));
    1250          24 :     PetscCall(BVInsertVec(matctx->V,i,matctx->t));
    1251          24 :     PetscCall(VecResetArray(matctx->t));
    1252          24 :     PetscCall(VecRestoreArrayRead(matctx->tg,&array));
    1253             :   }
    1254             : 
    1255           4 :   PetscCall(VecDuplicate(matctx->t,&matctx->Rv));
    1256           4 :   PetscCall(VecDuplicate(matctx->t,&matctx->Vi));
    1257           4 :   if (flg) PetscCall(PetscFree(A));
    1258           4 :   PetscFunctionReturn(PETSC_SUCCESS);
    1259             : }
    1260             : 
    1261           4 : static PetscErrorCode NRefSubcommDestroy(PEP pep,PEP_REFINE_EXPLICIT *matctx)
    1262             : {
    1263           4 :   PetscInt       i;
    1264             : 
    1265           4 :   PetscFunctionBegin;
    1266           4 :   PetscCall(VecScatterDestroy(&matctx->scatter_sub));
    1267          12 :   for (i=0;i<matctx->subc->n;i++) PetscCall(VecScatterDestroy(&matctx->scatter_id[i]));
    1268          16 :   for (i=0;i<pep->nmat;i++) PetscCall(MatDestroy(&matctx->A[i]));
    1269           4 :   if (pep->scheme==PEP_REFINE_SCHEME_EXPLICIT) {
    1270           6 :     for (i=0;i<matctx->subc->n;i++) PetscCall(VecScatterDestroy(&matctx->scatterp_id[i]));
    1271           2 :     PetscCall(VecDestroy(&matctx->tp));
    1272           2 :     PetscCall(VecDestroy(&matctx->tpg));
    1273           2 :     PetscCall(BVDestroy(&matctx->W));
    1274             :   }
    1275           4 :   PetscCall(PetscFree3(matctx->A,matctx->scatter_id,matctx->scatterp_id));
    1276           4 :   PetscCall(BVDestroy(&matctx->V));
    1277           4 :   PetscCall(VecDestroy(&matctx->t));
    1278           4 :   PetscCall(VecDestroy(&matctx->tg));
    1279           4 :   PetscCall(VecDestroy(&matctx->Rv));
    1280           4 :   PetscCall(VecDestroy(&matctx->Vi));
    1281           4 :   PetscFunctionReturn(PETSC_SUCCESS);
    1282             : }
    1283             : 
    1284          14 : PetscErrorCode PEPNewtonRefinement_TOAR(PEP pep,PetscScalar sigma,PetscInt *maxits,PetscReal *tol,PetscInt k,PetscScalar *S,PetscInt lds)
    1285             : {
    1286          14 :   PetscScalar         *H,*work,*dH,*fH,*dVS;
    1287          14 :   PetscInt            ldh,i,j,its=1,nmat=pep->nmat,nsubc=pep->npart,rds;
    1288          14 :   PetscBLASInt        k_,ld_,*p,info;
    1289          14 :   BV                  dV;
    1290          14 :   PetscBool           sinvert,flg;
    1291          14 :   PEP_REFINE_EXPLICIT *matctx=NULL;
    1292          14 :   Vec                 v;
    1293          14 :   Mat                 M,P;
    1294          14 :   PEP_REFINE_MATSHELL *ctx;
    1295             : 
    1296          14 :   PetscFunctionBegin;
    1297          14 :   PetscCall(PetscLogEventBegin(PEP_Refine,pep,0,0,0));
    1298          14 :   PetscCheck(k<=pep->n,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Multiple Refinement available only for invariant pairs of dimension smaller than n=%" PetscInt_FMT,pep->n);
    1299             :   /* the input tolerance is not being taken into account (by the moment) */
    1300          14 :   its = *maxits;
    1301          14 :   PetscCall(PetscMalloc3(k*k,&dH,nmat*k*k,&fH,k,&work));
    1302          14 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldh));
    1303          14 :   PetscCall(PetscMalloc1(2*k*k,&dVS));
    1304          14 :   PetscCall(STGetTransform(pep->st,&flg));
    1305          14 :   if (!flg && pep->st && pep->ops->backtransform) { /* BackTransform */
    1306          11 :     PetscCall(PetscBLASIntCast(k,&k_));
    1307          11 :     PetscCall(PetscBLASIntCast(ldh,&ld_));
    1308          11 :     PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert));
    1309          11 :     if (sinvert) {
    1310           7 :       PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
    1311           7 :       PetscCall(PetscMalloc1(k,&p));
    1312           7 :       PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
    1313           7 :       PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&k_,&k_,H,&ld_,p,&info));
    1314           7 :       SlepcCheckLapackInfo("getrf",info);
    1315           7 :       PetscCallBLAS("LAPACKgetri",LAPACKgetri_(&k_,H,&ld_,p,work,&k_,&info));
    1316           7 :       SlepcCheckLapackInfo("getri",info);
    1317           7 :       PetscCall(PetscFPTrapPop());
    1318           7 :       PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
    1319           7 :       pep->ops->backtransform = NULL;
    1320             :     }
    1321          11 :     if (sigma!=0.0) {
    1322           7 :       PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
    1323          43 :       for (i=0;i<k;i++) H[i+ldh*i] += sigma;
    1324           7 :       PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
    1325           7 :       pep->ops->backtransform = NULL;
    1326             :     }
    1327             :   }
    1328          14 :   if ((pep->scale==PEP_SCALE_BOTH || pep->scale==PEP_SCALE_SCALAR) && pep->sfactor!=1.0) {
    1329           0 :     PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
    1330           0 :     for (j=0;j<k;j++) {
    1331           0 :       for (i=0;i<k;i++) H[i+j*ldh] *= pep->sfactor;
    1332             :     }
    1333           0 :     PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
    1334           0 :     if (!flg) {
    1335             :       /* Restore original values */
    1336           0 :       for (i=0;i<pep->nmat;i++) {
    1337           0 :         pep->pbc[pep->nmat+i] *= pep->sfactor;
    1338           0 :         pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
    1339             :       }
    1340             :     }
    1341             :   }
    1342          14 :   if ((pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) && pep->Dr) {
    1343          28 :     for (i=0;i<k;i++) {
    1344          24 :       PetscCall(BVGetColumn(pep->V,i,&v));
    1345          24 :       PetscCall(VecPointwiseMult(v,v,pep->Dr));
    1346          24 :       PetscCall(BVRestoreColumn(pep->V,i,&v));
    1347             :     }
    1348             :   }
    1349          14 :   PetscCall(DSGetArray(pep->ds,DS_MAT_A,&H));
    1350             : 
    1351          14 :   PetscCall(NRefOrthogStep(pep,k,H,ldh,fH,S,lds));
    1352             :   /* check if H is in Schur form */
    1353         113 :   for (i=0;i<k-1;i++) {
    1354             : #if !defined(PETSC_USE_COMPLEX)
    1355          99 :     PetscCheck(H[i+1+i*ldh]==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Iterative Refinement requires the complex Schur form of the projected matrix");
    1356             : #else
    1357             :     PetscCheck(H[i+1+i*ldh]==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Iterative Refinement requires an upper triangular projected matrix");
    1358             : #endif
    1359             :   }
    1360          14 :   PetscCheck(nsubc<=k,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Number of subcommunicators should not be larger than the invariant pair dimension");
    1361          14 :   PetscCall(BVSetActiveColumns(pep->V,0,k));
    1362          14 :   PetscCall(BVDuplicateResize(pep->V,k,&dV));
    1363          14 :   if (pep->scheme!=PEP_REFINE_SCHEME_SCHUR) {
    1364           8 :     PetscCall(PetscMalloc1(1,&matctx));
    1365           8 :     if (nsubc>1) { /* splitting in subcommunicators */
    1366           4 :       matctx->subc = pep->refinesubc;
    1367           4 :       PetscCall(NRefSubcommSetup(pep,k,matctx,nsubc));
    1368           4 :     } else matctx->subc=NULL;
    1369             :   }
    1370             : 
    1371             :   /* Loop performing iterative refinements */
    1372          28 :   for (i=0;i<its;i++) {
    1373             :     /* Pre-compute the polynomial basis evaluated in H */
    1374          14 :     PetscCall(PEPEvaluateBasisforMatrix(pep,nmat,k,H,ldh,fH));
    1375          14 :     PetscCall(PEPNRefSetUp(pep,k,H,ldh,matctx,PetscNot(i)));
    1376             :     /* Solve the linear system */
    1377          14 :     PetscCall(PEPNRefForwardSubstitution(pep,k,S,lds,H,ldh,fH,dV,dVS,&rds,dH,k,pep->refineksp,matctx));
    1378             :     /* Update X (=V*S) and H, and orthogonalize [X;X*fH1;...;XfH(deg-1)] */
    1379          14 :     PetscCall(PEPNRefUpdateInvPair(pep,k,H,ldh,fH,dH,S,lds,dV,dVS,rds));
    1380             :   }
    1381          14 :   PetscCall(DSRestoreArray(pep->ds,DS_MAT_A,&H));
    1382          14 :   if (!flg && sinvert) PetscCall(PetscFree(p));
    1383          14 :   PetscCall(PetscFree3(dH,fH,work));
    1384          14 :   PetscCall(PetscFree(dVS));
    1385          14 :   PetscCall(BVDestroy(&dV));
    1386          14 :   switch (pep->scheme) {
    1387             :   case PEP_REFINE_SCHEME_EXPLICIT:
    1388          12 :     for (i=0;i<2;i++) PetscCall(MatDestroy(&matctx->E[i]));
    1389           4 :     PetscCall(PetscFree4(matctx->idxp,matctx->idxg,matctx->map0,matctx->map1));
    1390           4 :     PetscCall(VecDestroy(&matctx->tN));
    1391           4 :     PetscCall(VecDestroy(&matctx->ttN));
    1392           4 :     PetscCall(VecDestroy(&matctx->t1));
    1393           4 :     if (nsubc>1) PetscCall(NRefSubcommDestroy(pep,matctx));
    1394             :     else {
    1395           2 :       PetscCall(VecDestroy(&matctx->vseq));
    1396           2 :       PetscCall(VecScatterDestroy(&matctx->scatterctx));
    1397             :     }
    1398           4 :     PetscCall(PetscFree(matctx));
    1399           4 :     PetscCall(KSPGetOperators(pep->refineksp,&M,NULL));
    1400           4 :     PetscCall(MatDestroy(&M));
    1401             :     break;
    1402           4 :   case PEP_REFINE_SCHEME_MBE:
    1403           4 :     PetscCall(BVDestroy(&matctx->W));
    1404           4 :     PetscCall(BVDestroy(&matctx->Wt));
    1405           4 :     PetscCall(BVDestroy(&matctx->M2));
    1406           4 :     PetscCall(BVDestroy(&matctx->M3));
    1407           4 :     PetscCall(MatDestroy(&matctx->M1));
    1408           4 :     PetscCall(VecDestroy(&matctx->t));
    1409           4 :     PetscCall(PetscFree5(matctx->M4,matctx->w,matctx->wt,matctx->d,matctx->dt));
    1410           4 :     if (nsubc>1) PetscCall(NRefSubcommDestroy(pep,matctx));
    1411           4 :     PetscCall(PetscFree(matctx));
    1412             :     break;
    1413           6 :   case PEP_REFINE_SCHEME_SCHUR:
    1414           6 :     PetscCall(KSPGetOperators(pep->refineksp,&M,&P));
    1415           6 :     PetscCall(MatShellGetContext(M,&ctx));
    1416           6 :     PetscCall(PetscFree5(ctx->A,ctx->M4,ctx->pM4,ctx->work,ctx->fih));
    1417           6 :     PetscCall(MatDestroy(&ctx->M1));
    1418           6 :     PetscCall(BVDestroy(&ctx->M2));
    1419           6 :     PetscCall(BVDestroy(&ctx->M3));
    1420           6 :     PetscCall(BVDestroy(&ctx->W));
    1421           6 :     PetscCall(VecDestroy(&ctx->t));
    1422           6 :     PetscCall(PetscFree(ctx));
    1423           6 :     PetscCall(MatDestroy(&M));
    1424           6 :     PetscCall(MatDestroy(&P));
    1425             :     break;
    1426             :   }
    1427          14 :   PetscCall(PetscLogEventEnd(PEP_Refine,pep,0,0,0));
    1428          14 :   PetscFunctionReturn(PETSC_SUCCESS);
    1429             : }

Generated by: LCOV version 1.14