function G = lr_ye(A, tol, maxiter)
%
%	solve the rate matrix G and U by logrithmic reduction alg.
% 	using refined inverting algorithm 
%	by Qiang Ye, Apr. 1999 
%
%	A0, A1, A2 - input matrices. 
% 	iter_max - max. iteration 
%
%	U, G -  the output
% 	iter - iteration counts
%	error - inf norm of error 
%

A2 = A(:,:,1);
A1 = A(:,:,2);
A0 = A(:,:,3);

n = size(A0, 1); 
e = ones(n,1); 
error=ones(maxiter,3); 
G=A0; 
H=G; L=G;
U=G;  			% initialization 

U=A1; 
v=(A0+A2)*e; 
M=bgth(U, e, v, [A0, A2]);
H=M(:,1:n); L=M(:,(n+1):(2*n)); 
G=L;
T=H;

for iter = 1:maxiter, 
    iter
  error(iter,2)=norm((H+L)*e-e,inf); 
  U = H*L + L*H; 
  H=H^2;
  L=L^2; 

  v=(H+L)*e; 
  M=bgth(U, e, v, [H, L]); 
  H=M(:,1:n); L=M(:,(n+1):(2*n)); 

  U=T*L; 
  error1 = norm(U, inf);  
  G=G+U;
  error(iter,1) = norm(e-G*e, inf); 
  error(iter,3) = norm(A2+A1*G+A0*(G^2)-G, inf);
  res = error(iter,3)
  T=T*H; 
%  error1=norm(T, inf); 
%                to be comparable with standard impl, use norm(U) instead

  inc=min(min(U));
  if(inc < 0), 
    fprintf('nonmonotonic at iteration') 
    iter
    inc
  end

  if (error1 <= tol), 
     break
  end
end 

