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: "krylovschur"
      12             : 
      13             :    Method: thick-restarted Lanczos for Bethe-Salpeter pseudo-Hermitan matrices
      14             : 
      15             :    References:
      16             : 
      17             :        [1] M. Shao et al, "A structure preserving Lanczos algorithm for computing
      18             :            the optical absorption spectrum", SIAM J. Matrix Anal. App. 39(2), 2018.
      19             : 
      20             : */
      21             : #include <slepc/private/epsimpl.h>
      22             : #include "krylovschur.h"
      23             : 
      24         958 : static PetscErrorCode Orthog_Shao(Vec x,BV U,BV V,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
      25             : {
      26         958 :   PetscInt i;
      27             : 
      28         958 :   PetscFunctionBegin;
      29         958 :   PetscCall(BVSetActiveColumns(U,0,j));
      30         958 :   PetscCall(BVSetActiveColumns(V,0,j));
      31             :   /* c = real(V^* x) ; c2 = imag(U^* x)*1i */
      32             : #if defined(PETSC_USE_COMPLEX)
      33             :   PetscCall(BVDotVecBegin(V,x,c));
      34             :   PetscCall(BVDotVecBegin(U,x,c+j));
      35             :   PetscCall(BVDotVecEnd(V,x,c));
      36             :   PetscCall(BVDotVecEnd(U,x,c+j));
      37             : #else
      38         958 :   PetscCall(BVDotVec(V,x,c));
      39             : #endif
      40             : #if defined(PETSC_USE_COMPLEX)
      41             :   for (i=0; i<j; i++) {
      42             :     c[i] = PetscRealPart(c[i]);
      43             :     c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
      44             :   }
      45             : #endif
      46             :   /* x = x-U*c-V*c2 */
      47         958 :   PetscCall(BVMultVec(U,-1.0,1.0,x,c));
      48             : #if defined(PETSC_USE_COMPLEX)
      49             :   PetscCall(BVMultVec(V,-1.0,1.0,x,c+j));
      50             : #endif
      51             :   /* accumulate orthog coeffs into h */
      52       23516 :   for (i=0; i<2*j; i++) h[i] += c[i];
      53         958 :   PetscFunctionReturn(PETSC_SUCCESS);
      54             : }
      55             : 
      56             : /* Orthogonalize vector x against first j vectors in U and V
      57             : v is column j-1 of V */
      58         958 : static PetscErrorCode OrthogonalizeVector_Shao(Vec x,BV U,BV V,PetscInt j,Vec v,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
      59             : {
      60         958 :   PetscReal alpha;
      61         958 :   PetscInt  i,l;
      62             : 
      63         958 :   PetscFunctionBegin;
      64         958 :   PetscCall(PetscArrayzero(h,2*j));
      65             : 
      66             :   /* Local orthogonalization */
      67         958 :   l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
      68         958 :   PetscCall(VecDotRealPart(x,v,&alpha));
      69        3147 :   for (i=l; i<j-1; i++) h[i] = beta[i];
      70         958 :   h[j-1] = alpha;
      71             :   /* x = x-U(:,l:j-1)*h(l:j-1) */
      72         958 :   PetscCall(BVSetActiveColumns(U,l,j));
      73         958 :   PetscCall(BVMultVec(U,-1.0,1.0,x,h+l));
      74             : 
      75             :   /* Full orthogonalization */
      76         958 :   PetscCall(Orthog_Shao(x,U,V,j,h,h+2*j,breakdown));
      77         958 :   PetscFunctionReturn(PETSC_SUCCESS);
      78             : }
      79             : 
      80         173 : static PetscErrorCode EPSBSELanczos_Shao(EPS eps,BV U,BV V,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
      81             : {
      82         173 :   PetscInt       j,m = *M;
      83         173 :   Vec            v,x,y,w,f,g,vecs[2];
      84         173 :   Mat            H;
      85         173 :   IS             is[2];
      86         173 :   PetscReal      nrm;
      87         173 :   PetscScalar    *hwork,lhwork[100],gamma;
      88             : 
      89         173 :   PetscFunctionBegin;
      90         173 :   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
      91         173 :   else hwork = lhwork;
      92         173 :   PetscCall(STGetMatrix(eps->st,0,&H));
      93         173 :   PetscCall(MatNestGetISs(H,is,NULL));
      94             : 
      95             :   /* create work vectors */
      96         173 :   PetscCall(BVGetColumn(V,0,&v));
      97         173 :   PetscCall(VecDuplicate(v,&w));
      98         173 :   vecs[0] = v;
      99         173 :   vecs[1] = w;
     100         173 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
     101         173 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
     102         173 :   PetscCall(BVRestoreColumn(V,0,&v));
     103             : 
     104             :   /* Normalize initial vector */
     105         173 :   if (k==0) {
     106           7 :     if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
     107           7 :     PetscCall(BVGetColumn(U,0,&x));
     108           7 :     PetscCall(BVGetColumn(V,0,&y));
     109           7 :     PetscCall(VecCopy(x,w));
     110           7 :     PetscCall(VecConjugate(w));
     111           7 :     PetscCall(VecNestSetSubVec(f,0,x));
     112           7 :     PetscCall(VecNestSetSubVec(g,0,y));
     113           7 :     PetscCall(STApply(eps->st,f,g));
     114           7 :     PetscCall(VecDot(y,x,&gamma));
     115           7 :     nrm = PetscSqrtReal(PetscRealPart(gamma));
     116           7 :     PetscCall(VecScale(x,1.0/nrm));
     117           7 :     PetscCall(VecScale(y,1.0/nrm));
     118           7 :     PetscCall(BVRestoreColumn(U,0,&x));
     119           7 :     PetscCall(BVRestoreColumn(V,0,&y));
     120             :   }
     121             : 
     122        1131 :   for (j=k;j<m;j++) {
     123             :     /* j+1 columns (indexes 0 to j) have been computed */
     124         958 :     PetscCall(BVGetColumn(V,j,&v));
     125         958 :     PetscCall(BVGetColumn(U,j+1,&x));
     126         958 :     PetscCall(BVGetColumn(V,j+1,&y));
     127         958 :     PetscCall(VecCopy(v,w));
     128         958 :     PetscCall(VecConjugate(w));
     129         958 :     PetscCall(VecScale(w,-1.0));
     130         958 :     PetscCall(VecNestSetSubVec(f,0,v));
     131         958 :     PetscCall(VecNestSetSubVec(g,0,x));
     132         958 :     PetscCall(STApply(eps->st,f,g));
     133         958 :     PetscCall(OrthogonalizeVector_Shao(x,U,V,j+1,v,beta,k,hwork,breakdown));
     134         958 :     alpha[j] = PetscRealPart(hwork[j]);
     135         958 :     PetscCall(VecCopy(x,w));
     136         958 :     PetscCall(VecConjugate(w));
     137         958 :     PetscCall(VecNestSetSubVec(f,0,x));
     138         958 :     PetscCall(VecNestSetSubVec(g,0,y));
     139         958 :     PetscCall(STApply(eps->st,f,g));
     140         958 :     PetscCall(VecDot(x,y,&gamma));
     141         958 :     beta[j] = PetscSqrtReal(PetscRealPart(gamma));
     142         958 :     PetscCall(VecScale(x,1.0/beta[j]));
     143         958 :     PetscCall(VecScale(y,1.0/beta[j]));
     144         958 :     PetscCall(BVRestoreColumn(V,j,&v));
     145         958 :     PetscCall(BVRestoreColumn(U,j+1,&x));
     146         958 :     PetscCall(BVRestoreColumn(V,j+1,&y));
     147             :   }
     148         173 :   if (4*m > 100) PetscCall(PetscFree(hwork));
     149         173 :   PetscCall(VecDestroy(&w));
     150         173 :   PetscCall(VecDestroy(&f));
     151         173 :   PetscCall(VecDestroy(&g));
     152         173 :   PetscFunctionReturn(PETSC_SUCCESS);
     153             : }
     154             : 
     155           7 : static PetscErrorCode EPSComputeVectors_BSE_Shao(EPS eps)
     156             : {
     157           7 :   Mat         H;
     158           7 :   Vec         u1,v1;
     159           7 :   BV          U,V;
     160           7 :   IS          is[2];
     161           7 :   PetscInt    k;
     162           7 :   PetscScalar lambda;
     163             : 
     164           7 :   PetscFunctionBegin;
     165           7 :   PetscCall(STGetMatrix(eps->st,0,&H));
     166           7 :   PetscCall(MatNestGetISs(H,is,NULL));
     167           7 :   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     168          66 :   for (k=0; k<eps->nconv; k++) {
     169          59 :     PetscCall(BVGetColumn(U,k,&u1));
     170          59 :     PetscCall(BVGetColumn(V,k,&v1));
     171             :     /* approx eigenvector is [    (eigr[k]*u1+v1)]
     172             :                              [conj(eigr[k]*u1-v1)]  */
     173          59 :     lambda = eps->eigr[k];
     174          59 :     PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
     175          59 :     PetscCall(VecAYPX(u1,lambda,v1));
     176          59 :     PetscCall(VecAYPX(v1,-2.0,u1));
     177          59 :     PetscCall(VecConjugate(v1));
     178          59 :     PetscCall(BVRestoreColumn(U,k,&u1));
     179          59 :     PetscCall(BVRestoreColumn(V,k,&v1));
     180             :   }
     181           7 :   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     182             :   /* Normalize eigenvectors */
     183           7 :   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
     184           7 :   PetscCall(BVNormalize(eps->V,NULL));
     185           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     186             : }
     187             : 
     188        1858 : static PetscErrorCode Orthog_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool s,PetscBool *breakdown)
     189             : {
     190        1858 :   PetscInt i;
     191             : 
     192        1858 :   PetscFunctionBegin;
     193        1858 :   PetscCall(BVSetActiveColumns(U,0,j));
     194        1858 :   PetscCall(BVSetActiveColumns(HU,0,j));
     195        1858 :   if (s) {
     196         929 :     PetscCall(BVSetActiveColumns(V,0,j));
     197         929 :     PetscCall(BVSetActiveColumns(HV,0,j));
     198             :   } else {
     199         929 :     PetscCall(BVSetActiveColumns(V,0,j-1));
     200         929 :     PetscCall(BVSetActiveColumns(HV,0,j-1));
     201             :   }
     202             : #if defined(PETSC_USE_COMPLEX)
     203             :   PetscCall(BVDotVecBegin(HU,x,c));
     204             :   if (s || j>1) PetscCall(BVDotVecBegin(HV,x,c+j));
     205             :   PetscCall(BVDotVecEnd(HU,x,c));
     206             :   if (s || j>1) PetscCall(BVDotVecEnd(HV,x,c+j));
     207             : #else
     208        1858 :   if (s) PetscCall(BVDotVec(HU,x,c));
     209         929 :   else PetscCall(BVDotVec(HV,x,c+j));
     210             : #endif
     211       23572 :   for (i=0; i<j; i++) {
     212       21714 :     if (s) {   /* c1 = 2*real(HU^* x) ; c2 = 2*imag(HV^* x)*1i */
     213             : #if defined(PETSC_USE_COMPLEX)
     214             :       c[i] = PetscRealPart(c[i]);
     215             :       c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
     216             : #else
     217       10857 :       c[j+i] = 0.0;
     218             : #endif
     219             :     } else {   /* c1 = 2*imag(HU^* x)*1i ; c2 = 2*real(HV^* x) */
     220             : #if defined(PETSC_USE_COMPLEX)
     221             :       c[i] = PetscCMPLX(0.0,PetscImaginaryPart(c[i]));
     222             :       c[j+i] = PetscRealPart(c[j+i]);
     223             : #else
     224       10857 :       c[i] = 0.0;
     225             : #endif
     226             :     }
     227             :   }
     228             :   /* x = x-U*c1-V*c2 */
     229             : #if defined(PETSC_USE_COMPLEX)
     230             :   PetscCall(BVMultVec(U,-2.0,1.0,x,c));
     231             :   PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
     232             : #else
     233        1858 :   if (s) PetscCall(BVMultVec(U,-2.0,1.0,x,c));
     234         929 :   else PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
     235             : #endif
     236             :   /* accumulate orthog coeffs into h */
     237       45286 :   for (i=0; i<2*j; i++) h[i] += 2*c[i];
     238        1858 :   PetscFunctionReturn(PETSC_SUCCESS);
     239             : }
     240             : 
     241             : /* Orthogonalize vector x against first j vectors in U and V */
     242        1858 : static PetscErrorCode OrthogonalizeVector_Gruning(Vec x,BV U,BV V,BV HU,BV HV,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool s,PetscBool *breakdown)
     243             : {
     244        1858 :   PetscInt l,i;
     245        1858 :   Vec      u;
     246             : 
     247        1858 :   PetscFunctionBegin;
     248        1858 :   PetscCall(PetscArrayzero(h,4*j));
     249             : 
     250             :   /* Local orthogonalization */
     251        1858 :   if (s) {
     252         929 :     PetscCall(BVGetColumn(U,j-1,&u));
     253         929 :     PetscCall(VecAXPY(x,-*beta,u));
     254         929 :     PetscCall(BVRestoreColumn(U,j-1,&u));
     255         929 :     h[j-1] = *beta;
     256             :   } else {
     257         929 :     l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
     258        3002 :     for (i=l; i<j-1; i++) h[j+i] = beta[i];
     259             :     /* x = x-V(:,l:j-2)*h(l:j-2) */
     260         929 :     PetscCall(BVSetActiveColumns(V,l,j-1));
     261         929 :     PetscCall(BVMultVec(V,-1.0,1.0,x,h+j+l));
     262             :   }
     263             : 
     264             :   /* Full orthogonalization */
     265        1858 :   PetscCall(Orthog_Gruning(x,U,V,HU,HV,j,h,h+2*j,s,breakdown));
     266        1858 :   PetscFunctionReturn(PETSC_SUCCESS);
     267             : }
     268             : 
     269         165 : static PetscErrorCode EPSBSELanczos_Gruning(EPS eps,BV U,BV V,BV HU,BV HV,PetscReal *beta1,PetscReal *beta2,PetscInt k,PetscInt *M,PetscBool *breakdown)
     270             : {
     271         165 :   PetscInt       j,m = *M;
     272         165 :   Vec            v,x,y,w,f,g,vecs[2];
     273         165 :   Mat            H;
     274         165 :   IS             is[2];
     275         165 :   PetscReal      nrm;
     276         165 :   PetscScalar    *hwork,lhwork[100],dot;
     277             : 
     278         165 :   PetscFunctionBegin;
     279         165 :   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
     280         165 :   else hwork = lhwork;
     281         165 :   PetscCall(STGetMatrix(eps->st,0,&H));
     282         165 :   PetscCall(MatNestGetISs(H,is,NULL));
     283             : 
     284             :   /* create work vectors */
     285         165 :   PetscCall(BVGetColumn(V,0,&v));
     286         165 :   PetscCall(VecDuplicate(v,&w));
     287         165 :   vecs[0] = v;
     288         165 :   vecs[1] = w;
     289         165 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
     290         165 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
     291         165 :   PetscCall(BVRestoreColumn(V,0,&v));
     292             : 
     293             :   /* Normalize initial vector */
     294         165 :   if (k==0) {
     295           7 :     if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
     296             :     /* y = Hmult(v1,1) */
     297           7 :     PetscCall(BVGetColumn(U,k,&x));
     298           7 :     PetscCall(BVGetColumn(HU,k,&y));
     299           7 :     PetscCall(VecCopy(x,w));
     300           7 :     PetscCall(VecConjugate(w));
     301           7 :     PetscCall(VecNestSetSubVec(f,0,x));
     302           7 :     PetscCall(VecNestSetSubVec(g,0,y));
     303           7 :     PetscCall(STApply(eps->st,f,g));
     304             :     /* nrm = sqrt(2*real(u1'*y)); */
     305           7 :     PetscCall(VecDot(x,y,&dot));
     306           7 :     nrm = PetscSqrtReal(PetscRealPart(2*dot));
     307             :     /* U(:,j) = u1/nrm; */
     308             :     /* HU(:,j) = y/nrm; */
     309           7 :     PetscCall(VecScale(x,1.0/nrm));
     310           7 :     PetscCall(VecScale(y,1.0/nrm));
     311           7 :     PetscCall(BVRestoreColumn(U,k,&x));
     312           7 :     PetscCall(BVRestoreColumn(HU,k,&y));
     313             :   }
     314             : 
     315        1094 :   for (j=k;j<m;j++) {
     316             :     /* j+1 columns (indexes 0 to j) have been computed */
     317         929 :     PetscCall(BVGetColumn(HU,j,&x));
     318         929 :     PetscCall(BVGetColumn(V,j,&v));
     319         929 :     PetscCall(BVGetColumn(HV,j,&y));
     320         929 :     PetscCall(VecCopy(x,v));
     321         929 :     PetscCall(BVRestoreColumn(HU,j,&x));
     322             :     /* v = Orthogonalize HU(:,j) */
     323         929 :     PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,beta2,k,hwork,PETSC_FALSE,breakdown));
     324             :     /* y = Hmult(v,-1) */
     325         929 :     PetscCall(VecCopy(v,w));
     326         929 :     PetscCall(VecConjugate(w));
     327         929 :     PetscCall(VecScale(w,-1.0));
     328         929 :     PetscCall(VecNestSetSubVec(f,0,v));
     329         929 :     PetscCall(VecNestSetSubVec(g,0,y));
     330         929 :     PetscCall(STApply(eps->st,f,g));
     331             :     /* beta = sqrt(2*real(v'*y)); */
     332         929 :     PetscCall(VecDot(v,y,&dot));
     333         929 :     beta1[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
     334             :     /* V(:,j) = v/beta1; */
     335             :     /* HV(:,j) = y/beta1; */
     336         929 :     PetscCall(VecScale(v,1.0/beta1[j]));
     337         929 :     PetscCall(VecScale(y,1.0/beta1[j]));
     338         929 :     PetscCall(BVRestoreColumn(V,j,&v));
     339         929 :     PetscCall(BVRestoreColumn(HV,j,&y));
     340             : 
     341         929 :     PetscCall(BVGetColumn(HV,j,&x));
     342         929 :     PetscCall(BVGetColumn(U,j+1,&v));
     343         929 :     PetscCall(BVGetColumn(HU,j+1,&y));
     344         929 :     PetscCall(VecCopy(x,v));
     345         929 :     PetscCall(BVRestoreColumn(HV,j,&x));
     346             :     /* v = Orthogonalize HV(:,j) */
     347         929 :     PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,&beta1[j],k,hwork,PETSC_TRUE,breakdown));
     348             :     /* y = Hmult(v,1) */
     349         929 :     PetscCall(VecCopy(v,w));
     350         929 :     PetscCall(VecConjugate(w));
     351         929 :     PetscCall(VecNestSetSubVec(f,0,v));
     352         929 :     PetscCall(VecNestSetSubVec(g,0,y));
     353         929 :     PetscCall(STApply(eps->st,f,g));
     354             :     /* beta = sqrt(2*real(v'*y)); */
     355         929 :     PetscCall(VecDot(v,y,&dot));
     356         929 :     beta2[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
     357             :     /* U(:,j) = v/beta2; */
     358             :     /* HU(:,j) = y/beta2; */
     359         929 :     PetscCall(VecScale(v,1.0/beta2[j]));
     360         929 :     PetscCall(VecScale(y,1.0/beta2[j]));
     361         929 :     PetscCall(BVRestoreColumn(U,j+1,&v));
     362         929 :     PetscCall(BVRestoreColumn(HU,j+1,&y));
     363             :   }
     364         165 :   if (4*m > 100) PetscCall(PetscFree(hwork));
     365         165 :   PetscCall(VecDestroy(&w));
     366         165 :   PetscCall(VecDestroy(&f));
     367         165 :   PetscCall(VecDestroy(&g));
     368         165 :   PetscFunctionReturn(PETSC_SUCCESS);
     369             : }
     370             : 
     371           7 : static PetscErrorCode EPSComputeVectors_BSE_Gruning(EPS eps)
     372             : {
     373           7 :   Mat         H;
     374           7 :   Vec         u1,v1;
     375           7 :   BV          U,V;
     376           7 :   IS          is[2];
     377           7 :   PetscInt    k;
     378             : 
     379           7 :   PetscFunctionBegin;
     380           7 :   PetscCall(STGetMatrix(eps->st,0,&H));
     381           7 :   PetscCall(MatNestGetISs(H,is,NULL));
     382           7 :   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     383             :   /* approx eigenvector [x1] is [     u1+v1       ]
     384             :                         [x2]    [conj(u1)-conj(v1)]  */
     385          68 :   for (k=0; k<eps->nconv; k++) {
     386          61 :     PetscCall(BVGetColumn(U,k,&u1));
     387          61 :     PetscCall(BVGetColumn(V,k,&v1));
     388             :     /* x1 = u1 + v1 */
     389          61 :     PetscCall(VecAXPY(u1,1.0,v1));
     390             :     /* x2 = conj(u1) - conj(v1) = conj(u1 - v1) = conj((u1 + v1) - 2*v1) */
     391          61 :     PetscCall(VecAYPX(v1,-2.0,u1));
     392          61 :     PetscCall(VecConjugate(v1));
     393          61 :     PetscCall(BVRestoreColumn(U,k,&u1));
     394          61 :     PetscCall(BVRestoreColumn(V,k,&v1));
     395             :   }
     396           7 :   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     397             :   /* Normalize eigenvectors */
     398           7 :   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
     399           7 :   PetscCall(BVNormalize(eps->V,NULL));
     400           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     401             : }
     402             : 
     403         929 : static PetscErrorCode Orthog_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
     404             : {
     405         929 :   PetscInt          i;
     406         929 :   Mat               MX,MY,MXl,MYl;
     407         929 :   Vec               c1,c2,hxl,hyl,hz;
     408         929 :   PetscScalar       *cx1,*cx2;
     409         929 :   PetscMPIInt       len;
     410             : 
     411         929 :   PetscFunctionBegin;
     412         929 :   PetscCall(PetscMPIIntCast(j,&len));
     413         929 :   PetscCall(BVSetActiveColumns(X,0,j));
     414         929 :   PetscCall(BVSetActiveColumns(Y,0,j));
     415             :   /* BVTDotVec does not exist yet, implemented via MatMult operations */
     416         929 :   PetscCall(BVGetMat(X,&MX));
     417         929 :   PetscCall(BVGetMat(Y,&MY));
     418         929 :   PetscCall(MatDenseGetLocalMatrix(MX,&MXl));
     419         929 :   PetscCall(MatDenseGetLocalMatrix(MY,&MYl));
     420         929 :   PetscCall(MatCreateVecs(MXl,&c1,&hyl));
     421         929 :   PetscCall(MatCreateVecs(MXl,&c2,&hxl));
     422             : 
     423             :   /* c1 =  X^* hx - Y^* hy
     424             :    * c2 = -Y^T hx + X^T hy */
     425             : 
     426         929 :   PetscCall(VecGetLocalVector(hx,hxl));
     427         929 :   PetscCall(VecGetLocalVector(hy,hyl));
     428         929 :   PetscCall(VecDuplicate(hx,&hz));
     429             :   /* c1 = -(Y^* hy) */
     430         929 :   PetscCall(MatMultHermitianTranspose(MYl,hyl,c1));
     431         929 :   PetscCall(VecScale(c1,-1));
     432             :   /* c1 = c1 + X^* hx */
     433         929 :   PetscCall(MatMultHermitianTransposeAdd(MXl,hxl,c1,c1));
     434         929 :   PetscCall(VecGetArray(c1,&cx1));
     435        2787 :   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx1,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
     436         929 :   PetscCall(VecRestoreArray(c1,&cx1));
     437             :   /* c2 = -(Y^T hx) */
     438         929 :   PetscCall(MatMultTranspose(MYl,hxl,c2));
     439         929 :   PetscCall(VecScale(c2,-1));
     440             :   /* c2 = c2 + X^T hy */
     441         929 :   PetscCall(MatMultTransposeAdd(MXl,hyl,c2,c2));
     442         929 :   PetscCall(VecGetArray(c2,&cx2));
     443        2787 :   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
     444         929 :   PetscCall(VecRestoreArray(c2,&cx2));
     445         929 :   PetscCall(VecRestoreLocalVector(hx,hxl));
     446         929 :   PetscCall(VecRestoreLocalVector(hy,hyl));
     447             : 
     448             :   /* accumulate orthog coeffs into h */
     449         929 :   PetscCall(VecGetArrayRead(c1,(const PetscScalar**)&cx1));
     450         929 :   PetscCall(VecGetArrayRead(c2,(const PetscScalar**)&cx2));
     451       11916 :   for (i=0; i<j; i++) h[i] += cx1[i];
     452       11916 :   for (i=0; i<j; i++) h[i+j] += cx2[i];
     453         929 :   PetscCall(VecRestoreArrayRead(c1,(const PetscScalar**)&cx1));
     454         929 :   PetscCall(VecRestoreArrayRead(c2,(const PetscScalar**)&cx2));
     455             : 
     456             :   /* u = hx - X c1 - conj(Y) c2 */
     457             : 
     458             :   /* conj(Y) c2 */
     459         929 :   PetscCall(VecConjugate(c2));
     460         929 :   PetscCall(VecGetLocalVector(hz,hxl));
     461         929 :   PetscCall(MatMult(MYl,c2,hxl));
     462         929 :   PetscCall(VecConjugate(hxl));
     463             :   /* X c1 */
     464         929 :   PetscCall(MatMultAdd(MXl,c1,hxl,hxl));
     465         929 :   PetscCall(VecRestoreLocalVector(hz,hxl));
     466         929 :   PetscCall(VecAXPY(hx,-1,hz));
     467             : 
     468         929 :   PetscCall(BVRestoreMat(X,&MX));
     469         929 :   PetscCall(BVRestoreMat(Y,&MY));
     470         929 :   PetscCall(VecDestroy(&c1));
     471         929 :   PetscCall(VecDestroy(&c2));
     472         929 :   PetscCall(VecDestroy(&hxl));
     473         929 :   PetscCall(VecDestroy(&hyl));
     474         929 :   PetscCall(VecDestroy(&hz));
     475         929 :   PetscFunctionReturn(PETSC_SUCCESS);
     476             : }
     477             : 
     478             : /* Orthogonalize vector x against first j vectors in X and Y */
     479         929 : static PetscErrorCode OrthogonalizeVector_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
     480             : {
     481         929 :   PetscInt    l,i;
     482         929 :   PetscScalar alpha,alpha1,alpha2;
     483         929 :   Vec         x,y;
     484             : 
     485         929 :   PetscFunctionBegin;
     486         929 :   PetscCall(PetscArrayzero(h,2*j));
     487             : 
     488             :   /* Local orthogonalization */
     489         929 :   l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
     490             :   /* alpha = X(:,j-1)'*hx-Y(:,j-1)'*hy */
     491         929 :   PetscCall(BVGetColumn(X,j-1,&x));
     492         929 :   PetscCall(BVGetColumn(Y,j-1,&y));
     493         929 :   PetscCall(VecDotBegin(hx,x,&alpha1));
     494         929 :   PetscCall(VecDotBegin(hy,y,&alpha2));
     495         929 :   PetscCall(VecDotEnd(hx,x,&alpha1));
     496         929 :   PetscCall(VecDotEnd(hy,y,&alpha2));
     497         929 :   PetscCall(BVRestoreColumn(X,j-1,&x));
     498         929 :   PetscCall(BVRestoreColumn(Y,j-1,&y));
     499         929 :   alpha = alpha1-alpha2;
     500             :   /* Store coeffs into h */
     501        3088 :   for (i=l; i<j-1; i++) h[i] = h[j+i] = beta[i];
     502         929 :   h[j-1] = alpha;
     503         929 :   h[2*j-1] = alpha-1.0;
     504             :   /* Orthogonalize: hx = hx - X(:,l:j-1)*h1 - conj(Y(:,l:j-1))*h2 */
     505         929 :   PetscCall(BVSetActiveColumns(X,l,j));
     506         929 :   PetscCall(BVSetActiveColumns(Y,l,j));
     507         929 :   PetscCall(BVMultVec(X,-1.0,1.0,hx,h+l));
     508         929 :   PetscCall(VecConjugate(hx));
     509        4017 :   for (i=j+l; i<2*j; i++) h[j+i] = PetscConj(h[i]);
     510         929 :   PetscCall(BVMultVec(Y,-1.0,1.0,hx,h+2*j+l));
     511         929 :   PetscCall(VecConjugate(hx));
     512             : 
     513             :   /* Full orthogonalization */
     514         929 :   PetscCall(VecCopy(hx,hy));
     515         929 :   PetscCall(VecConjugate(hy));
     516         929 :   PetscCall(Orthog_ProjectedBSE(hx,hy,X,Y,j,h,h+2*j,breakdown));
     517         929 :   PetscFunctionReturn(PETSC_SUCCESS);
     518             : }
     519             : 
     520         171 : static PetscErrorCode EPSBSELanczos_ProjectedBSE(EPS eps,BV X,BV Y,Vec v,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
     521             : {
     522         171 :   PetscInt       j,m = *M;
     523         171 :   Vec            u,x,y,w,f,g,vecs[2];
     524         171 :   Mat            H;
     525         171 :   IS             is[2];
     526         171 :   PetscReal      nrm;
     527         171 :   PetscScalar    *hwork,lhwork[100],gamma;
     528             : 
     529         171 :   PetscFunctionBegin;
     530         171 :   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
     531         171 :   else hwork = lhwork;
     532         171 :   PetscCall(STGetMatrix(eps->st,0,&H));
     533         171 :   PetscCall(MatNestGetISs(H,is,NULL));
     534             : 
     535             :   /* create work vectors */
     536         171 :   PetscCall(BVGetColumn(Y,0,&u));
     537         171 :   PetscCall(VecDuplicate(u,&w));
     538         171 :   vecs[0] = u;
     539         171 :   vecs[1] = w;
     540         171 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
     541         171 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
     542         171 :   PetscCall(BVRestoreColumn(Y,0,&u));
     543             : 
     544             :   /* Normalize initial vector */
     545         171 :   if (k==0) {
     546           6 :     if (eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,0));
     547           6 :     PetscCall(BVGetColumn(X,0,&x));
     548             :     /* v = Hmult(u,1) */
     549           6 :     PetscCall(BVGetColumn(Y,0,&y));
     550           6 :     PetscCall(VecCopy(x,w));
     551           6 :     PetscCall(VecConjugate(w));
     552           6 :     PetscCall(VecNestSetSubVec(f,0,x));
     553           6 :     PetscCall(VecNestSetSubVec(g,0,y));
     554           6 :     PetscCall(STApply(eps->st,f,g));
     555             :     /* nrm = sqrt(real(u'v)) */
     556           6 :     PetscCall(VecDot(y,x,&gamma));
     557           6 :     nrm = PetscSqrtReal(PetscRealPart(gamma));
     558             :     /* u = u /(nrm*2) */
     559           6 :     PetscCall(VecScale(x,1.0/(2.0*nrm)));
     560             :     /* v = v /(nrm*2) */
     561           6 :     PetscCall(VecScale(y,1.0/(2.0*nrm)));
     562           6 :     PetscCall(VecCopy(y,v));
     563             :     /* X(:,1) = (u+v) */
     564           6 :     PetscCall(VecAXPY(x,1,y));
     565             :     /* Y(:,1) = conj(u-v) */
     566           6 :     PetscCall(VecAYPX(y,-2,x));
     567           6 :     PetscCall(VecConjugate(y));
     568           6 :     PetscCall(BVRestoreColumn(X,0,&x));
     569           6 :     PetscCall(BVRestoreColumn(Y,0,&y));
     570             :   }
     571             : 
     572        1100 :   for (j=k;j<m;j++) {
     573             :     /* j+1 columns (indexes 0 to j) have been computed */
     574         929 :     PetscCall(BVGetColumn(X,j+1,&x));
     575         929 :     PetscCall(BVGetColumn(Y,j+1,&y));
     576             :     /* u = Hmult(v,-1)*/
     577         929 :     PetscCall(VecCopy(v,w));
     578         929 :     PetscCall(VecConjugate(w));
     579         929 :     PetscCall(VecScale(w,-1.0));
     580         929 :     PetscCall(VecNestSetSubVec(f,0,v));
     581         929 :     PetscCall(VecNestSetSubVec(g,0,x));
     582         929 :     PetscCall(STApply(eps->st,f,g));
     583             :     /* hx = (u+v) */
     584         929 :     PetscCall(VecCopy(x,y));
     585         929 :     PetscCall(VecAXPY(x,1,v));
     586             :     /* hy = conj(u-v) */
     587         929 :     PetscCall(VecAXPY(y,-1,v));
     588         929 :     PetscCall(VecConjugate(y));
     589             :     /* [u,cd] = orthog(hx,hy,X(:,1:j),Y(:,1:j),opt)*/
     590         929 :     PetscCall(OrthogonalizeVector_ProjectedBSE(x,y,X,Y,j+1,beta,k,hwork,breakdown));
     591             :     /* alpha(j) = real(cd(j))-1/2 */
     592         929 :     alpha[j] = 2*(PetscRealPart(hwork[j]) - 0.5);
     593             :     /* v = Hmult(u,1) */
     594         929 :     PetscCall(VecCopy(x,w));
     595         929 :     PetscCall(VecConjugate(w));
     596         929 :     PetscCall(VecNestSetSubVec(f,0,x));
     597         929 :     PetscCall(VecNestSetSubVec(g,0,y));
     598         929 :     PetscCall(STApply(eps->st,f,g));
     599             :     /* nrm = sqrt(real(u'*v)) */
     600             :     /* beta(j) = nrm */
     601         929 :     PetscCall(VecDot(x,y,&gamma));
     602         929 :     beta[j] = 2.0*PetscSqrtReal(PetscRealPart(gamma));
     603             :     /* u = u/(nrm*2) */
     604         929 :     PetscCall(VecScale(x,1.0/beta[j]));
     605             :     /* v = v/(nrm*2) */
     606         929 :     PetscCall(VecScale(y,1.0/beta[j]));
     607         929 :     PetscCall(VecCopy(y,v));
     608             :     /* X(:,j+1) = (u+v) */
     609         929 :     PetscCall(VecAXPY(x,1,y));
     610             :     /* Y(:,j+1) = conj(u-v) */
     611         929 :     PetscCall(VecAYPX(y,-2,x));
     612         929 :     PetscCall(VecConjugate(y));
     613             :     /* Restore */
     614         929 :     PetscCall(BVRestoreColumn(X,j+1,&x));
     615         929 :     PetscCall(BVRestoreColumn(Y,j+1,&y));
     616             :   }
     617         171 :   if (4*m > 100) PetscCall(PetscFree(hwork));
     618         171 :   PetscCall(VecDestroy(&w));
     619         171 :   PetscCall(VecDestroy(&f));
     620         171 :   PetscCall(VecDestroy(&g));
     621         171 :   PetscFunctionReturn(PETSC_SUCCESS);
     622             : }
     623             : 
     624           6 : static PetscErrorCode EPSComputeVectors_BSE_ProjectedBSE(EPS eps)
     625             : {
     626           6 :   Mat         H;
     627           6 :   Vec         x1,y1,cx1,cy1;
     628           6 :   BV          X,Y;
     629           6 :   IS          is[2];
     630           6 :   PetscInt    k;
     631           6 :   PetscScalar delta1,delta2,lambda;
     632             : 
     633           6 :   PetscFunctionBegin;
     634           6 :   PetscCall(STGetMatrix(eps->st,0,&H));
     635           6 :   PetscCall(MatNestGetISs(H,is,NULL));
     636           6 :   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&X,&Y));
     637           6 :   PetscCall(BVCreateVec(X,&cx1));
     638           6 :   PetscCall(BVCreateVec(Y,&cy1));
     639          58 :   for (k=0; k<eps->nconv; k++) {
     640             :     /* approx eigenvector is [ delta1*x1 + delta2*conj(y1) ]
     641             :                              [ delta1*y1 + delta2*conj(x1) ] */
     642          52 :     lambda = eps->eigr[k];
     643          52 :     PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
     644          52 :     delta1 = lambda+1.0;
     645          52 :     delta2 = lambda-1.0;
     646          52 :     PetscCall(BVGetColumn(X,k,&x1));
     647          52 :     PetscCall(BVGetColumn(Y,k,&y1));
     648          52 :     PetscCall(VecCopy(x1,cx1));
     649          52 :     PetscCall(VecCopy(y1,cy1));
     650          52 :     PetscCall(VecConjugate(cx1));
     651          52 :     PetscCall(VecConjugate(cy1));
     652          52 :     PetscCall(VecScale(x1,delta1));
     653          52 :     PetscCall(VecScale(y1,delta1));
     654          52 :     PetscCall(VecAXPY(x1,delta2,cy1));
     655          52 :     PetscCall(VecAXPY(y1,delta2,cx1));
     656          52 :     PetscCall(BVRestoreColumn(X,k,&x1));
     657          52 :     PetscCall(BVRestoreColumn(Y,k,&y1));
     658             :   }
     659           6 :   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&X,&Y));
     660           6 :   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
     661             :   /* Normalize eigenvector */
     662           6 :   PetscCall(BVNormalize(eps->V,NULL));
     663           6 :   PetscCall(VecDestroy(&cx1));
     664           6 :   PetscCall(VecDestroy(&cy1));
     665           6 :   PetscFunctionReturn(PETSC_SUCCESS);
     666             : }
     667             : 
     668          20 : PetscErrorCode EPSSetUp_KrylovSchur_BSE(EPS eps)
     669             : {
     670          20 :   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
     671          20 :   PetscBool       flg,sinvert;
     672          20 :   PetscInt        nev;
     673             : 
     674          20 :   PetscFunctionBegin;
     675          20 :   PetscCheck((eps->problem_type==EPS_BSE),PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Problem type should be BSE");
     676          20 :   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_BALANCE,PETSC_TRUE," with BSE structure");
     677          20 :   if (eps->nev==0 && eps->stop!=EPS_STOP_THRESHOLD) eps->nev = 1;
     678          20 :   nev = (eps->nev+1)/2;
     679          20 :   PetscCall(EPSSetDimensions_Default(eps,&nev,&eps->ncv,&eps->mpd));
     680          20 :   PetscCheck(eps->ncv<=nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
     681          34 :   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
     682             : 
     683          20 :   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STSHIFT,""));
     684          20 :   PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Krylov-Schur BSE only supports shift and shift-and-invert ST");
     685          20 :   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
     686          20 :   if (!eps->which) {
     687          18 :     if (sinvert) eps->which = EPS_TARGET_MAGNITUDE;
     688          15 :     else eps->which = EPS_SMALLEST_MAGNITUDE;
     689             :   }
     690             : 
     691          20 :   if (!ctx->keep) ctx->keep = 0.5;
     692          20 :   PetscCall(STSetStructured(eps->st,PETSC_TRUE));
     693             : 
     694          20 :   PetscCall(EPSAllocateSolution(eps,1));
     695          20 :   switch (ctx->bse) {
     696           7 :     case EPS_KRYLOVSCHUR_BSE_SHAO:
     697           7 :       eps->ops->solve = EPSSolve_KrylovSchur_BSE_Shao;
     698           7 :       eps->ops->computevectors = EPSComputeVectors_BSE_Shao;
     699           7 :       PetscCall(DSSetType(eps->ds,DSHEP));
     700           7 :       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
     701           7 :       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
     702           7 :       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
     703             :       break;
     704           7 :     case EPS_KRYLOVSCHUR_BSE_GRUNING:
     705           7 :       eps->ops->solve = EPSSolve_KrylovSchur_BSE_Gruning;
     706           7 :       eps->ops->computevectors = EPSComputeVectors_BSE_Gruning;
     707           7 :       PetscCall(DSSetType(eps->ds,DSSVD));
     708           7 :       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
     709           7 :       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
     710           7 :       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
     711             :       break;
     712           6 :     case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
     713           6 :       eps->ops->solve = EPSSolve_KrylovSchur_BSE_ProjectedBSE;
     714           6 :       eps->ops->computevectors = EPSComputeVectors_BSE_ProjectedBSE;
     715           6 :       PetscCall(DSSetType(eps->ds,DSHEP));
     716           6 :       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
     717           6 :       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
     718           6 :       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
     719             :       break;
     720           0 :     default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
     721             :   }
     722          20 :   PetscFunctionReturn(PETSC_SUCCESS);
     723             : }
     724             : 
     725           7 : PetscErrorCode EPSSolve_KrylovSchur_BSE_Shao(EPS eps)
     726             : {
     727           7 :   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
     728           7 :   PetscInt        i,k,l,ld,nv,nconv=0,nevsave;
     729           7 :   Mat             H,Q;
     730           7 :   BV              U,V;
     731           7 :   IS              is[2];
     732           7 :   PetscReal       *a,*b,beta;
     733           7 :   PetscBool       breakdown=PETSC_FALSE;
     734             : 
     735           7 :   PetscFunctionBegin;
     736           7 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     737             : 
     738             :   /* Extract matrix blocks */
     739           7 :   PetscCall(STGetMatrix(eps->st,0,&H));
     740           7 :   PetscCall(MatNestGetISs(H,is,NULL));
     741             : 
     742             :   /* Get the split bases */
     743           7 :   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     744             : 
     745           7 :   nevsave  = eps->nev;
     746           7 :   eps->nev = (eps->nev+1)/2;
     747           7 :   l = 0;
     748             : 
     749             :   /* Restart loop */
     750           7 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     751         173 :     eps->its++;
     752             : 
     753             :     /* Compute an nv-step Lanczos factorization */
     754         173 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     755         173 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     756         173 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
     757         173 :     b = a + ld;
     758         173 :     PetscCall(EPSBSELanczos_Shao(eps,U,V,a,b,eps->nconv+l,&nv,&breakdown));
     759         173 :     beta = b[nv-1];
     760         173 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
     761         173 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     762         173 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
     763         173 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     764             : 
     765             :     /* Solve projected problem */
     766         173 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     767         173 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     768         173 :     PetscCall(DSUpdateExtraRow(eps->ds));
     769         173 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     770             : 
     771             :     /* Check convergence */
     772        2535 :     for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
     773         173 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
     774         173 :     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
     775         173 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     776         173 :     nconv = k;
     777             : 
     778             :     /* Update l */
     779         173 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
     780         166 :     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
     781         173 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     782         173 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     783             : 
     784         173 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
     785         166 :       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
     786             :       /* Prepare the Rayleigh quotient for restart */
     787         166 :       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
     788             :     }
     789             :     /* Update the corresponding vectors
     790             :        U(:,idx) = U*Q(:,idx),  V(:,idx) = V*Q(:,idx) */
     791         173 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
     792         173 :     PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
     793         173 :     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
     794         173 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
     795             : 
     796         173 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
     797         166 :       PetscCall(BVCopyColumn(eps->V,nv,k+l));
     798         166 :       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
     799           1 :         eps->ncv = eps->mpd+k;
     800           1 :         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     801           1 :         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
     802           1 :         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     803           7 :         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
     804           1 :         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
     805           1 :         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     806             :       }
     807             :     }
     808         173 :     eps->nconv = k;
     809         180 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
     810             :   }
     811             : 
     812           7 :   eps->nev = nevsave;
     813             : 
     814           7 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     815           7 :   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     816           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     817             : }
     818             : 
     819             : /*
     820             :    EPSConvergence_Gruning - convergence check based on SVDKrylovConvergence().
     821             : */
     822         165 : static PetscErrorCode EPSConvergence_Gruning(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
     823             : {
     824         165 :   PetscInt       k,marker,ld;
     825         165 :   PetscReal      *alpha,*beta,resnorm;
     826         165 :   PetscBool      extra;
     827             : 
     828         165 :   PetscFunctionBegin;
     829         165 :   *kout = 0;
     830         165 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     831         165 :   PetscCall(DSGetExtraRow(eps->ds,&extra));
     832         165 :   PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
     833         165 :   marker = -1;
     834         165 :   if (eps->trackall) getall = PETSC_TRUE;
     835         165 :   PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     836         165 :   beta = alpha + ld;
     837         226 :   for (k=kini;k<kini+nits;k++) {
     838         226 :     resnorm = PetscAbsReal(beta[k]);
     839         226 :     PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
     840         226 :     if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
     841         226 :     if (marker!=-1 && !getall) break;
     842             :   }
     843         165 :   PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     844         165 :   if (marker!=-1) k = marker;
     845         165 :   *kout = k;
     846         165 :   PetscFunctionReturn(PETSC_SUCCESS);
     847             : }
     848             : 
     849           7 : PetscErrorCode EPSSolve_KrylovSchur_BSE_Gruning(EPS eps)
     850             : {
     851           7 :   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
     852           7 :   PetscInt        i,k,l,ld,nv,nconv=0,nevsave;
     853           7 :   Mat             H,Q,Z;
     854           7 :   BV              U,V,HU,HV;
     855           7 :   IS              is[2];
     856           7 :   PetscReal       *d1,*d2,beta;
     857           7 :   PetscBool       breakdown=PETSC_FALSE;
     858             : 
     859           7 :   PetscFunctionBegin;
     860           7 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     861             : 
     862             :   /* Extract matrix blocks */
     863           7 :   PetscCall(STGetMatrix(eps->st,0,&H));
     864           7 :   PetscCall(MatNestGetISs(H,is,NULL));
     865             : 
     866             :   /* Get the split bases */
     867           7 :   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     868             : 
     869             :   /* Create HU and HV */
     870           7 :   PetscCall(BVDuplicate(U,&HU));
     871           7 :   PetscCall(BVDuplicate(V,&HV));
     872             : 
     873           7 :   nevsave  = eps->nev;
     874           7 :   eps->nev = (eps->nev+1)/2;
     875           7 :   l = 0;
     876             : 
     877             :   /* Restart loop */
     878           7 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     879         165 :     eps->its++;
     880             : 
     881             :     /* Compute an nv-step Lanczos factorization */
     882         165 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     883         165 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     884         165 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&d1));
     885         165 :     d2 = d1 + ld;
     886         165 :     PetscCall(EPSBSELanczos_Gruning(eps,U,V,HU,HV,d1,d2,eps->nconv+l,&nv,&breakdown));
     887         165 :     beta = d1[nv-1];
     888         165 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&d1));
     889             : 
     890             :     /* Compute SVD */
     891         165 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     892         165 :     PetscCall(DSSVDSetDimensions(eps->ds,nv));
     893         165 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
     894         165 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     895             : 
     896         165 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     897         165 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     898         165 :     PetscCall(DSUpdateExtraRow(eps->ds));
     899         165 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     900             : 
     901             :     /* Check convergence */
     902         165 :     PetscCall(EPSConvergence_Gruning(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
     903         165 :     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
     904         165 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     905         165 :     nconv = k;
     906             : 
     907             :     /* Update l */
     908         165 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
     909         158 :     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
     910         165 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     911         165 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     912             : 
     913         165 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
     914         158 :       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
     915             :       /* Prepare the Rayleigh quotient for restart */
     916         158 :       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
     917             :     }
     918             :     /* Update the corresponding vectors
     919             :        U(:,idx) = U*Q(:,idx),  V(:,idx) = V*Z(:,idx) */
     920         165 :     PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Q));
     921         165 :     PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Z));
     922         165 :     PetscCall(BVMultInPlace(U,Z,eps->nconv,k+l));
     923         165 :     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
     924         165 :     PetscCall(BVMultInPlace(HU,Z,eps->nconv,k+l));
     925         165 :     PetscCall(BVMultInPlace(HV,Q,eps->nconv,k+l));
     926         165 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Q));
     927         165 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Z));
     928             : 
     929         165 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
     930         158 :       PetscCall(BVCopyColumn(U,nv,k+l));
     931         158 :       PetscCall(BVCopyColumn(HU,nv,k+l));
     932         158 :       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
     933           1 :         eps->ncv = eps->mpd+k;
     934           1 :         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     935           1 :         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
     936           1 :         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     937           1 :         PetscCall(BVResize(HU,eps->ncv+1,PETSC_TRUE));
     938           1 :         PetscCall(BVResize(HV,eps->ncv+1,PETSC_TRUE));
     939           7 :         for (i=nv;i<eps->ncv;i++) { eps->perm[i] = i; eps->eigi[i] = 0.0; }
     940           1 :         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
     941           1 :         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     942             :       }
     943             :     }
     944         165 :     eps->nconv = k;
     945         172 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
     946             :   }
     947             : 
     948           7 :   eps->nev = nevsave;
     949             : 
     950           7 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     951           7 :   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     952             : 
     953           7 :   PetscCall(BVDestroy(&HU));
     954           7 :   PetscCall(BVDestroy(&HV));
     955           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     956             : }
     957             : 
     958           6 : PetscErrorCode EPSSolve_KrylovSchur_BSE_ProjectedBSE(EPS eps)
     959             : {
     960           6 :   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
     961           6 :   PetscInt        i,k,l,ld,nv,nconv=0,nevsave;
     962           6 :   Mat             H,Q;
     963           6 :   Vec             v;
     964           6 :   BV              U,V;
     965           6 :   IS              is[2];
     966           6 :   PetscReal       *a,*b,beta;
     967           6 :   PetscBool       breakdown=PETSC_FALSE;
     968             : 
     969           6 :   PetscFunctionBegin;
     970           6 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     971             : 
     972             :   /* Extract matrix blocks */
     973           6 :   PetscCall(STGetMatrix(eps->st,0,&H));
     974           6 :   PetscCall(MatNestGetISs(H,is,NULL));
     975             : 
     976             :   /* Get the split bases */
     977           6 :   PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     978             : 
     979             :   /* Vector v */
     980           6 :   PetscCall(BVCreateVec(V,&v));
     981             : 
     982           6 :   nevsave  = eps->nev;
     983           6 :   eps->nev = (eps->nev+1)/2;
     984           6 :   l = 0;
     985             : 
     986             :   /* Restart loop */
     987           6 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     988         171 :     eps->its++;
     989             : 
     990             :     /* Compute an nv-step Lanczos factorization */
     991         171 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     992         171 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     993         171 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
     994         171 :     b = a + ld;
     995         171 :     PetscCall(EPSBSELanczos_ProjectedBSE(eps,U,V,v,a,b,eps->nconv+l,&nv,&breakdown));
     996         171 :     beta = b[nv-1];
     997         171 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
     998         171 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     999         171 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
    1000         171 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
    1001             : 
    1002             :     /* Solve projected problem */
    1003         171 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
    1004         171 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
    1005         171 :     PetscCall(DSUpdateExtraRow(eps->ds));
    1006         171 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
    1007             : 
    1008             :     /* Check convergence */
    1009        2501 :     for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
    1010         171 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
    1011         171 :     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
    1012         171 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
    1013         171 :     nconv = k;
    1014             : 
    1015             :     /* Update l */
    1016         171 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
    1017         165 :     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
    1018         171 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
    1019         171 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
    1020             : 
    1021         171 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
    1022         165 :       PetscCheck(!breakdown,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Breakdown in BSE Krylov-Schur (beta=%g)",(double)beta);
    1023             :       /* Prepare the Rayleigh quotient for restart */
    1024         165 :       PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
    1025             :     }
    1026             :     /* Update the corresponding vectors
    1027             :        U(:,idx) = U*Q(:,idx),  V(:,idx) = V*Q(:,idx) */
    1028         171 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
    1029         171 :     PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
    1030         171 :     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
    1031         171 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
    1032             : 
    1033         171 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
    1034         165 :       PetscCall(BVCopyColumn(eps->V,nv,k+l));
    1035         165 :       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
    1036           1 :         eps->ncv = eps->mpd+k;
    1037           1 :         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
    1038           1 :         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
    1039           1 :         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
    1040           7 :         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
    1041           1 :         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
    1042           1 :         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
    1043             :       }
    1044             :     }
    1045         171 :     eps->nconv = k;
    1046         177 :     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
    1047             :   }
    1048             : 
    1049           6 :   eps->nev = nevsave;
    1050             : 
    1051           6 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
    1052           6 :   PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
    1053             : 
    1054           6 :   PetscCall(VecDestroy(&v));
    1055           6 :   PetscFunctionReturn(PETSC_SUCCESS);
    1056             : }