function [x, error, iter] = livolant(x, b, atv, errtol, maxit, varargin)
% Livolant-accelerated linear equation solver
% function [x, error, total_iters] = livolant(x0, b, atv, errtol, maxit, varargin)
% input parameters:
% x      initial iterate
% b      right hand side
% atv    character name of a matrix-vector product routine returning x+(b-Ax)
%        when x is input. The format for atv is "function x = atv(x,b,p1,p2,...)"
%        where p1 and p2 are optional parameters
% errtol relative residual reduction factor
% maxit  maximum number of iterations
% varargin optional parameters (p1,p2,...) for atv
% output parameters:
% x      solution of the linear system
% error  vector of residual norms for the history of the iteration
% iter   number of iterations
% (c) 2007 Alain Hebert, Ecole Polytechnique de Montreal
  icl1=3; icl2=3;
  n=length(b);
  errtol=errtol*norm(b);
  rho=Inf; error=[ ];
  old2=zeros(n);
  iter=0;
  while((rho > errtol) && (iter < maxit))
    iter=iter+1;
    old1=old2; old2=x;
    r=feval(atv,x,b,varargin{:})-x;
    rho=norm(r); error=[error;rho];
    x=x+r;
    if mod(iter-1,icl1+icl2) >= icl1
      r1=old2-old1; r2=x-old2;
      dmu=-sum(r1.*(r2-r1))/sum((r2-r1).*(r2-r1));
      if dmu > 0
        x=old2+dmu.*r2;
        old2=old1+dmu.*r1;
      end
    end
  end