function [lam, Z, res, hsty, info]=...
         xeLOBP4dCG(K, M, Z0, V, nrmKM, Nd, shift, tol, Ne, opts)
%
% xeLOBP4dCG: Expertly extended Locally Optimal Block Preconditioned 
% 4d Conjugate Gradient
%
% This function computes the (k0+1)st to (k0+Nd)th smallest positive 
% eigenvalues of
%
%    H=[ 0  K ], where K and M are n-by-n and positive definite and
%      [ M  0 ]  one of them is definite,
%
% where k0 is the number of columns in V which contains the y-components 
% of the k0 approximate eigenvectors associated with the first k0 nonegative 
% eigenvalues of H. 
%
% The eigenvalue problem is the same as A - lambda B, where
%
%         [ M  0 ]       [ 0  I ]
%  A =    [ 0  K ] , B = [ I  0 ]
%
%
%  Algorithm description can be found in the following reference
%
%  Zhaojun Bai, Ren-cang Li and Wenwei Lin,  Linear response
%  eigenvalue problem solved by extended locally optimal 
%  preconditioned conjugate gradient methods, Science China 
%  Mathematics Vol.59, pp.1443-1460, 2016.
%
%
% Inputs
% ======
%
%       K      array (n-by-n), SPD
%       M      array (n-by-n), SPD
%       Z0     array (2n-by-nb) whose columns span an approximate 
%              invariant subspace associated with (k0+1)st to (k0+k)th 
%              smallest positive eigenvalues
%       V      n-by-k0, whose columns contains the y-components of 
%              the k0 approximate eigenvectors associated with the first 
%              k0 nonegative eigenvalues of H. 
%              Possibly k0=0, i.e., U is empty array.
%       nrmKM  2-by-1, nrmKM(1): estimate of ||K||
%                      nrmKM(2): estimate of ||M||
%       Nd     integer, # of eigenpairs needed to be computed by this call
%       shift  real for shifting away known eigenvalues
%       tol    tolerance for testing convergence 
%       Ne     integer, the order of the Krylov subspace for building 
%              the extended search space Ne>=2. 
%              set Ne=2 for the usual xLOBP4dCG
%       opts   options ...
%              opts.precond=1,  use 1st type preconditioner 
%                               [A+[0 0; 0 shift*V*V']]^{-1}
%                   1) opts.precond_one=1: direct solver
%                   2) opts.precond_one=2: 
%                       [A+[0 0; 0 shift*V*V']]^{-1} implemented by CG 
%                   3) opts.precond_one=3: 
%                       use incomplete Cholesky of K and M to construct
%                       a preconditioner for [A+[0 0; 0 shift*V*V']]^{-1} 
%                       implemented by CG -- not implemented yet
%              opts.precond=2,  use 2nd type preconditioner; suppose A-sgm B 
%                               approx LDL'.
%                   1) opts.precond_two=1: (LL')^{-1} -- not implemented yet
%                   2) opts.precond_two=2: [LL'+[0 0; 0 shift*V*V']]^{-1}  
%                                         -- not implemented yet
%              opts.precond=3,  use 3rd type preconditioner 
%                        [A+[0 0; 0 shift*V*V']]-sgm B is SPD use linear CG
%
%              opts.precond_sgm  the shift parameter to build a preconditioner
%
% Outputs
% =======  
%
%        lam   computed eigenvalue
%        Z     corresponding eigenvectors
%        res   residual error  ||Hz - lam z||
%        hsty  struct 
%              hsty.eig -- row i contains the eigenvalue approximations 
%                          of ith iteration
%              hsty.res -- row i contains normalized residuals 
%                          ||H zi - lam_i zi||/(||H||*||z_i||+|lam_i|*||zi||)
%                          for ith iteration
%        info  struct
%              info.itn number of iterations
%              info.kc  number of converged eigenpairs by this call
%
% Note: The code will always run blocksize nb LOBP4dCG. Naturally 
% if one eigenpair is determined converged, we append a working block 
% by the next approximation in the process.
%
%---------------------------------------------------------------
% RCL and ZJB, ... 09/04/2013 --- 11/18/16
%---------------------------------------------------------------

n=max(size(K)); nb=size(Z0,2);  
k0=size(V,2);

kc=k0;  % # of converged eigenpairs (including pass-ins)
Z=[];

lam=[]; res=[];    % converged eigenvalues and corresponding residuals

maxitn=round(min(0.1*n,200)*(Nd/nb));

nrmH=max(nrmKM);

if opts.precond==1
   if opts.precond_one==1
      Rk=chol(K); Rm=chol(M);
   elseif opts.precond_one==2
      CGopts1.nitn=10;
      CGopts1.tol=1e-2;
      CGopts1.met=2;
      CGopts2.nitn=10;
      CGopts2.tol=1e-2;
      CGopts2.met=2;
      CG_M=0;
      CGx0=zeros(n,1);
      P=zeros(n,nb); Q=zeros(n,nb);
   elseif opts.precond_one==3
      % TBI
   end
elseif opts.precond==2
   sgm=opts.precond_sgm;
   % TBI
elseif opts.precond==3
   sgm=opts.precond_sgm;
   CGopts1.nitn=10;
   CGopts1.tol=1e-2;
   CGopts1.met=2;
%   CG_M=0;
   CGz0=zeros(2*n,1);
   Z=zeros(2*n,nb);
end  

[m,n]=size(K);
% initial approximations
X=Z0(n+1:2*n,:); 
Y=Z0(1:n,:);  
   

if k0==0, 

   if opts.precond==1
      if opts.precond_one==1
         KiV=[];
      elseif opts.precond_one==2
         CGopts1.update = 0;
         CGopts2.update = 0;
      elseif opts.precond_one==3
             % TBI
      end
   elseif opts.precond==2
             % TBI
   elseif opts.precond==3
      CGopts1.update = 0;
   end    
   
else
   if opts.precond==1
      if opts.precond_one==1
         % [K+shift*V*V']^{-1} = K^{-1}-shift*K^{-1}*V*[I+shift*V'*K^{-1}*V]^{-1}*V'*K^{-1}
         KiV=Rk\(Rk'\V);
         T=eye(k0)+shift*V'*KiV; invT=inv(T);
         KiViT=KiV*(shift*invT);
      elseif opts.precond_one==2
         CGopts1.update = 1;
         CGopts1.shift =shift;
         CGopts1.V = V;
         CGopts2.update = 0;
      elseif opts.precond_one==3
             % TBI
      end
   elseif opts.precond==2
             % TBI
   elseif opts.precond==3
       CGopts1.update = 1;
       CGopts1.shift =shift;
       CGopts1.V = V;
   end 
   
end

KX=K*X;
if k0>0,
   KX=KX+V*(shift*(V'*X)); 
end
XKX=X'*KX; colXKX=chol(XKX); invXKX=inv(colXKX); X=X*invXKX; KX=KX*invXKX; % X'*shiftedK*X=I
MY=M*Y;
YMY=Y'*MY; colYMY=chol(YMY); invYMY=inv(colYMY); Y=Y*invYMY; MY=MY*invYMY; % Y'*M*Y=I

W=X'*Y; W=0.5*(W+W');
RX=KX*W-Y; RY=MY*W-X;
   
% Preconditioning
if k0 == 0,

   if opts.precond==1
      if opts.precond_one==1
         % P=K\RX, Q=M\RY
         P=Rk\(Rk'\RX); Q=Rm\(Rm'\RY);
      elseif opts.precond_one==2
         % P=K\RX, Q=M\RY
         for i=1:nb,
             [P(:,i), error, iter, flag] = LinCG(K, CGx0, RX(:,i), CG_M, CGopts1); 
             [Q(:,i), error, iter, flag] = LinCG(M, CGx0, RY(:,i), CG_M, CGopts2); 
         end
      elseif opts.precond_one==3
             % TBI
      end
   elseif opts.precond==2
             % TBI
   elseif opts.precond==3
      for i=1:nb,
          [Z(:,i), error, iter, flag] = LinCG4LR(K, M, sgm, CGz0, [RY(:,i); RX(:,i)], CGopts1); 
          P=Z(n+1:2*n,:); Q=Z(1:n,:);
      end
   end 

else

   if opts.precond==1
      if opts.precond_one==1
         % P=[K+shift*V*V']\RX, Q=M\RY         
         KinvRX=Rk\(Rk'\RX);
         P=KinvRX-KiViT*(KiV'*RX);
         Q=Rm\(Rm'\RY);
      elseif opts.precond_one==2
         % P=[K+shift*V*V']\RX, Q=M\RY
         for i=1:nb,
             [P(:,i), error, iter, flag] = LinCG(K, CGx0, RX(:,i), CG_M, CGopts1); 
             [Q(:,i), error, iter, flag] = LinCG(M, CGx0, RY(:,i), CG_M, CGopts2); 
         end
      elseif opts.precond_one==3
             % TBI
      end
   elseif opts.precond==2
             % TBI
   elseif opts.precond==3
      for i=1:nb,
          [Z(:,i), error, iter, flag] = LinCG4LR(K, M, sgm, CGz0, [RY(:,i); RX(:,i)], CGopts1); 
          P=Z(n+1:2*n,:); Q=Z(1:n,:);
      end
   end 
   
end
   
% expand X; note ||X||_K=1
Xw=MGSgS(K, shift, V, P, X); Xw=[X, Xw];

% expand Y; note ||Y||_M=1
Yw=MGSg(M, Q, Y); Yw=[Y, Yw];

[U1,S1,V1] = svd(Yw'*Xw); [s,idx]=sort(diag(S1),'descend'); ns=length(s); %s=s(1:nb);
X=Xw*V1(:,idx(1:nb)); nX=size(Xw,2); XmX=Xw(:,nb+1:nX)*V1(nb+1:nX,idx(1:nb));
Y=Yw*U1(:,idx(1:nb)); nY=size(Yw,2); YmY=Yw(:,nb+1:nY)*U1(nb+1:nY,idx(1:nb));
rho=ones(1,ns)./(s(1:ns))';

KX=K*X; 
if k0>0,
   KX=KX+V*(shift*(V'*X));
end
MY=M*Y;

RX=KX-Y*diag(rho(1:nb)); RY=MY-X*diag(rho(1:nb));      % Residuals

nrmZ=sum(abs(X))+sum(abs(Y)); 
nrmR=(sum(abs(RX))+sum(abs(RY)))./(nrmH*nrmZ+rho(1:nb).*nrmZ);   % Residual errors in all-components
%hsty.eig=[hsty.eig; rho]; hsty.res=[hsty.res; nrmR];
hsty.eig=rho(1:nb); hsty.res=nrmR;

% Convergence test
cvindx=find(nrmR<=tol); kk=length(cvindx); 
cvindx_c=find(nrmR>tol);
% # of converged eigenpairs is kk.  Theorectically, all nb pairs could coverge at once. 
% But that's unlikely, we assume they don't, i.e., kk<nb.
if kk>0
   lam=[lam; (rho(cvindx))']; res=[res nrmR(cvindx)]; 
   Z=[Z, [Y(:,cvindx); X(:,cvindx)]];
   V=[V, Y(:,cvindx)]; kc=kc+kk;    
      
   if opts.precond==1
      if opts.precond_one==1
         % [K+shift*V*V']^{-1} = K^{-1}-shift*K^{-1}*V*[I+shift*V'*K^{-1}*V]^{-1}*V'*K^{-1}
         KiV=[KiV, Rk\(Rk'\Y(:,cvindx))];
         T=eye(kc)+shift*V'*KiV; invT=inv(T);
         KiViT=KiV*(shift*invT);
      elseif opts.precond_one==2
         CGopts1.update = 1;
         CGopts1.shift =shift;
         CGopts1.V = V;
      elseif opts.precond_one==3
             % TBI
      end
   elseif opts.precond==2
             % TBI
   elseif opts.precond==3
       CGopts1.update = 1;
       CGopts1.shift =shift;
       CGopts1.V = V;
   end 
   
   % make working block to nb again
   Xa=Xw*V1(:,idx(nb+1:nb+kk));           % Xa should have K-orthonormal columns
   Ya=Yw*U1(:,idx(nb+1:nb+kk));           % Ya should have M-orthonormal columns
   KXa=K*Xa+V*(shift*(V'*Xa));
   MYa=M*Ya;
   RXa=KXa-Ya*diag(rho(nb+1:nb+kk));
   RYa=MYa-Xa*diag(rho(nb+1:nb+kk));
   
   RX=[RX(:,cvindx_c) RXa]; RY=[RY(:,cvindx_c) RYa];
   rho=rho([cvindx_c nb+1:nb+kk]); 
   X=[X(:,cvindx_c) Xa]; Y=[Y(:,cvindx_c) Ya];
   KX=[KX(:,cvindx_c) KXa];               % Part of this KX does not include the updates for the ones that have just converged.
                                          % Probably doesn't matter because, they will be included in the
                                          % very next iteration and forward.
   MY=[MY(:,cvindx_c) MYa];
else
   rho=rho(1:nb);
end

itn=1;
   
while itn < maxitn && kc<k0+Nd,
    
    XKX=X'*KX; % It should be I, but as precaution ...
    colXKX=chol(XKX); invXKX=inv(colXKX); X=X*invXKX;  % X'*shiftedK*X=I
 
    YMY=Y'*MY; % It should be I, but as precaution ...
    colYMY=chol(YMY); invYMY=inv(colYMY); Y=Y*invYMY;  % Y'*M*Y=I
    
    % Build search space. Usually Ne is very small, 2, 3, 4, or 5. So we
    % will use brute force computations
    Xw=[]; Yw=[];
    for ie=1:Ne-1       

        % Preconditioning  
        if kc == 0,            
           if opts.precond==1    
              if opts.precond_one==1
                 % P=K\RX, Q=M\RY
                 P=Rk\(Rk'\RX); Q=Rm\(Rm'\RY);
              elseif opts.precond_one==2
                 % P=K\RX, Q=M\RY
                 for i=1:nb,
                     [P(:,i), error, iter, flag] = LinCG(K, CGx0, RX(:,i), CG_M, CGopts1); 
                     [Q(:,i), error, iter, flag] = LinCG(M, CGx0, RY(:,i), CG_M, CGopts2); 
                 end
              elseif opts.precond_one==3
                     % TBI
              end
           elseif opts.precond==2
                     % TBI
           elseif opts.precond==3
              for i=1:nb,
                  [Z(:,i), error, iter, flag] = LinCG4LR(K, M, sgm, CGz0, [RY(:,i); RX(:,i)], CGopts1); 
                  P=Z(n+1:2*n,:); Q=Z(1:n,:);
              end
           end               
        else % kc>0          
           if opts.precond==1
              if opts.precond_one==1
                 % P=[K+shift*V*V']\RX, Q=M\RY         
                 KinvRX=Rk\(Rk'\RX);
                 P=KinvRX-KiViT*(KiV'*RX);
                 Q=Rm\(Rm'\RY);
              elseif opts.precond_one==2
                 % P=[K+shift*V*V']\RX, Q=M\RY
                 for i=1:nb,
                     [P(:,i), error, iter, flag] = LinCG(K, CGx0, RX(:,i), CG_M, CGopts1); 
                     [Q(:,i), error, iter, flag] = LinCG(M, CGx0, RY(:,i), CG_M, CGopts2); 
                 end
              elseif opts.precond_one==3
                    % TBI
              end
           elseif opts.precond==2
                    % TBI
           elseif opts.precond==3
              for i=1:nb,
                  [Z(:,i), error, iter, flag] = LinCG4LR(K, M, sgm, CGz0, [RY(:,i); RX(:,i)], CGopts1); 
                  P=Z(n+1:2*n,:); Q=Z(1:n,:);
              end
           end             
        end           
           
        Xw=[Xw P]; Yw=[Yw Q];
        if ie<Ne-1
           KX=K*P;    
           if kc>0,
              KX=KX+V*(shift*(V'*P));
           end
           MY=M*Q;   
           RX=KX-Q*diag(rho); RY=MY-P*diag(rho);
        end
         
    end        
           
    Xn=MGSgS(K, shift, V, [Xw, XmX], X); Xw=[X, Xn];
    Yn=MGSg(M, [Yw, YmY], Y);            Yw=[Y, Yn];
    
    [U1,S1,V1] = svd(Yw'*Xw); [s,idx]=sort(diag(S1),'descend'); ns=length(s); %s=s(1:nb);
    X=Xw*V1(:,idx(1:nb)); nX=size(Xw,2); XmX=Xw(:,nb+1:nX)*V1(nb+1:nX,idx(1:nb));
    Y=Yw*U1(:,idx(1:nb)); nY=size(Yw,2); YmY=Yw(:,nb+1:nY)*U1(nb+1:nY,idx(1:nb));
    rho=ones(1,ns)./(s(1:ns))';
    
    KX=K*X;  
    if kc>0,
       KX=KX+V*(shift*(V'*X));
    end
    MY=M*Y;
    
    RX=KX-Y*diag(rho(1:nb)); RY=MY-X*diag(rho(1:nb));      % Residuals
    
    nrmZ=sum(abs(X))+sum(abs(Y)); 
    nrmR=(sum(abs(RX))+sum(abs(RY)))./(nrmH*nrmZ+rho(1:nb).*nrmZ);   % Residual errors in all-components
    hsty.eig=[hsty.eig; rho(1:nb)]; hsty.res=[hsty.res; nrmR];   
    
    % Convergence test
    cvindx=find(nrmR<=tol); kk=length(cvindx); 
    cvindx_c=find(nrmR>tol);
    % # of converged eigenpairs is kk.  Theorectically, all nb pairs could coverge at once. 
    % But that's unlikely, we assume they don't, i.e., kk<nb.
    if kk>0
       lam=[lam; (rho(cvindx))']; res=[res nrmR(cvindx)]; 
       Z=[Z, [Y(:,cvindx); X(:,cvindx)]];
       V=[V, Y(:,cvindx)]; kc=kc+kk;    
          
       if opts.precond==1
          if opts.precond_one==1
             % [K+shift*V*V']^{-1} = K^{-1}-shift*K^{-1}*V*[I+shift*V'*K^{-1}*V]^{-1}*V'*K^{-1}
             KiV=[KiV, Rk\(Rk'\Y(:,cvindx))];
             T=eye(kc)+shift*V'*KiV; invT=inv(T);
             KiViT=KiV*(shift*invT);
          elseif opts.precond_one==2
             CGopts1.update = 1;
             CGopts1.shift =shift;
             CGopts1.V = V;
          elseif opts.precond_one==3
                 % TBI
          end
       elseif opts.precond==2
                 % TBI
       elseif opts.precond==3
           CGopts1.update = 1;
           CGopts1.shift =shift;
           CGopts1.V = V;
       end 
       
       % make working block to nb again
       Xa=Xw*V1(:,idx(nb+1:nb+kk));           % Xa should have K-orthonormal columns
       Ya=Yw*U1(:,idx(nb+1:nb+kk));           % Ya should have M-orthonormal columns
       KXa=K*Xa+V*(shift*(V'*Xa));
       MYa=M*Ya;
       RXa=KXa-Ya*diag(rho(nb+1:nb+kk));
       RYa=MYa-Xa*diag(rho(nb+1:nb+kk));
       
       RX=[RX(:,cvindx_c) RXa]; RY=[RY(:,cvindx_c) RYa];
       %rho=rho(kk+1:kk+nb); 
       rho=rho([cvindx_c nb+1:nb+kk]); 
       X=[X(:,cvindx_c) Xa]; Y=[Y(:,cvindx_c) Ya];
       KX=[KX(:,cvindx_c) KXa];               % Part of this KX does not include the updates for the ones that have just converged.
                                              % Probably doesn't matter because, they will be included in the
                                              % very next iteration and forward.
       MY=[MY(:,cvindx_c) MYa];
    else
       rho=rho(1:nb);
    end
    
    itn=itn+1; 
  
end

info.itn = itn; info.kc=kc-k0; Z=[Z [Y; X]]; 
