function [x,flag,relres,iter,resvec] = lsqr(A,b,tol,maxit,M1,M2,x0,varargin)
%LSQR    LSQR Implementation of Conjugate Gradients on the Normal Equations.
%   X = LSQR(A,B) attempts to solve the system of linear equations A*X=B
%   for X if A is consistent, otherwise it attempts to solve the least squares
%   solution X that minimizes norm(B-A*X).  The M-by-N coefficient matrix A
%   need not be square but the right hand side column vector B must have
%   length M.  A may be a function AFUN such that AFUN(X) returns A*X and
%   AFUN(X,'transp') returns A'*X.
%
%   LSQR(A,B,TOL) specifies the tolerance of the method.  If TOL is []
%   then LSQR uses the default, 1e-6.
%
%   LSQR(A,B,TOL,MAXIT) specifies the maximum number of iterations.  If
%   MAXIT is [] then LSQR uses the default, min([M,N,20]).
%
%   LSQR(A,B,TOL,MAXIT,M1) and LSQR(A,B,TOL,MAXIT,M1,M2) use N-by-N
%   preconditioner M or M = M1*M2 and effectively solve the system
%   A*inv(M)*Y = B for Y, where X = M*Y.  If M is [] then a preconditioner is
%   not applied.  M may be a function MFUN such that MFUN(X) returns M\X and
%   MFUN(X,'transp') returns M'\X.
%
%   LSQR(A,B,TOL,MAXIT,M1,M2,X0) specifies the N-by-1 initial guess.  If X0 is
%   [] then LSQR uses the default, an all zero vector.
%
%   LSQR(AFUN,B,TOL,MAXIT,M1FUN,M2FUN,X0,P1,P2,...) passes parameters P1,P2,...
%   to AFUN: AFUN(X,P1,P2,...) and AFUN(X,P1,P2,...,'transp') and
%   similarly to the preconditioner functions M1FUN and M2FUN.
%    
%   [X,FLAG] = LSQR(A,B,TOL,MAXIT,M1,M2,X0) also returns a convergence FLAG:
%    0 LSQR converged to the desired tolerance TOL within MAXIT iterations.
%    1 LSQR iterated MAXIT times but did not converge.
%    2 preconditioner M was ill-conditioned.
%    3 LSQR stagnated (two consecutive iterates were the same).
%    4 one of the scalar quantities calculated during LSQR became too
%      small or too large to continue computing.
%
%   [X,FLAG,RELRES] = LSQR(A,B,TOL,MAXIT,M1,M2,X0) also returns an estimate of
%   the relative residual NORM(B-A*X)/NORM(B).  If FLAG is 0, RELRES <= TOL.
%
%   [X,FLAG,RELRES,ITER] = LSQR(A,B,TOL,MAXIT,M1,M2,X0) also returns the
%   iteration number at which X was computed: 0 <= ITER <= MAXIT.
%
%   [X,FLAG,RELRES,ITER,RESVEC] = LSQR(A,B,TOL,MAXIT,M1,M2,X0) also returns a
%   vector of the residual norm estimates at each iteration, including
%   NORM(B-A*X0).
%   
%   Example:
%   n = 100; on = ones(n,1); A = spdiags([-2*on 4*on -on],-1:1,n,n);
%   b = sum(A,2); tol = 1e-8; maxit = 15;
%   M1 = spdiags([on/(-2) on],-1:0,n,n); M2 = spdiags([4*on -on],0:1,n,n);
%   x = lsqr(A,b,tol,maxit,M1,M2,[]);
%   Use this matrix-vector product function
%      function y = afun(x,n,transp_flag)
%      if (nargin > 2) & strcmp(transp_flag,'transp')
%         y = 4 * x;
%         y(1:n-1) = y(1:n-1) - 2 * x(2:n);
%         y(2:n) = y(2:n) - x(1:n-1);
%      else
%         y = 4 * x;
%         y(2:n) = y(2:n) - 2 * x(1:n-1);
%         y(1:n-1) = y(1:n-1) - x(2:n);
%      end
%
%   as input to LSQR
%      x1 = lsqr(@afun,b,tol,maxit,M1,M2,[],n);
%
%   See also BICG, BICGSTAB, CGS, GMRES, MINRES, PCG, QMR, SYMMLQ, @.

%   Penny Anderson, 1996.
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 1.3 $ $Date: 2000/06/02 17:59:24 $

% Check for an acceptable number of input arguments
if nargin < 2
    error('Not enough input arguments.');
end

% Determine whether A is a matrix, a string expression,
% the name of a function or an inline object.
[atype,afun,afcnstr] = iterchk(A);
if isequal(atype,'matrix')
    % Check matrix and right hand side vector inputs have appropriate sizes
    [m,n] = size(A);
    if ~isequal(size(b),[m,1])
        es = sprintf(['Right hand side must be a column vector of' ...
                ' length %d to match the coefficient matrix.'],m);
        error(es);
    end
else
    m = size(b,1);
    if (size(b,2) ~= 1)
        error('Right hand side must be a column vector.');
    end
    % How many columns does A have?  Same as entries of x = A'*b.
    x = iterapp(afun,atype,afcnstr,b,varargin{:},'transp');
    n = size(x,1);
end

% Assign default values to unspecified parameters
if nargin < 3 | isempty(tol)
    tol = 1e-6;
end
if nargin < 4 | isempty(maxit)
    maxit = min([m,n,20]);
end

% Check for all zero right hand side vector => all zero solution
n2b = norm(b);                      % Norm of rhs vector, b
if (n2b == 0)                       % if    rhs vector is all zeros
    x = zeros(n,1);                  % then  solution is all zeros
    flag = 0;                        % a valid solution has been obtained
    relres = 0;                      % the relative residual is actually 0/0
    iter = 0;                        % no iterations need be performed
    resvec = 0;                      % resvec(1) = norm(b-A*x) = norm(0)      
    if (nargout < 2)
        itermsg('lsqr',tol,maxit,0,flag,iter,NaN);
    end
    return
end

if ((nargin >= 5) & ~isempty(M1))
    existM1 = 1;
    [m1type,m1fun,m1fcnstr] = iterchk(M1);
    if isequal(m1type,'matrix')
        if ~isequal(size(M1),[n,n])
            es = sprintf(['Preconditioner must be a square matrix' ...
                    ' of size %d to match the problem size.'],n);
            error(es);
        end      
    end   
else
    existM1 = 0;
    m1type = 'matrix';
end

if ((nargin >= 6) & ~isempty(M2))
    existM2 = 1;
    [m2type,m2fun,m2fcnstr] = iterchk(M2);
    if isequal(m2type,'matrix')
        if ~isequal(size(M2),[n,n])
            es = sprintf(['Preconditioner must be a square matrix' ...
                    ' of size %d to match the problem size.'],n);
            error(es);
        end
    end
else
    existM2 = 0;
    m2type = 'matrix';
end

if ((nargin >= 7) & ~isempty(x0))
    if ~isequal(size(x0),[n,1])
        es = sprintf(['Initial guess must be a column vector of' ...
                ' length %d to match the problem size.'],n);
        error(es);
    else
        x = x0;
    end
else
    x = zeros(n,1);
end

if ((nargin > 7) & isequal(atype,'matrix') & ...
        isequal(m1type,'matrix') & isequal(m2type,'matrix'))
    error('Too many input arguments.');
end

% Set up for the method
flag = 1;
xmin = x;                          % Iterate which has minimal residual so far
imin = 0;                          % Iteration at which xmin was computed
tolb = tol * n2b;                  % Relative tolerance
if isequal(atype,'matrix')
   u = b - A * x;
else
   u = b - iterapp(afun,atype,afcnstr,x,varargin{:});
end
beta = norm(u);
normr = beta;                      % Norm of residual

if (normr <= tolb)                 % Initial guess is a good enough solution
    flag = 0;
    relres = normr / n2b;
    iter = 0;
    resvec = normr;
    if (nargout < 2)
        itermsg('lsqr',tol,maxit,0,flag,iter,relres);
    end
    return
end

u = u / beta;
c = 1;
s = 0;
phibar = beta;
d = zeros(n,1);
if isequal(atype,'matrix')
    v = A' * u;
else
    v = iterapp(afun,atype,afcnstr,u,varargin{:},'transp');
end
if existM1
    if isequal(m1type,'matrix')
        v = M1' \ v;
    else
        v = iterapp(m1fun,m1type,m1fcnstr,v,varargin{:},'transp');
    end
    if isinf(norm(v,inf))
        flag = 2;
        relres = normr/n2b;
        iter = 0;
        resvec = normr;
        if nargout < 2
            itermsg('lsqr',tol,maxit,0,flag,iter,relres);
        end
        return
    end
end
if existM2
    if isequal(m2type,'matrix')
        v = M2' \ v;
    else
        v = iterapp(m2fun,m2type,m2fcnstr,v,varargin{:},'transp');
    end
    if isinf(norm(v,inf))
        flag = 2;
        relres = normr/n2b;
        iter = 0;
        resvec = normr;
        if nargout < 2
            itermsg('lsqr',tol,maxit,0,flag,iter,relres);
        end
        return
    end
end
alpha = norm(v);
v = v / alpha;

resvec = zeros(maxit+1,1);         % Preallocate vector for norm of residuals
resvec(1) = normr;                 % resvec(1) = norm(b-A*x0)
normrmin = normr;                  % Norm of residual from xmin

stag = 0;                          % stagnation of the method

% loop over maxit iterations (unless convergence or failure)

for i = 1 : maxit
    if existM2
        if isequal(m2type,'matrix')
            z = M2 \ v;
        else
            z = iterapp(m2fun,m2type,m2fcnstr,v,varargin{:});
        end
        if isinf(norm(z,inf))
            flag = 2;
            break
        end
    else
        z = v;
    end
    if existM1
        if isequal(m1type,'matrix')
            z = M1 \ z;
        else
            z = iterapp(m1fun,m1type,m1fcnstr,z,varargin{:});
        end
        if isinf(norm(z,inf))
            flag = 2;
            break
        end
    end
    if isequal(atype,'matrix')
        u = A * z - alpha * u;
    else
        u = iterapp(afun,atype,afcnstr,z,varargin{:}) - alpha * u;
    end
    beta = norm(u);
    u = u / beta;
    thet = - s * alpha;
    rhot = c * alpha;
    rho = sqrt(rhot^2 + beta^2);
    c = rhot / rho;
    s = - beta / rho;
    phi = c * phibar;
    if (phi == 0)                  % stagnation of the method
        stag = 1;
    end
    phibar = s * phibar;    
    d = (z - thet * d) / rho;
    
    % Check for stagnation of the method
    if (stag == 0)
        stagtest = zeros(n,1);
        ind = (x ~= 0);
        stagtest(ind) = d(ind) ./ x(ind);
        stagtest(~ind & d ~= 0) = Inf;
        if norm(stagtest,inf) < eps
            stag = 1;
        end
    end

    if normr <= tolb                 % check for convergence
        flag = 0;
        iter = i;
        break
    end
    
    if stag == 1
        flag = 3;
        break
    end
    
    if normr < normrmin              % update minimal norm quantities
        normrmin = normr;
        xmin = x;
        imin = i;
    end
    
    x = x + phi * d;
    normr = abs(s) * normr;
    resvec(i+1) = normr;
    if isequal(atype,'matrix')
        vt = A' * u;
    else
        vt = iterapp(afun,atype,afcnstr,u,varargin{:},'transp');
    end
    if existM1
        if isequal(m1type,'matrix')
            vt = M1' \ vt;
        else
            vt = iterapp(m1fun,m1type,m1fcnstr,vt,varargin{:},'transp');
        end
        if isinf(norm(vt,inf))
            flag = 2;
            break
        end
    end
    if existM2
        if isequal(m2type,'matrix')
            vt = M2' \ vt;
        else
            vt = iterapp(m2fun,m2type,m2fcnstr,vt,varargin{:},'transp');
        end
        if isinf(norm(vt,inf))
            flag = 2;
            break
        end
    end
    v = vt - beta * v;
    alpha = norm(v);
    v = v / alpha;
    
end                                % for i = 1 : maxit

% returned solution is first with minimal residual
if flag == 0
    relres = normr / n2b;
else
    x = xmin;
    iter = imin;
    relres = normrmin / n2b;
end

% truncate the zeros from resvec
if flag <= 1 | flag == 3
    resvec = resvec(1:i+1);
else
    resvec = resvec(1:i);
end

% only display a message if the output flag is not used
if nargout < 2
    itermsg('lsqr',tol,maxit,i,flag,iter,relres);
end
