function [x, phist] = CGLS (G, d, x_ini, x_bar, tol)

% ------------------------------------------------------------------------
% Conjugate Gradients with early termination
% >> adopted from A. Nemirovski code
%
% solves the problem min ||Gx - d||
%
% INPUT:
%   G     -- operator matrix
%   d     -- rhs (observed data)
%   x_ini -- initial point
%   x_bar -- true solution, used for relative accuracy calculations
%   tol.f -- residual
%   tol.G -- error in G
%   tol.d -- relative error in d
%
% OUTPUT:
%   x     -- regularized solution
%   phist -- history points [norm(x) norm(res) rel.accuracy]
%
% ------------------------------------------------------------------------
% Developed by Oleg Grodzevich as a part of Master of Mathematics Thesis,
% University of Waterloo, Combinatorics and Optimization department.
%
% E-mail: illinar@mindon.net
% ------------------------------------------------------------------------

% ------------------------------------------------------------------------
% Initialization
% ------------------------------------------------------------------------

x       = x_ini;                % solution
res_lsp = G*x-d;                % residual of LS problem
res_neq = G'*res_lsp;           % residual of Normal Equations problem
rho     = res_neq;
err_d   = norm(d)*tol.d;        % error in rhs
phist   = [];                   % points history
nx_bar  = norm(x_bar);

% info header
disp([' ']);
disp(['   +------------------------------------------------------------------------------------+']);
disp([sprintf('   | CGLS method: %4d x %4d matrix, relative rhs error: %12.4e %17c|', size(G,1), size(G,2), tol.d, ' ')]);
disp(['   +------+------------------------+------------------------+---------------------------+']);
disp(['   |   ## |                norm(x) |             norm(Gx-d) |     relative accuracy [%] |']);
disp(['   +------+------------------------+------------------------+---------------------------+']);

% ------------------------------------------------------------------------
% Loop
% ------------------------------------------------------------------------
for itcount = 1:4*size(G,2),

  s         = G*res_neq;
  s2        = full(s'*s);
  alpha     = full(res_lsp'*s)/s2;

  % update solution and residual
  x         = x-alpha*res_neq;
  res_lsp   = res_lsp-alpha*s;

  % keep track of norm of the residual and solution, also check the accuracy
  norm_res  = norm(res_lsp);
  norm_x    = norm(x);
  acc_x     = norm(x-x_bar)/nx_bar;

  % save current point
  pt        = [norm_x norm_res acc_x];
  phist     = [pt ; phist];

  % plot point on the L-curve
  hold on; loglog(norm_x,norm_res,'bo'); hold off;

  % information
  disp([sprintf('   | %4d | %22.4e | %22.4e | %25.2f |', itcount, pt(1), pt(2), pt(3)*100)]);
  disp(['   +------+------------------------+------------------------+---------------------------+']);
 
  % stopping criteria
  if norm_res <= tol.f*(tol.G*norm_x+err_d), break; end;

  % final update
  sgm       = (s'*G)';
  rho       = rho-alpha*sgm;
  beta      = full(rho'*sgm)/s2;
  res_neq   = beta*res_neq-rho;
end
