function [Z,lyap_residual] = LR_ADI(A,B,opts)
% function [Z,lyap_residual] = LR_ADI(A,B,opts)
% Function that approximately solves the Lyapunov equation
% 
%  A*X + X*A^T + B*B^T = 0,    A stable
%
% by the LR-ADI method.
%
% INPUT
% A: coefficient matrix
% B: low-rank factor of the rhs
% opts: structure containing different parameters
%       opts.maxit_j: maximum number of iterations for LR-ADI.
%       opts.p: vector containing the shifts to be used
%       opts.tol_outer: outer tolerance for the Lyapunov residual


if not(isfield(opts,'inner_solve'))
    opts.inner_solve='direct';
end
j_max=opts.maxit_j;
tol_outer=opts.tol_outer;
p=opts.p;

% initial residual norm     
nrmb=norm(B'*B,'fro');  
theta=nrmb;
[n,sh]=size(B);
I=speye(n);

res_inner=0;
u=0;
lyap_residual=zeros(1,j_max);
lyap_residual(1)=nrmb;

Ssol=zeros(n,(j_max+1)*sh);
Upsilon=B;
j=1;

% compute once and for all the LU factorization of A or a preconditioner
% for the iterative solution of A
if norm(A-A',1)<1e-14
    A_sym=1;
    if ~strcmp(opts.inner_solve,'direct')
        params_innersolver.solver='bcg';
        params_innersolver.tol_res=1e-10;
        params_innersolver.max_iterations = 500;
        params_innersolver.verbose=0;   
    end
else
    A_sym=0;
    if ~strcmp(opts.inner_solve,'direct')
        opts.inner_solve_flag='gmres';
        params_innersolver.solver='bgmres';
        [params_innersolver.Q2,params_innersolver.precond_right2]=schur(opts.P2);
        params_innersolver.tol_res=1e-10;
        params_innersolver.max_iterations = 500;   
        params_innersolver.verbose=0;
    end
end


while j<=j_max

    AA=A+p(j)*I;
        
    % set the
    if A_sym
        if strcmp(opts.inner_solve,'iterative')
            if not(isfield(opts,'P'))
               L = ichol(-AA);
            else
               L=chol(-(opts.P+p(j)*I),'lower');
            end
            % relaxed inner residual tolerance from Kuerschner and Freitag:
            % eq: 46b - 47
            u=u-2*real(p(j))*norm(Ssol(:,(j-1)*sh+1:j*sh),'fro')*res_inner;  
            params_innersolver.tol_res=(j*tol_outer/j_max-2*u)/(sqrt(theta)*4);
            params_innersolver.precond_left=L;
            params_innersolver.precond_right=L';       
       end
    else
        if strcmp(opts.inner_solve,'iterative')
            params_innersolver.p=p(j);
            params_innersolver.prec_shifted=1;
            % relaxed inner residual tolerance from Kuerschner and Freitag:
            % eq: 46b - 47
            u=u-2*real(p(j))*norm(Ssol(:,(j-1)*sh+1:j*sh),'fro')*res_inner;  
            params_innersolver.tol_res=(j*tol_outer/j_max-2*u)/(sqrt(theta)*4);       
        end
    end

    % solve the current shifted linear system
    if strcmp(opts.inner_solve,'direct')
        S_temp = AA\Upsilon;
    else
        [S_temp,out] = block_solve(-AA,-Upsilon,params_innersolver); 
        res_inner=out.res_norm(end);
    end

  
    Upsilon_temp=Upsilon-2*real(p(j))*S_temp;
    
    % avoid further computation for complex shifts
    if imag(p(j))~=0
        beta=real(p(j))/imag(p(j));
        S_temp1=sqrt(2*(beta^2+1))*imag(S_temp);
        Upsilon_temp=Upsilon-4*real(p(j))*(real(S_temp)+beta*imag(S_temp));
        S_temp=sqrt(2)*(real(S_temp)+beta*imag(S_temp));
        Ssol(:,(j-1)*sh+1:j*sh)=S_temp;
        j=j+1;
        Ssol(:,(j-1)*sh+1:j*sh)=S_temp1; 
    else
        Ssol(:,(j-1)*sh+1:j*sh)=S_temp;      
    end
          
    Upsilon=Upsilon_temp;
            
    lyap_residual(j)=norm(Upsilon'*Upsilon,'fro');
   
    % check if the Lyapunov residual norm is sufficiently small
    if lyap_residual(j)<= nrmb*tol_outer
      % Done
      lyap_residual=lyap_residual(1:j)/nrmb;
      Z=Ssol(:,1:j*sh)*(kron(sqrt(-2*diag(real(p(1:j)))),eye(sh)));
      return
    end  
    j=j+1;
end



end
