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        1434 : static PetscErrorCode Orthog_Shao(Vec x,BV U,BV V,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
      25             : {
      26        1434 :   PetscInt i;
      27             : 
      28        1434 :   PetscFunctionBegin;
      29        1434 :   PetscCall(BVSetActiveColumns(U,0,j));
      30        1434 :   PetscCall(BVSetActiveColumns(V,0,j));
      31             :   /* c = real(V^* x) ; c2 = imag(U^* x)*1i */
      32             : #if defined(PETSC_USE_COMPLEX)
      33        1434 :   PetscCall(BVDotVecBegin(V,x,c));
      34        1434 :   PetscCall(BVDotVecBegin(U,x,c+j));
      35        1434 :   PetscCall(BVDotVecEnd(V,x,c));
      36        1434 :   PetscCall(BVDotVecEnd(U,x,c+j));
      37             : #else
      38             :   PetscCall(BVDotVec(V,x,c));
      39             : #endif
      40             : #if defined(PETSC_USE_COMPLEX)
      41       19376 :   for (i=0; i<j; i++) {
      42       17942 :     c[i] = PetscRealPart(c[i]);
      43       17942 :     c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
      44             :   }
      45             : #endif
      46             :   /* x = x-U*c-V*c2 */
      47        1434 :   PetscCall(BVMultVec(U,-1.0,1.0,x,c));
      48             : #if defined(PETSC_USE_COMPLEX)
      49        1434 :   PetscCall(BVMultVec(V,-1.0,1.0,x,c+j));
      50             : #endif
      51             :   /* accumulate orthog coeffs into h */
      52       37318 :   for (i=0; i<2*j; i++) h[i] += c[i];
      53        1434 :   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        1434 : static PetscErrorCode OrthogonalizeVector_Shao(Vec x,BV U,BV V,PetscInt j,Vec v,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
      59             : {
      60        1434 :   PetscReal alpha;
      61        1434 :   PetscInt  i,l;
      62             : 
      63        1434 :   PetscFunctionBegin;
      64        1434 :   PetscCall(PetscArrayzero(h,2*j));
      65             : 
      66             :   /* Local orthogonalization */
      67        1434 :   l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
      68        1434 :   PetscCall(VecDotRealPart(x,v,&alpha));
      69        4920 :   for (i=l; i<j-1; i++) h[i] = beta[i];
      70        1434 :   h[j-1] = alpha;
      71             :   /* x = x-U(:,l:j-1)*h(l:j-1) */
      72        1434 :   PetscCall(BVSetActiveColumns(U,l,j));
      73        1434 :   PetscCall(BVMultVec(U,-1.0,1.0,x,h+l));
      74             : 
      75             :   /* Full orthogonalization */
      76        1434 :   PetscCall(Orthog_Shao(x,U,V,j,h,h+2*j,breakdown));
      77        1434 :   PetscFunctionReturn(PETSC_SUCCESS);
      78             : }
      79             : 
      80         260 : static PetscErrorCode EPSBSELanczos_Shao(EPS eps,BV U,BV V,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
      81             : {
      82         260 :   PetscInt       j,m = *M;
      83         260 :   Vec            v,x,y,w,f,g,vecs[2];
      84         260 :   Mat            H;
      85         260 :   IS             is[2];
      86         260 :   PetscReal      nrm;
      87         260 :   PetscScalar    *hwork,lhwork[100],gamma;
      88             : 
      89         260 :   PetscFunctionBegin;
      90         260 :   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
      91         260 :   else hwork = lhwork;
      92         260 :   PetscCall(STGetMatrix(eps->st,0,&H));
      93         260 :   PetscCall(MatNestGetISs(H,is,NULL));
      94             : 
      95             :   /* create work vectors */
      96         260 :   PetscCall(BVGetColumn(V,0,&v));
      97         260 :   PetscCall(VecDuplicate(v,&w));
      98         260 :   vecs[0] = v;
      99         260 :   vecs[1] = w;
     100         260 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
     101         260 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
     102         260 :   PetscCall(BVRestoreColumn(V,0,&v));
     103             : 
     104             :   /* Normalize initial vector */
     105         260 :   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        1694 :   for (j=k;j<m;j++) {
     123             :     /* j+1 columns (indexes 0 to j) have been computed */
     124        1434 :     PetscCall(BVGetColumn(V,j,&v));
     125        1434 :     PetscCall(BVGetColumn(U,j+1,&x));
     126        1434 :     PetscCall(BVGetColumn(V,j+1,&y));
     127        1434 :     PetscCall(VecCopy(v,w));
     128        1434 :     PetscCall(VecConjugate(w));
     129        1434 :     PetscCall(VecScale(w,-1.0));
     130        1434 :     PetscCall(VecNestSetSubVec(f,0,v));
     131        1434 :     PetscCall(VecNestSetSubVec(g,0,x));
     132        1434 :     PetscCall(STApply(eps->st,f,g));
     133        1434 :     PetscCall(OrthogonalizeVector_Shao(x,U,V,j+1,v,beta,k,hwork,breakdown));
     134        1434 :     alpha[j] = PetscRealPart(hwork[j]);
     135        1434 :     PetscCall(VecCopy(x,w));
     136        1434 :     PetscCall(VecConjugate(w));
     137        1434 :     PetscCall(VecNestSetSubVec(f,0,x));
     138        1434 :     PetscCall(VecNestSetSubVec(g,0,y));
     139        1434 :     PetscCall(STApply(eps->st,f,g));
     140        1434 :     PetscCall(VecDot(x,y,&gamma));
     141        1434 :     beta[j] = PetscSqrtReal(PetscRealPart(gamma));
     142        1434 :     PetscCall(VecScale(x,1.0/beta[j]));
     143        1434 :     PetscCall(VecScale(y,1.0/beta[j]));
     144        1434 :     PetscCall(BVRestoreColumn(V,j,&v));
     145        1434 :     PetscCall(BVRestoreColumn(U,j+1,&x));
     146        1434 :     PetscCall(BVRestoreColumn(V,j+1,&y));
     147             :   }
     148         260 :   if (4*m > 100) PetscCall(PetscFree(hwork));
     149         260 :   PetscCall(VecDestroy(&w));
     150         260 :   PetscCall(VecDestroy(&f));
     151         260 :   PetscCall(VecDestroy(&g));
     152         260 :   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          57 :   for (k=0; k<eps->nconv; k++) {
     169          50 :     PetscCall(BVGetColumn(U,k,&u1));
     170          50 :     PetscCall(BVGetColumn(V,k,&v1));
     171             :     /* approx eigenvector is [    (eigr[k]*u1+v1)]
     172             :                              [conj(eigr[k]*u1-v1)]  */
     173          50 :     lambda = eps->eigr[k];
     174          50 :     PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
     175          50 :     PetscCall(VecAYPX(u1,lambda,v1));
     176          50 :     PetscCall(VecAYPX(v1,-2.0,u1));
     177          50 :     PetscCall(VecConjugate(v1));
     178          50 :     PetscCall(BVRestoreColumn(U,k,&u1));
     179          50 :     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        2814 : 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        2814 :   PetscInt i;
     191             : 
     192        2814 :   PetscFunctionBegin;
     193        2814 :   PetscCall(BVSetActiveColumns(U,0,j));
     194        2814 :   PetscCall(BVSetActiveColumns(HU,0,j));
     195        2814 :   if (s) {
     196        1407 :     PetscCall(BVSetActiveColumns(V,0,j));
     197        1407 :     PetscCall(BVSetActiveColumns(HV,0,j));
     198             :   } else {
     199        1407 :     PetscCall(BVSetActiveColumns(V,0,j-1));
     200        1407 :     PetscCall(BVSetActiveColumns(HV,0,j-1));
     201             :   }
     202             : #if defined(PETSC_USE_COMPLEX)
     203        2814 :   PetscCall(BVDotVecBegin(HU,x,c));
     204        2814 :   if (s || j>1) PetscCall(BVDotVecBegin(HV,x,c+j));
     205        2814 :   PetscCall(BVDotVecEnd(HU,x,c));
     206        2814 :   if (s || j>1) PetscCall(BVDotVecEnd(HV,x,c+j));
     207             : #else
     208             :   if (s) PetscCall(BVDotVec(HU,x,c));
     209             :   else PetscCall(BVDotVec(HV,x,c+j));
     210             : #endif
     211       38144 :   for (i=0; i<j; i++) {
     212       35330 :     if (s) {   /* c1 = 2*real(HU^* x) ; c2 = 2*imag(HV^* x)*1i */
     213             : #if defined(PETSC_USE_COMPLEX)
     214       17665 :       c[i] = PetscRealPart(c[i]);
     215       17665 :       c[j+i] = PetscCMPLX(0.0,PetscImaginaryPart(c[j+i]));
     216             : #else
     217             :       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       17665 :       c[i] = PetscCMPLX(0.0,PetscImaginaryPart(c[i]));
     222       17665 :       c[j+i] = PetscRealPart(c[j+i]);
     223             : #else
     224             :       c[i] = 0.0;
     225             : #endif
     226             :     }
     227             :   }
     228             :   /* x = x-U*c1-V*c2 */
     229             : #if defined(PETSC_USE_COMPLEX)
     230        2814 :   PetscCall(BVMultVec(U,-2.0,1.0,x,c));
     231        2814 :   PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
     232             : #else
     233             :   if (s) PetscCall(BVMultVec(U,-2.0,1.0,x,c));
     234             :   else PetscCall(BVMultVec(V,-2.0,1.0,x,c+j));
     235             : #endif
     236             :   /* accumulate orthog coeffs into h */
     237       73474 :   for (i=0; i<2*j; i++) h[i] += 2*c[i];
     238        2814 :   PetscFunctionReturn(PETSC_SUCCESS);
     239             : }
     240             : 
     241             : /* Orthogonalize vector x against first j vectors in U and V */
     242        2814 : 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        2814 :   PetscInt l,i;
     245        2814 :   Vec      u;
     246             : 
     247        2814 :   PetscFunctionBegin;
     248        2814 :   PetscCall(PetscArrayzero(h,4*j));
     249             : 
     250             :   /* Local orthogonalization */
     251        2814 :   if (s) {
     252        1407 :     PetscCall(BVGetColumn(U,j-1,&u));
     253        1407 :     PetscCall(VecAXPY(x,-*beta,u));
     254        1407 :     PetscCall(BVRestoreColumn(U,j-1,&u));
     255        1407 :     h[j-1] = *beta;
     256             :   } else {
     257        1407 :     l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
     258        4842 :     for (i=l; i<j-1; i++) h[j+i] = beta[i];
     259             :     /* x = x-V(:,l:j-2)*h(l:j-2) */
     260        1407 :     PetscCall(BVSetActiveColumns(V,l,j-1));
     261        1407 :     PetscCall(BVMultVec(V,-1.0,1.0,x,h+j+l));
     262             :   }
     263             : 
     264             :   /* Full orthogonalization */
     265        2814 :   PetscCall(Orthog_Gruning(x,U,V,HU,HV,j,h,h+2*j,s,breakdown));
     266        2814 :   PetscFunctionReturn(PETSC_SUCCESS);
     267             : }
     268             : 
     269         255 : 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         255 :   PetscInt       j,m = *M;
     272         255 :   Vec            v,x,y,w,f,g,vecs[2];
     273         255 :   Mat            H;
     274         255 :   IS             is[2];
     275         255 :   PetscReal      nrm;
     276         255 :   PetscScalar    *hwork,lhwork[100],dot;
     277             : 
     278         255 :   PetscFunctionBegin;
     279         255 :   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
     280         255 :   else hwork = lhwork;
     281         255 :   PetscCall(STGetMatrix(eps->st,0,&H));
     282         255 :   PetscCall(MatNestGetISs(H,is,NULL));
     283             : 
     284             :   /* create work vectors */
     285         255 :   PetscCall(BVGetColumn(V,0,&v));
     286         255 :   PetscCall(VecDuplicate(v,&w));
     287         255 :   vecs[0] = v;
     288         255 :   vecs[1] = w;
     289         255 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
     290         255 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
     291         255 :   PetscCall(BVRestoreColumn(V,0,&v));
     292             : 
     293             :   /* Normalize initial vector */
     294         255 :   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        1662 :   for (j=k;j<m;j++) {
     316             :     /* j+1 columns (indexes 0 to j) have been computed */
     317        1407 :     PetscCall(BVGetColumn(HU,j,&x));
     318        1407 :     PetscCall(BVGetColumn(V,j,&v));
     319        1407 :     PetscCall(BVGetColumn(HV,j,&y));
     320        1407 :     PetscCall(VecCopy(x,v));
     321        1407 :     PetscCall(BVRestoreColumn(HU,j,&x));
     322             :     /* v = Orthogonalize HU(:,j) */
     323        1407 :     PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,beta2,k,hwork,PETSC_FALSE,breakdown));
     324             :     /* y = Hmult(v,-1) */
     325        1407 :     PetscCall(VecCopy(v,w));
     326        1407 :     PetscCall(VecConjugate(w));
     327        1407 :     PetscCall(VecScale(w,-1.0));
     328        1407 :     PetscCall(VecNestSetSubVec(f,0,v));
     329        1407 :     PetscCall(VecNestSetSubVec(g,0,y));
     330        1407 :     PetscCall(STApply(eps->st,f,g));
     331             :     /* beta = sqrt(2*real(v'*y)); */
     332        1407 :     PetscCall(VecDot(v,y,&dot));
     333        1407 :     beta1[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
     334             :     /* V(:,j) = v/beta1; */
     335             :     /* HV(:,j) = y/beta1; */
     336        1407 :     PetscCall(VecScale(v,1.0/beta1[j]));
     337        1407 :     PetscCall(VecScale(y,1.0/beta1[j]));
     338        1407 :     PetscCall(BVRestoreColumn(V,j,&v));
     339        1407 :     PetscCall(BVRestoreColumn(HV,j,&y));
     340             : 
     341        1407 :     PetscCall(BVGetColumn(HV,j,&x));
     342        1407 :     PetscCall(BVGetColumn(U,j+1,&v));
     343        1407 :     PetscCall(BVGetColumn(HU,j+1,&y));
     344        1407 :     PetscCall(VecCopy(x,v));
     345        1407 :     PetscCall(BVRestoreColumn(HV,j,&x));
     346             :     /* v = Orthogonalize HV(:,j) */
     347        1407 :     PetscCall(OrthogonalizeVector_Gruning(v,U,V,HU,HV,j+1,&beta1[j],k,hwork,PETSC_TRUE,breakdown));
     348             :     /* y = Hmult(v,1) */
     349        1407 :     PetscCall(VecCopy(v,w));
     350        1407 :     PetscCall(VecConjugate(w));
     351        1407 :     PetscCall(VecNestSetSubVec(f,0,v));
     352        1407 :     PetscCall(VecNestSetSubVec(g,0,y));
     353        1407 :     PetscCall(STApply(eps->st,f,g));
     354             :     /* beta = sqrt(2*real(v'*y)); */
     355        1407 :     PetscCall(VecDot(v,y,&dot));
     356        1407 :     beta2[j] = PetscSqrtReal(PetscRealPart(2*dot)); //FIXME Check beta != 0?
     357             :     /* U(:,j) = v/beta2; */
     358             :     /* HU(:,j) = y/beta2; */
     359        1407 :     PetscCall(VecScale(v,1.0/beta2[j]));
     360        1407 :     PetscCall(VecScale(y,1.0/beta2[j]));
     361        1407 :     PetscCall(BVRestoreColumn(U,j+1,&v));
     362        1407 :     PetscCall(BVRestoreColumn(HU,j+1,&y));
     363             :   }
     364         255 :   if (4*m > 100) PetscCall(PetscFree(hwork));
     365         255 :   PetscCall(VecDestroy(&w));
     366         255 :   PetscCall(VecDestroy(&f));
     367         255 :   PetscCall(VecDestroy(&g));
     368         255 :   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          57 :   for (k=0; k<eps->nconv; k++) {
     386          50 :     PetscCall(BVGetColumn(U,k,&u1));
     387          50 :     PetscCall(BVGetColumn(V,k,&v1));
     388             :     /* x1 = u1 + v1 */
     389          50 :     PetscCall(VecAXPY(u1,1.0,v1));
     390             :     /* x2 = conj(u1) - conj(v1) = conj(u1 - v1) = conj((u1 + v1) - 2*v1) */
     391          50 :     PetscCall(VecAYPX(v1,-2.0,u1));
     392          50 :     PetscCall(VecConjugate(v1));
     393          50 :     PetscCall(BVRestoreColumn(U,k,&u1));
     394          50 :     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        1434 : static PetscErrorCode Orthog_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscScalar *h,PetscScalar *c,PetscBool *breakdown)
     404             : {
     405        1434 :   PetscInt          i;
     406        1434 :   Mat               MX,MY,MXl,MYl;
     407        1434 :   Vec               c1,c2,hxl,hyl,hz;
     408        1434 :   PetscScalar       *cx1,*cx2;
     409        1434 :   PetscMPIInt       len;
     410             : 
     411        1434 :   PetscFunctionBegin;
     412        1434 :   PetscCall(PetscMPIIntCast(j,&len));
     413        1434 :   PetscCall(BVSetActiveColumns(X,0,j));
     414        1434 :   PetscCall(BVSetActiveColumns(Y,0,j));
     415             :   /* BVTDotVec does not exist yet, implemented via MatMult operations */
     416        1434 :   PetscCall(BVGetMat(X,&MX));
     417        1434 :   PetscCall(BVGetMat(Y,&MY));
     418        1434 :   PetscCall(MatDenseGetLocalMatrix(MX,&MXl));
     419        1434 :   PetscCall(MatDenseGetLocalMatrix(MY,&MYl));
     420        1434 :   PetscCall(MatCreateVecs(MXl,&c1,&hyl));
     421        1434 :   PetscCall(MatCreateVecs(MXl,&c2,&hxl));
     422             : 
     423             :   /* c1 =  X^* hx - Y^* hy
     424             :    * c2 = -Y^T hx + X^T hy */
     425             : 
     426        1434 :   PetscCall(VecGetLocalVector(hx,hxl));
     427        1434 :   PetscCall(VecGetLocalVector(hy,hyl));
     428        1434 :   PetscCall(VecDuplicate(hx,&hz));
     429             :   /* c1 = -(Y^* hy) */
     430        1434 :   PetscCall(MatMultHermitianTranspose(MYl,hyl,c1));
     431        1434 :   PetscCall(VecScale(c1,-1));
     432             :   /* c1 = c1 + X^* hx */
     433        1434 :   PetscCall(MatMultHermitianTransposeAdd(MXl,hxl,c1,c1));
     434        1434 :   PetscCall(VecGetArray(c1,&cx1));
     435        4302 :   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx1,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
     436        1434 :   PetscCall(VecRestoreArray(c1,&cx1));
     437             :   /* c2 = -(Y^T hx) */
     438        1434 :   PetscCall(MatMultTranspose(MYl,hxl,c2));
     439        1434 :   PetscCall(VecScale(c2,-1));
     440             :   /* c2 = c2 + X^T hy */
     441        1434 :   PetscCall(MatMultTransposeAdd(MXl,hyl,c2,c2));
     442        1434 :   PetscCall(VecGetArray(c2,&cx2));
     443        4302 :   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,cx2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)hx)));
     444        1434 :   PetscCall(VecRestoreArray(c2,&cx2));
     445        1434 :   PetscCall(VecRestoreLocalVector(hx,hxl));
     446        1434 :   PetscCall(VecRestoreLocalVector(hy,hyl));
     447             : 
     448             :   /* accumulate orthog coeffs into h */
     449        1434 :   PetscCall(VecGetArrayRead(c1,(const PetscScalar**)&cx1));
     450        1434 :   PetscCall(VecGetArrayRead(c2,(const PetscScalar**)&cx2));
     451       19952 :   for (i=0; i<j; i++) h[i] += cx1[i];
     452       19952 :   for (i=0; i<j; i++) h[i+j] += cx2[i];
     453        1434 :   PetscCall(VecRestoreArrayRead(c1,(const PetscScalar**)&cx1));
     454        1434 :   PetscCall(VecRestoreArrayRead(c2,(const PetscScalar**)&cx2));
     455             : 
     456             :   /* u = hx - X c1 - conj(Y) c2 */
     457             : 
     458             :   /* conj(Y) c2 */
     459        1434 :   PetscCall(VecConjugate(c2));
     460        1434 :   PetscCall(VecGetLocalVector(hz,hxl));
     461        1434 :   PetscCall(MatMult(MYl,c2,hxl));
     462        1434 :   PetscCall(VecConjugate(hxl));
     463             :   /* X c1 */
     464        1434 :   PetscCall(MatMultAdd(MXl,c1,hxl,hxl));
     465        1434 :   PetscCall(VecRestoreLocalVector(hz,hxl));
     466        1434 :   PetscCall(VecAXPY(hx,-1,hz));
     467             : 
     468        1434 :   PetscCall(BVRestoreMat(X,&MX));
     469        1434 :   PetscCall(BVRestoreMat(Y,&MY));
     470        1434 :   PetscCall(VecDestroy(&c1));
     471        1434 :   PetscCall(VecDestroy(&c2));
     472        1434 :   PetscCall(VecDestroy(&hxl));
     473        1434 :   PetscCall(VecDestroy(&hyl));
     474        1434 :   PetscCall(VecDestroy(&hz));
     475        1434 :   PetscFunctionReturn(PETSC_SUCCESS);
     476             : }
     477             : 
     478             : /* Orthogonalize vector x against first j vectors in X and Y */
     479        1434 : static PetscErrorCode OrthogonalizeVector_ProjectedBSE(Vec hx,Vec hy,BV X,BV Y,PetscInt j,PetscReal *beta,PetscInt k,PetscScalar *h,PetscBool *breakdown)
     480             : {
     481        1434 :   PetscInt    l,i;
     482        1434 :   PetscScalar alpha,alpha1,alpha2;
     483        1434 :   Vec         x,y;
     484             : 
     485        1434 :   PetscFunctionBegin;
     486        1434 :   PetscCall(PetscArrayzero(h,2*j));
     487             : 
     488             :   /* Local orthogonalization */
     489        1434 :   l = j==k+1?0:j-2;  /* 1st column to orthogonalize against */
     490             :   /* alpha = X(:,j-1)'*hx-Y(:,j-1)'*hy */
     491        1434 :   PetscCall(BVGetColumn(X,j-1,&x));
     492        1434 :   PetscCall(BVGetColumn(Y,j-1,&y));
     493        1434 :   PetscCall(VecDotBegin(hx,x,&alpha1));
     494        1434 :   PetscCall(VecDotBegin(hy,y,&alpha2));
     495        1434 :   PetscCall(VecDotEnd(hx,x,&alpha1));
     496        1434 :   PetscCall(VecDotEnd(hy,y,&alpha2));
     497        1434 :   PetscCall(BVRestoreColumn(X,j-1,&x));
     498        1434 :   PetscCall(BVRestoreColumn(Y,j-1,&y));
     499        1434 :   alpha = alpha1-alpha2;
     500             :   /* Store coeffs into h */
     501        5007 :   for (i=l; i<j-1; i++) h[i] = h[j+i] = beta[i];
     502        1434 :   h[j-1] = alpha;
     503        1434 :   h[2*j-1] = alpha-1.0;
     504             :   /* Orthogonalize: hx = hx - X(:,l:j-1)*h1 - conj(Y(:,l:j-1))*h2 */
     505        1434 :   PetscCall(BVSetActiveColumns(X,l,j));
     506        1434 :   PetscCall(BVSetActiveColumns(Y,l,j));
     507        1434 :   PetscCall(BVMultVec(X,-1.0,1.0,hx,h+l));
     508        1434 :   PetscCall(VecConjugate(hx));
     509        6441 :   for (i=j+l; i<2*j; i++) h[j+i] = PetscConj(h[i]);
     510        1434 :   PetscCall(BVMultVec(Y,-1.0,1.0,hx,h+2*j+l));
     511        1434 :   PetscCall(VecConjugate(hx));
     512             : 
     513             :   /* Full orthogonalization */
     514        1434 :   PetscCall(VecCopy(hx,hy));
     515        1434 :   PetscCall(VecConjugate(hy));
     516        1434 :   PetscCall(Orthog_ProjectedBSE(hx,hy,X,Y,j,h,h+2*j,breakdown));
     517        1434 :   PetscFunctionReturn(PETSC_SUCCESS);
     518             : }
     519             : 
     520         259 : static PetscErrorCode EPSBSELanczos_ProjectedBSE(EPS eps,BV X,BV Y,Vec v,PetscReal *alpha,PetscReal *beta,PetscInt k,PetscInt *M,PetscBool *breakdown)
     521             : {
     522         259 :   PetscInt       j,m = *M;
     523         259 :   Vec            u,x,y,w,f,g,vecs[2];
     524         259 :   Mat            H;
     525         259 :   IS             is[2];
     526         259 :   PetscReal      nrm;
     527         259 :   PetscScalar    *hwork,lhwork[100],gamma;
     528             : 
     529         259 :   PetscFunctionBegin;
     530         259 :   if (4*m > 100) PetscCall(PetscMalloc1(4*m,&hwork));
     531         259 :   else hwork = lhwork;
     532         259 :   PetscCall(STGetMatrix(eps->st,0,&H));
     533         259 :   PetscCall(MatNestGetISs(H,is,NULL));
     534             : 
     535             :   /* create work vectors */
     536         259 :   PetscCall(BVGetColumn(Y,0,&u));
     537         259 :   PetscCall(VecDuplicate(u,&w));
     538         259 :   vecs[0] = u;
     539         259 :   vecs[1] = w;
     540         259 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&f));
     541         259 :   PetscCall(VecCreateNest(PetscObjectComm((PetscObject)eps),2,is,vecs,&g));
     542         259 :   PetscCall(BVRestoreColumn(Y,0,&u));
     543             : 
     544             :   /* Normalize initial vector */
     545         259 :   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        1693 :   for (j=k;j<m;j++) {
     573             :     /* j+1 columns (indexes 0 to j) have been computed */
     574        1434 :     PetscCall(BVGetColumn(X,j+1,&x));
     575        1434 :     PetscCall(BVGetColumn(Y,j+1,&y));
     576             :     /* u = Hmult(v,-1)*/
     577        1434 :     PetscCall(VecCopy(v,w));
     578        1434 :     PetscCall(VecConjugate(w));
     579        1434 :     PetscCall(VecScale(w,-1.0));
     580        1434 :     PetscCall(VecNestSetSubVec(f,0,v));
     581        1434 :     PetscCall(VecNestSetSubVec(g,0,x));
     582        1434 :     PetscCall(STApply(eps->st,f,g));
     583             :     /* hx = (u+v) */
     584        1434 :     PetscCall(VecCopy(x,y));
     585        1434 :     PetscCall(VecAXPY(x,1,v));
     586             :     /* hy = conj(u-v) */
     587        1434 :     PetscCall(VecAXPY(y,-1,v));
     588        1434 :     PetscCall(VecConjugate(y));
     589             :     /* [u,cd] = orthog(hx,hy,X(:,1:j),Y(:,1:j),opt)*/
     590        1434 :     PetscCall(OrthogonalizeVector_ProjectedBSE(x,y,X,Y,j+1,beta,k,hwork,breakdown));
     591             :     /* alpha(j) = real(cd(j))-1/2 */
     592        1434 :     alpha[j] = 2*(PetscRealPart(hwork[j]) - 0.5);
     593             :     /* v = Hmult(u,1) */
     594        1434 :     PetscCall(VecCopy(x,w));
     595        1434 :     PetscCall(VecConjugate(w));
     596        1434 :     PetscCall(VecNestSetSubVec(f,0,x));
     597        1434 :     PetscCall(VecNestSetSubVec(g,0,y));
     598        1434 :     PetscCall(STApply(eps->st,f,g));
     599             :     /* nrm = sqrt(real(u'*v)) */
     600             :     /* beta(j) = nrm */
     601        1434 :     PetscCall(VecDot(x,y,&gamma));
     602        1434 :     beta[j] = 2.0*PetscSqrtReal(PetscRealPart(gamma));
     603             :     /* u = u/(nrm*2) */
     604        1434 :     PetscCall(VecScale(x,1.0/beta[j]));
     605             :     /* v = v/(nrm*2) */
     606        1434 :     PetscCall(VecScale(y,1.0/beta[j]));
     607        1434 :     PetscCall(VecCopy(y,v));
     608             :     /* X(:,j+1) = (u+v) */
     609        1434 :     PetscCall(VecAXPY(x,1,y));
     610             :     /* Y(:,j+1) = conj(u-v) */
     611        1434 :     PetscCall(VecAYPX(y,-2,x));
     612        1434 :     PetscCall(VecConjugate(y));
     613             :     /* Restore */
     614        1434 :     PetscCall(BVRestoreColumn(X,j+1,&x));
     615        1434 :     PetscCall(BVRestoreColumn(Y,j+1,&y));
     616             :   }
     617         259 :   if (4*m > 100) PetscCall(PetscFree(hwork));
     618         259 :   PetscCall(VecDestroy(&w));
     619         259 :   PetscCall(VecDestroy(&f));
     620         259 :   PetscCall(VecDestroy(&g));
     621         259 :   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          52 :   for (k=0; k<eps->nconv; k++) {
     640             :     /* approx eigenvector is [ delta1*x1 + delta2*conj(y1) ]
     641             :                              [ delta1*y1 + delta2*conj(x1) ] */
     642          46 :     lambda = eps->eigr[k];
     643          46 :     PetscCall(STBackTransform(eps->st,1,&lambda,&eps->eigi[k]));
     644          46 :     delta1 = lambda+1.0;
     645          46 :     delta2 = lambda-1.0;
     646          46 :     PetscCall(BVGetColumn(X,k,&x1));
     647          46 :     PetscCall(BVGetColumn(Y,k,&y1));
     648          46 :     PetscCall(VecCopy(x1,cx1));
     649          46 :     PetscCall(VecCopy(y1,cy1));
     650          46 :     PetscCall(VecConjugate(cx1));
     651          46 :     PetscCall(VecConjugate(cy1));
     652          46 :     PetscCall(VecScale(x1,delta1));
     653          46 :     PetscCall(VecScale(y1,delta1));
     654          46 :     PetscCall(VecAXPY(x1,delta2,cy1));
     655          46 :     PetscCall(VecAXPY(y1,delta2,cx1));
     656          46 :     PetscCall(BVRestoreColumn(X,k,&x1));
     657          46 :     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         260 :     eps->its++;
     752             : 
     753             :     /* Compute an nv-step Lanczos factorization */
     754         260 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     755         260 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     756         260 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
     757         260 :     b = a + ld;
     758         260 :     PetscCall(EPSBSELanczos_Shao(eps,U,V,a,b,eps->nconv+l,&nv,&breakdown));
     759         260 :     beta = b[nv-1];
     760         260 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
     761         260 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     762         260 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
     763         260 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     764             : 
     765             :     /* Solve projected problem */
     766         260 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     767         260 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     768         260 :     PetscCall(DSUpdateExtraRow(eps->ds));
     769         260 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     770             : 
     771             :     /* Check convergence */
     772        4006 :     for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
     773         260 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
     774         260 :     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
     775         260 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     776         260 :     nconv = k;
     777             : 
     778             :     /* Update l */
     779         260 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
     780         253 :     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
     781         260 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     782         260 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     783             : 
     784         260 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
     785         253 :       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         253 :       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         260 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
     792         260 :     PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
     793         260 :     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
     794         260 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
     795             : 
     796         260 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
     797         253 :       PetscCall(BVCopyColumn(eps->V,nv,k+l));
     798         253 :       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
     799           2 :         eps->ncv = eps->mpd+k;
     800           2 :         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     801           2 :         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
     802           2 :         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     803          14 :         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
     804           2 :         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
     805           2 :         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     806             :       }
     807             :     }
     808         260 :     eps->nconv = k;
     809         267 :     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         255 : static PetscErrorCode EPSConvergence_Gruning(EPS eps,PetscBool getall,PetscInt kini,PetscInt nits,PetscInt *kout)
     823             : {
     824         255 :   PetscInt       k,marker,ld;
     825         255 :   PetscReal      *alpha,*beta,resnorm;
     826         255 :   PetscBool      extra;
     827             : 
     828         255 :   PetscFunctionBegin;
     829         255 :   *kout = 0;
     830         255 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     831         255 :   PetscCall(DSGetExtraRow(eps->ds,&extra));
     832         255 :   PetscCheck(extra,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Only implemented for DS with extra row");
     833         255 :   marker = -1;
     834         255 :   if (eps->trackall) getall = PETSC_TRUE;
     835         255 :   PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&alpha));
     836         255 :   beta = alpha + ld;
     837         305 :   for (k=kini;k<kini+nits;k++) {
     838         305 :     resnorm = PetscAbsReal(beta[k]);
     839         305 :     PetscCall((*eps->converged)(eps,eps->eigr[k],eps->eigi[k],resnorm,&eps->errest[k],eps->convergedctx));
     840         305 :     if (marker==-1 && eps->errest[k] >= eps->tol) marker = k;
     841         305 :     if (marker!=-1 && !getall) break;
     842             :   }
     843         255 :   PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&alpha));
     844         255 :   if (marker!=-1) k = marker;
     845         255 :   *kout = k;
     846         255 :   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         255 :     eps->its++;
     880             : 
     881             :     /* Compute an nv-step Lanczos factorization */
     882         255 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     883         255 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     884         255 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&d1));
     885         255 :     d2 = d1 + ld;
     886         255 :     PetscCall(EPSBSELanczos_Gruning(eps,U,V,HU,HV,d1,d2,eps->nconv+l,&nv,&breakdown));
     887         255 :     beta = d1[nv-1];
     888         255 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&d1));
     889             : 
     890             :     /* Compute SVD */
     891         255 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     892         255 :     PetscCall(DSSVDSetDimensions(eps->ds,nv));
     893         255 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
     894         255 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     895             : 
     896         255 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     897         255 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     898         255 :     PetscCall(DSUpdateExtraRow(eps->ds));
     899         255 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     900             : 
     901             :     /* Check convergence */
     902         255 :     PetscCall(EPSConvergence_Gruning(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,&k));
     903         255 :     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
     904         255 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     905         255 :     nconv = k;
     906             : 
     907             :     /* Update l */
     908         255 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
     909         248 :     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
     910         255 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     911         255 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     912             : 
     913         255 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
     914         248 :       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         248 :       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         255 :     PetscCall(DSGetMat(eps->ds,DS_MAT_U,&Q));
     921         255 :     PetscCall(DSGetMat(eps->ds,DS_MAT_V,&Z));
     922         255 :     PetscCall(BVMultInPlace(U,Z,eps->nconv,k+l));
     923         255 :     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
     924         255 :     PetscCall(BVMultInPlace(HU,Z,eps->nconv,k+l));
     925         255 :     PetscCall(BVMultInPlace(HV,Q,eps->nconv,k+l));
     926         255 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_U,&Q));
     927         255 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_V,&Z));
     928             : 
     929         255 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
     930         248 :       PetscCall(BVCopyColumn(U,nv,k+l));
     931         248 :       PetscCall(BVCopyColumn(HU,nv,k+l));
     932         248 :       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
     933           2 :         eps->ncv = eps->mpd+k;
     934           2 :         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
     935           2 :         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
     936           2 :         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
     937           2 :         PetscCall(BVResize(HU,eps->ncv+1,PETSC_TRUE));
     938           2 :         PetscCall(BVResize(HV,eps->ncv+1,PETSC_TRUE));
     939          14 :         for (i=nv;i<eps->ncv;i++) { eps->perm[i] = i; eps->eigi[i] = 0.0; }
     940           2 :         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
     941           2 :         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     942             :       }
     943             :     }
     944         255 :     eps->nconv = k;
     945         262 :     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         259 :     eps->its++;
     989             : 
     990             :     /* Compute an nv-step Lanczos factorization */
     991         259 :     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
     992         259 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     993         259 :     PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
     994         259 :     b = a + ld;
     995         259 :     PetscCall(EPSBSELanczos_ProjectedBSE(eps,U,V,v,a,b,eps->nconv+l,&nv,&breakdown));
     996         259 :     beta = b[nv-1];
     997         259 :     PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
     998         259 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
     999         259 :     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
    1000         259 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
    1001             : 
    1002             :     /* Solve projected problem */
    1003         259 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
    1004         259 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
    1005         259 :     PetscCall(DSUpdateExtraRow(eps->ds));
    1006         259 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
    1007             : 
    1008             :     /* Check convergence */
    1009        4091 :     for (i=0;i<nv;i++) eps->eigr[i] = PetscSqrtReal(PetscRealPart(eps->eigr[i]));
    1010         259 :     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
    1011         259 :     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
    1012         259 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
    1013         259 :     nconv = k;
    1014             : 
    1015             :     /* Update l */
    1016         259 :     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
    1017         253 :     else l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
    1018         259 :     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
    1019         259 :     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
    1020             : 
    1021         259 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
    1022         253 :       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         253 :       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         259 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
    1029         259 :     PetscCall(BVMultInPlace(U,Q,eps->nconv,k+l));
    1030         259 :     PetscCall(BVMultInPlace(V,Q,eps->nconv,k+l));
    1031         259 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
    1032             : 
    1033         259 :     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
    1034         253 :       PetscCall(BVCopyColumn(eps->V,nv,k+l));
    1035         253 :       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
    1036           2 :         eps->ncv = eps->mpd+k;
    1037           2 :         PetscCall(BVRestoreSplitRows(eps->V,is[0],is[1],&U,&V));
    1038           2 :         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
    1039           2 :         PetscCall(BVGetSplitRows(eps->V,is[0],is[1],&U,&V));
    1040          14 :         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
    1041           2 :         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
    1042           2 :         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
    1043             :       }
    1044             :     }
    1045         259 :     eps->nconv = k;
    1046         265 :     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             : }