function [Z1,Z2,time_tot,time_res]=SKSM_Sylv_cTri(A,B,C1,C2,m,tol)
%[Z1,Z2,time_tot,time_res]=SKSM_Sylv_cTri(A,B,C1,C2,m,tol)
% Solution of the Sylvester equation 
%
% A*X+X*B+C1*C2^T=0,    A,B sym neg.def. and symmetric, C1,C2 tall
%
% by the (standard) block Krylov subspace method (SKSM) with Galerkin condition,
% with enhanced computational strategy for the residual norm computation.
%
% Reference manuscript:
%
% Davide Palitta and Valeria Simoncini
% Computationally enhanced projection methods for symmetric Sylvester and Lyapunov matrix equations
% February 2016, ArXiv: 1602.05033
%
% email: davide.palitta3@unibo.it, valeria.simoncini@unibo.it
%
% REQUIREMENT: the subroutines eig_dstevr and tridiag_dsbtrd which make use of LAPACK and C-BLAS subroutines. 
% 
%
% INPUT
% A, B: sym negative def. matrix  
% C1, C2: rhs matrix 
% m: max number of iterations
% tol: final residual accuracy (relative residual F-norm)
%
% OUTPUT
% Z1, Z2: tall matrices such that Z1*Z2^T is the approx solution (Z is truncated to sing.val >1e-12)
% time_tot: total cputime
% time_res: total time requested to compute the residual norm
%
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS 
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR 
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER 
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN 
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 

% Check the input
if norm(A-A','fro')>1e-12 
    error('The coefficient matrix must be symmetric')
end

if norm(B-B','fro')>1e-12 
    error('The coefficient matrix must be symmetric')
end

if size(A)~=size(B)
    error('The current implementation deals with the case size(A)=size(B) only')
end

tic;
[n,s]=size(C1); 
nrmb=norm(C1,'fro')*norm(C2,'fro');

% Orthonormalize C1
[U(1:n,1:s),beta]=qr(C1,0);

% Orthonormalize C2
[V(1:n,1:s),betaB]=qr(C2,0);

beta2=beta*betaB';

% Preallocation
TA=zeros((m+1)*s,m*s);
TTA=zeros(s+1,m*s);
TB=zeros((m+1)*s,m*s);
TTB=zeros(s+1,m*s);

time_res=zeros(m,1);
time_subtot=0;
        
t1=toc;
for j=1:m
    tic
    jms=(j-1)*s+1;j1s=(j+1)*s;js=j*s;js1=js+1;
   
    %sequence in A
    %new basis block (double modified gram-schmidt)
    Up = A*U(1:n,end-s+1:end);
  
    for l=1:2
        K=min([j, 3]);
        for kk=1:K
              k1=(kk-1)*s+1; k2=kk*s;
              coeff= U(1:n,k1:k2)'*Up;
              TA((j-K+kk-1)*s+1:(j-K+kk)*s,jms:js) = TA((j-K+kk-1)*s+1:(j-K+kk)*s,jms:js)+ coeff; 
              Up = Up - U(1:n,k1:k2)*coeff; 
        end
    end
    
  
    if (j<2)
        [U(1:n,s+1:2*s),TA(js1:j1s,jms:js)]=qr(Up,0);
        
        % prepare LLA for dstrd
        for i=1:s
            TTA(end-i+1,i:s)=diag(TA(jms:js,jms:js),i-1);
        end
        LLA=TTA(:,1:j*s);
        
    else
        % move the blocks
        if j>2
            U(1:n,1:2*s)=U(1:n,s+1:3*s);
        end
        [U(1:n,2*s+1:3*s),TA(js1:j1s,jms:js)]=qr(Up,0);
    
        % prepare LLA for dstrd
        HHA=zeros(s+1,s);
        for i=1:s+1
            HHA(end-i+1,1:s)=diag(TA(js-s+2-i:js-i+1,js-s+1:js));
        end
        TTA(:,(j-1)*s+1:j*s)=HHA;
        LLA=TTA(:,1:j*s);
    
    end
    
    %sequence in B
    %new basis block (double modified gram-schmidt)
    Vp = B*V(1:n,end-s+1:end);
  
    for l=1:2
        K=min([j, 3]);
        for kk=1:K
              k1=(kk-1)*s+1; k2=kk*s;
              coeff= V(1:n,k1:k2)'*Vp;
              TB((j-K+kk-1)*s+1:(j-K+kk)*s,jms:js) = TB((j-K+kk-1)*s+1:(j-K+kk)*s,jms:js)+ coeff; 
              Vp = Vp - V(1:n,k1:k2)*coeff; 
        end
    end
    
  
    if (j<2)
        [V(1:n,s+1:2*s),TB(js1:j1s,jms:js)]=qr(Vp,0);
        
        % prepare LLB for dstrd
        for i=1:s
            TTB(end-i+1,i:s)=diag(TB(jms:js,jms:js),i-1);
        end
        LLB=TTB(:,1:j*s);
        
    else
        % move the blocks
        if j>2
            V(1:n,1:2*s)=V(1:n,s+1:3*s);
        end
        [V(1:n,2*s+1:3*s),TB(js1:j1s,jms:js)]=qr(Vp,0);
    
        % prepare LLA for dstrd
        HHB=zeros(s+1,s);
        for i=1:s+1
            HHB(end-i+1,1:s)=diag(TB(js-s+2-i:js-i+1,js-s+1:js));
        end
                
        TTB(:,(j-1)*s+1:j*s)=HHB;
        LLB=TTB(:,1:j*s);
    
    end
    
    t2=toc;
        
    tic
    % take the ((j+1)s,js) block of TA and TB
    ll = TA(js1:j1s,jms:js)';             
    ll2 = TB(js1:j1s,jms:js)';             

    % tridiagonalize LLA and LLb by dsbtrd
    [D,E,P]=tridiag_dsbtrd(LLA);
    [DB,EB,PB]=tridiag_dsbtrd(LLB);
  
    
    % compute the eigendecompositions of the obtained tridiagonal matrices by dstevr
    [Q,lambda]=eig_dstevr(D,E);
    [QB,lambdaB]=eig_dstevr(DB,EB);
    
    
    % Compute the residual norm
    S=((P(1:s,:)*Q)'*beta2*(PB(1:s,:)*QB))./(lambda*ones(1,js)+ones(js,1)*lambdaB');
    W1=(P(end-s+1:end,:)*Q)'*ll;
    W2=(PB(end-s+1:end,:)*QB)'*ll2;
    
    res=sqrt(norm(S*W2,'fro')^2+norm(S'*W1,'fro')^2);

    
    time_res(j)=toc;

    fprintf('It: %3d, Relative Res: %10.5e\n', j, res/nrmb)
    time_subtot=time_subtot+t2;
    
    if (res/nrmb<tol) 
         break
    end
    

end
tic
it=j;


% Compute the full eigendecomposition of T and the matrix Y
TA(1:js,1:js)=(TA(1:js,1:js)+TA(1:js,1:js)')/2;
TB(1:js,1:js)=(TB(1:js,1:js)+TB(1:js,1:js)')/2;
[QA,lambda]=eig(full(TA(1:js,1:js)));
lambdaA=diag(lambda);
Q1=QA(1:s,:);
[QB,lambda]=eig(full(TB(1:js,1:js)));
lambdaB=diag(lambda);
Q2=QB(1:s,:);

% Compute the solution of the last projected equation
Y=-(Q1'*beta2*Q2)./(lambdaA*ones(1,js)+ones(js,1)*lambdaB');

% Truncation strategy
[uY,sY,vY]=svd(full(Y)); [sY,id]=sort(abs(diag(sY)));
sY=flipud(sY);
uY=uY(:,id(end:-1:1));
vY=vY(:,id(end:-1:1));
is=sum(abs(sY)>1e-12);
Y01 = QA*(uY(:,1:is)*diag(sqrt(sY(1:is))));
Y02= QB*(vY(:,1:is)*diag(sqrt(sY(1:is))));
t4=toc;

% Two-pass

U=[];
V=[];
TB=zeros(j1s,js);

tic
[U(1:n,1:s),~]=qr(C1,0);
Z1=U(1:n,1:s)*Y01(1:s,:);

[V(1:n,1:s),~]=qr(C2,0);
Z2=V(1:n,1:s)*Y02(1:s,:);

%separe j=1,2 from the other j's to avoid "if" within the "for" loop
j=1;
j1s=2*s;
js1=s+1;
% sequence in A
UpA = A*U(1:n,end-s+1:end);
U(1:n,s+1:2*s)=(UpA-U(1:n,end-s+1:end)*TA((j-1)*s+1:j*s,(j-1)*s+1:j*s))/TA(j*s+1:(j+1)*s,(j-1)*s+1:j*s);
Z1=Z1+U(:,end-s+1:end)*Y01(js1:j1s,:);

j=2;
j1s=3*s;
js1=2*s+1;
%sequence in A
UpA = A*U(1:n,end-s+1:end);
U(1:n,2*s+1:3*s)=(UpA-U(1:n,1:2*s)*TA((j-2)*s+1:j*s,(j-1)*s+1:j*s))/TA(j*s+1:(j+1)*s,(j-1)*s+1:j*s);
Z1=Z1+U(:,end-s+1:end)*Y01(js1:j1s,:);

for j=3:it-1
    j1s=(j+1)*s;
    js1=j*s+1;
    
    % sequence in A
    UpA = A*U(1:n,end-s+1:end);
    U(1:n,1:2*s)=U(1:n,s+1:3*s);
    U(1:n,2*s+1:3*s)=(UpA-U(1:n,1:2*s)*TA((j-2)*s+1:j*s,(j-1)*s+1:j*s))/TA(j*s+1:(j+1)*s,(j-1)*s+1:j*s);
    Z1=Z1+U(:,end-s+1:end)*Y01(js1:j1s,:);
end
 % sequence in B
 for j=1:it-1
    % tic
     if j<=2
         K=j;
     else
         K=3;
     end
     jms=(j-1)*s+1;j1s=(j+1)*s;js=j*s;js1=js+1;
     
     
   % Sequence in B
     UpB = B*V(1:n,end-s+1:end);
     
      for l=1:2
         for kk=1:K
               k1=(kk-1)*s+1; k2=kk*s;
               gamma= V(1:n,k1:k2)'*UpB;
               TB((j-K+kk-1)*s+1:(j-K+kk)*s,jms:js) = TB((j-K+kk-1)*s+1:(j-K+kk)*s,jms:js)+ gamma; % U(1:n,k1:k2)'*Up;
               UpB = UpB - V(1:n,k1:k2)*gamma; %H(k1:k2,jms:js);
         end
     end
     
%   
     if (j<2)
         [V(1:n,s+1:2*s),TB(js1:j1s,jms:js)]=qr(UpB,0);
     else
         % moving the blocks
         if j>2
             V(1:n,1:2*s)=V(1:n,s+1:3*s);
         end
         [V(1:n,2*s+1:3*s),TB(js1:j1s,jms:js)]=qr(UpB,0);
     end
    Z2=Z2+V(:,end-s+1:end)*Y02(js1:j1s,:);
 end
time_twopass=toc;
fprintf('\n time_twopass: %10.5e\n',time_twopass)

time_subtot=time_subtot+t4+t1+time_twopass;
time_res=sum(time_res);
time_tot=time_subtot+time_res;    
