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