minres1.m

function [x,flag,iter,relres,relAres,...
          Anorm,Acond,xnorm,Axnorm,resvec,Aresvec] = ...
	minres1(A,b,rtol,maxit,M,shift,maxxnorm,Acondlim,show)
%MINRES1: min-length solution to symmetric (possibly singular) Ax=b or min||Ax-b||.
%										%   This is the MINRES-QLP code with the QLP part ripped out.
%
%   X = MINRES1(A,B) solves the system of linear equations A*X=B
%   or the least-squares problem min norm(B-A*X) if A is singular.
%   The N-by-N matrix A must be symmetric or Hermitian, but need not be
%   positive definite or nonsingular.  It may be double or single.
%   The rhs vector B must have length N.  It may be real or complex,
%   double or single.
%
%   X = MINRES1(AFUN,B) accepts a function handle AFUN instead of
%   the matrix A.  Y = AFUN(X) returns the matrix-vector product Y=A*X.
%   In all of the following syntaxes, A can be replaced by AFUN.
%
%   X = MINRES1('AFUN',B) accepts the name of a function handle AFUN
%   instead of the matrix A.  Y = AFUN(X) returns the matrix-vector product
%   Y=A*X. In all of the following syntaxes, A can be replaced by AFUN.
%
%   X = MINRES1(A,B,RTOL) specifies a stopping tolerance.
%   If RTOL=[] or is absent, a default value is used.
%   (Similarly for all later input parameters.)
%   Default RTOL=eps.
%
%   X = MINRES1(A,B,RTOL,MAXIT)
%   specifies the maximum number of iterations.  Default MAXIT=4N.
%
%   X = MINRES1(A,B,RTOL,MAXIT,M)
%   uses a matrix M as preconditioner.  M must be positive definite
%   and symmetric or Hermitian.  It may be a function handle MFUN
%   such that Y=MFUN(X) returns Y=M\X.
%   If M=[], a preconditioner is not applied.
%
%   X = MINRES1(A,B,RTOL,MAXIT,M,SHIFT)
%   solves (A - SHIFT*I)X = B, or the corresponding least-squares problem
%   if (A - SHIFT*I) is singular, where SHIFT is a real scalar.
%   Default SHIFT=0.
%
%   X = MINRES1(A,B,RTOL,MAXIT,M,SHIFT,MAXXNORM,ACONDLIM)
%   specifies three parameters associated with singular or
%   ill-conditioned systems (A - SHIFT*I)*X = B.
%
%   MAXXNORM is an upper bound on NORM(X).
%   Default MAXXNORM=1e7.
%
%   ACONDLIM is an upper bound on ACOND, an estimate of COND(A).
%   Default ACONDLIM=1e15.
%
%   X = MINRES1(A,B,RTOL,MAXIT,M,SHIFT,MAXXNORM,ACONDLIM,SHOW)
%   specifies the printing option.
%   If SHOW=true,  an iteration log will be output.
%   If SHOW=false, the log is suppressed.
%   Default SHOW=true.
%
%
%   [X,FLAG] = MINRES1(A,B,...) returns a convergence FLAG:
%   -1 (beta2=0)  B and X are eigenvectors of (A - SHIFT*I).
%    0 (beta1=0)  B = 0.  The exact solution is X = 0.
%    1 X solves the compatible (possibly singular) system (A - SHIFT*I)X = B
%      to the desired tolerance:
%         RELRES = RNORM / (ANORM*XNORM + NORM(B)) <= RTOL,
%      where
%              R = B - (A - SHIFT*I)X and RNORM = norm(R).
%    2 X solves the incompatible (singular) system (A - SHIFT*I)X = B
%      to the desired tolerance:
%         RELARES = ARNORM / (ANORM * RNORM) <= RTOL,
%      where
%              AR = (A - SHIFT*I)R and ARNORM = NORM(AR).
%    3 Same as 1 with RTOL = EPS.
%    4 Same as 2 with RTOL = EPS.
%    5 X converged to an eigenvector of (A - SHIFT*I).
%    6 XNORM exceeded MAXXNORM.
%    7 ACOND exceeded ACONDLIM.
%    8 MAXIT iterations were performed before one of the previous
%      conditions was satisfied.
%    9 The system appears to be exactly singular.  XNORM does not
%      yet exceed MAXXNORM, but would if further iterations were
%      performed.
%
%    [X,FLAG,ITER] = MINRES1(A,B,...) returns ITER,  the
%    number of iterations performed.
%
%    [X,FLAG,ITER,RELRES,RELARES] = MINRES1(A,B,...)
%    returns relative residuals for (A - SHIFT*I)X = B and the
%    associated least-squares problem.  RELRES and RELARES are
%    defined above in the description of FLAG.
%
%    [X,FLAG,ITER,RELRES,RELARES,ANORM,ACOND,XNORM,AXNORM] =
%       MINRES1(A,B,...) returns
%       ANORM,  an estimate of the 2-norm of A-SHIFT*I.
%       ACOND,  an estimate of COND(A-SHIFT*I,2).
%       XNORM,  a recurred estimate of NORM(X).
%       AXNORM, a recurred estimate of NORM((A-SHIFT*I)X)
%
%    [X,FLAG,ITER,RELRES,RELARES,ANORM,ACOND,XNORM,AXNORM,...
%       RESVEC,ARESVEC] = MINRES1(A,B,...) returns
%       RESVEC,  a vector of estimates of NORM(R) at each iteration,
%                including NORM(B) as the first entry.
%       ARESVEC, a vector of estimates of NORM((A-SHIFT*I)R) at each
%                iteration, including NORM((A-SHIFT*I)B) as the first entry.
%       RESVEC and ARESVEC have length ITER+1.
%
%
% EXAMPLE 1:
%   n = 100;                                  e = ones(n,1);
%   A = spdiags([-2*e 4*e -2*e],-1:1,n,n);    M = spdiags(4*e,0,n,n);
%   b = sum(A,2);             rtol = 1e-10;   maxit = 50;
%   x = minres1(A,b,rtol,maxit,M);
%
%   Alternatively, use this matrix-vector product function:
%     function y = Afun(x,n)
%       y = 4*x;
%       y(2:n)   = y(2:n)   - 2*x(1:n-1);
%       y(1:n-1) = y(1:n-1) - 2*x(2:n);
%   as input to minres1:
%   n = 100;
%   A = @(x)Afun(x,n);
%   x = minres1(A,b,rtol,maxit,M);
%
% EXAMPLE 2: A is Laplacian on a 50 by 50 grid, singular and indefinite.
%   n = 50;  N = n^2;  e = ones(n,1);
%   B = spdiags([e e e], -1:1, n, n);
%   A = sparse([],[],[],N,N,(3*n-2)^2);
%   for i=1:n
%     A((i-1)*n+1:i*n,(i-1)*n+1:i*n) = B;
%     if i*n+1 < n*n,   A(i*n+1:(i+1)*n,(i-1)*n+1:i*n)     = B; end
%     if (i-2)*n+1 > 0, A((i-2)*n+1:(i-1)*n,(i-1)*n+1:i*n) = B; end
%   end
%   b = sum(A,2);   rtol   = 1e-5;    shift = 0;    maxxnorm = 1e2;
%   M = [];         Acondlim = [];    show  = true;
%   x = minres1(A,b,rtol,N,M,shift,maxxnorm,Acondlim,show);
%
% EXAMPLE 3: A is diagonal, singular and indefinite.
%   h = 1;  a = -10; b = -a; n = 2*b/h + 1;
%   A = spdiags((a:h:b)', 0, n, n);
%   b = ones(n,1);  rtol   = 1e-6;    shift = 0;    maxxnorm = 1e2;
%   M = [];         Acondlim = [];    show  = true;
%   x = minres1(A,b,rtol,N,M,shift,maxxnorm,Acondlim,show);
%
% EXAMPLE 4: Use this matrix-vector product function in file \Workouts\Afun.m:
%       %-----------------------------------%
%         function y = Afun(x, N)
%
%         n = length(x);
%
%         e = ones(N,1);
%         T = spdiags([e e e], -1:1, N, N);
%
%         M = n/N;
%         y = zeros(n,1);
%
%         for i=1:M
%           i1 = (i-2)*N + 1;
%           i2 = (i-1)*N + 1;
%           i3 = i    *N + 1;
%           i4 = i3 + N - 1;
%           if (i1 >=1 && i3 <= n)
%             y(i2:(i3-1)) = T * (x(i1:(i2-1)) + x(i2:(i3-1)) + x(i3:i4));
%           elseif (i1 <= 0 && i3 <= n)
%             y(i2:(i3-1)) = T * (               x(i2:(i3-1)) + x(i3:i4));
%           elseif (i1 >=1 && i3 > n)
%             y(i2:(i3-1)) = T * (x(i1:(i2-1)) + x(i2:(i3-1))            );
%           end
%         end
%       %-----------------------------------%
% Call MINRES1 using the anonymous function handle defined by Afun.m:
%       b = rand(400,1);
%       x = minres1(@(x,N)Afun(x,20), b);
%
% NOTE: The first x is vector output. The 'x''s in the first input are
% abstraction for the input of fun. Note how specific values for additional
% input parameters of Afun (such as N) is passed into the function.
%
% Alternatively, we can define the function handle
%       g = @(x,N)Afun(x,20);
% and then call MINRES1 using
%       x = minres1(g, b);
% or
%       x = minres1('@(x,N)Afun(x,20)', b);
%
% See also BICG, BICGSTAB, BICGSTABL, CGS, GMRES, LSQR, PCG, QMR, SYMMLQ,
% TFQMR, CHOLINC, FUNCTION_HANDLE.
% Also MINRES, SYMMLQ, LSQR, CGLS downloadable from
%   http://www.stanford.edu/group/SOL/software.html
%
% REFERENCES:
% Sou-Cheng T. Choi and Michael A. Saunders, 
% ALGORITHM: MINRES-QLP for Singular Symmetric and Hermitian Linear
% Equations and Least-Squares Problems, to appear in ACM Transactions on
% Mathematical Software.
%
% Sou-Cheng T. Choi, Christopher C. Paige, and Michael A. Saunders,
% MINRES-QLP: A Krylov Subspace Method for Indefinite or Singular Symmetric
% Systems, SIAM Journal on Scientific Computing, Vol. 33, No. 4, August
% 2011, pp. 1810--1836.
%
% Sou-Cheng T. Choi's PhD Dissertation, Stanford University, 2006:
%   http://www.stanford.edu/group/SOL/dissertations.html
%
% CURRENT / FUTURE RELEASES of minresqlp:
%   Version 2:  
%      http://code.google.com/p/minres-qlp/
%      http://www.mathworks.com/matlabcentral/fileexchange
%   Version 1:  
%      http://code.google.com/p/minres-qlp/
%      http://www.stanford.edu/group/SOL/download.html
%   Other implementations in Fortran 90/95, Python: 
%      http://code.google.com/p/minres-qlp/
%

% 
%
% MODIFICATION HISTORY:
%   28 Jun 2013: Second version for MATLAB Central.
%   27 Jun 2013: (1) Fixed iteration log:
%                (a) The heading came out ok every 20 lines initially, but
%                stopped after itn 80. 
%                (b) Subsequent output was for itns 109, 119, ... rather
%                than 110, 120,...
%                (2) Introduced local variable "likeLS", which is true if
%                Ax = b looks more like a least-squares system.
%                (3) Fixed an error in minresxxxM(). Moved two lines that
%                involve R into the else statement. Added code for handling
%                preconditioner as a function handle.
%                (4) Added debug statements.
%   28 Jul 2012: (1) Fixed a bug in pnorm_0 used in Anorm recurrence
%                relation so that {Anorm(k)} are monotonic increasing
%                underestimates of ||A||_2.
%                (2) Fixed SymOrtho to ensure the 2x2 Hermitian reflectors
%                are orthonormal.
%   20 Jul 2012: (1) Changed default RTOL to 1e-15 and default MAXIT to 4N.
%                (2) Initalized a few local variables to zeros.
%   12 May 2006: Created MINRESQLPs.m from research file minresqlp35.m.
%   19 Apr 2010: Help formatted like Matlab routines. Added output Aresvec.
%   25 Apr 2010: Eliminated unreferenced variables. Shortened
%                resvec and Aresvec from maxit to (iter+1).
%                The final rnorm, Arnorm, xnorm are computed directly.
%   26 Apr 2010: Special tests for beta2=0 and/or alfa1=0 when no
%                preconditioning.
%   02 May 2010: First version for
%                http://www.stanford.edu/group/SOL/software.html.
%
% KNOWN BUGS:
%   DD Mmm YYYY: ---
%
% NOTES:
%   DD Mmm YYYY: ---
%
% AUTHORS: Sou-Cheng (Terrya) Choi, CI, University of Chicago
%          Michael Saunders, SOL, Stanford University
%
% COPYRIGHT NOTICE:
%
%   This is Copyrighted Material. The software is COPYRIGHTED by the
%   original authors.
%
% COPYLEFT NOTICE:
%
%   Permission is granted to make and distribute verbatim copies of this
%   file, provided that the entire file is copied **together** as a
%   **unit**.
%
%   The purpose of this permission notice is to allow you to make copies
%   of the software and distribute them to others, for free or for a fee,
%   subject to the constraint that you maintain the entire file here
%   as a unit.  This enables people you give the software to be aware
%   of its origins, to ask questions of us by e-mail, to request
%   improvements, obtain later releases, and so forth.
%
%   If you seek permission to copy and distribute translations of this
%   software into another language, please e-mail a specific request to
%   saunders@stanford.edu and scchoi@stanford.edu.
%
%   If you seek permission to excerpt a **part** of the software library,
%   for example to appear in a scientific publication, please e-mail a
%   specific request to saunders@stanford.edu and scchoi@stanford.edu.
%
% COMMENTS?
%
%   Email sctchoi@uchicago.edu and saunders@stanford.edu 
%


%%  Check inputs and set default values.
debug  = true;
n      = length(b);
nin    = nargin;
precon = true;
if nin <  2 || ~logical(exist('A'       ,'var')) || isempty(A) ...
            || ~logical(exist('b'       ,'var')) || isempty(b)       ,...
                          error('Please provide at least non-empty input for A and b');  
end
if nin <  3 || ~logical(exist('rtol'    ,'var')) || isempty(rtol)    , rtol     = eps;   end
if nin <  4 || ~logical(exist('maxit'   ,'var')) || isempty(maxit)   , maxit    = 4*n;   end
if nin <  5 || ~logical(exist('M'       ,'var')) || isempty(M)       , precon   = false; end
if nin <  6 || ~logical(exist('shift'   ,'var')) || isempty(shift)   , shift    = 0;     end
if nin <  7 || ~logical(exist('maxxnorm','var')) || isempty(maxxnorm), maxxnorm = 1e7;   end
if nin <  8 || ~logical(exist('Acondlim','var')) || isempty(Acondlim), Acondlim = 1e15;  end
if nin <  9 || ~logical(exist('show'    ,'var')) || isempty(show)    , show     = true;  end
%if nin< 10 || ~logical(exist('disable' ,'var')) || isempty(disable) , disable  = false; end

if nargout> 11
  resvec  = zeros(maxit+1,1);
  Aresvec = zeros(maxit+1,1);
else
  resvec  = [];
  Aresvec = [];
end

if (ischar(A))
  if (~isempty(A))
    A = str2func(A);
  else
    error('Empty string for function name A');
  end
end
if (precon && ischar(M))
  if (~isempty(M))
    M = str2func(M);
  else
    error('Empty string for function name M');
  end
end

%% Set up {beta1, p, v} for the first Lanczos vector v1.
r2    = full(b);    % r2    = b
r3    = r2;         % r3    = b
beta1 = norm(r2);   % beta1 = norm(b)
if precon
  r3    = minresxxxM(M,r2);   % M*r3  = b
  beta1 = r3'*r2;             % beta1 = b'*inv(M)*b
  if beta1 < 0
    error('"M" appears to be indefinite.');
  else
    beta1 = sqrt(beta1);
  end
end

%% Initialize other quantities.
flag0    = -2;     flag     = flag0;
iter     = 0;
lines    = 1;      headlines= 20;
beta     = 0;      tau      = 0;          taul     = 0;      phi      = beta1;
betan    = beta1;  gmin     = 0;          cs       = -1;     sn       = 0;
cr1      = -1;     sr1      = 0;          cr2      = -1;     sr2      = 0;
dltan    = 0;      eplnn    = 0;          gama     = 0;      gamal    = 0;
eta      = 0;      etal     = 0;          etal2    = 0;  % BvdS: initial value of gamal2 overwritten 
vepln    = 0;      veplnl   = 0;          veplnl2  = 0;      ul3      = 0;
ul2      = 0;      ul       = 0;          rnorm    = betan;  % BvdS:  initial value of u overwritten     
xnorm    = 0;      xl2norm  = 0;          Axnorm   = 0;
Anorm    = 0;      Acond    = 1;
gminl    = 0;
relres   = rnorm / (beta1 + 1e-50);       % Safeguard for beta1 = 0
relresl  = 0;
relAresl = 0;
x        = zeros(n,1);
xl2      = x;
w        = x;
wl       = x;
% BvdS:  initial value of wl2 overwritten; wl2 may be section-local.
r1       = x;

if ~isempty(resvec)
  resvec(1)  = real(beta1);
end

%% print header if show
first = 'Enter minres1.  ';
last  = 'Exit minres1.  ';
msg=[' beta2 = 0.  b and x are eigenvectors                   '   % -1
     ' beta1 = 0.  The exact solution is  x = 0               '   %  0
     ' A solution to Ax = b found, given rtol                 '   %  1
     ' Min-length solution for singular LS problem, given rtol'   %  2
     ' A solution to Ax = b found, given eps                  '   %  3
     ' Min-length solution for singular LS problem, given eps '   %  4
     ' x has converged to an eigenvector                      '   %  5
     ' xnorm has exceeded maxxnorm                            '   %  6
     ' Acond has exceeded Acondlim                            '   %  7
     ' The iteration limit was reached                        '   %  8
     ' Least-squares problem but no converged solution yet    ']; %  9
head = '    iter     rnorm     Arnorm   Compatible     LS        Anorm      Acond      xnorm';
if show
  fprintf('\n%s%s', first, 'Min-length solution of symmetric (A-sI)x = b or min ||(A-sI)x - b||')
  fprintf('\nn      =%7g   ||b||    =%10.3e   shift    =%10.3e   rtol     =%10.3e',...
             n,            beta1,             shift,             rtol)
  fprintf('\nmaxit  =%7g   maxxnorm =%10.3e   Acondlim =%10.3e',...
             maxit,        maxxnorm,          Acondlim)
  fprintf('\nprecon =%7g\n' , precon)
  fprintf('\n%s\n', head)
end

if beta1==0, flag = 0; end   % b = 0 => x = 0.  We will skip the main loop.

%% Main iteration
while flag == flag0 && iter < maxit

%% Lanczos
  iter  = iter + 1;
  betal = beta;        beta = real(betan);
  v = r3*(1/beta);     r3 = minresxxxA(A,v);
  if shift ~= 0,       r3 = r3 - shift*v;         end
  if iter  >  1,       r3 = r3 - (beta/betal)*r1; end

  alfa = real(r3'*v);  %%%%% Allow for Hermitian A.  Must get real alfa here.
  r3   = r3 - (alfa/beta)*r2;     r1 = r2;     r2 = r3;

  if ~precon
    betan = norm(r3);
    if iter == 1       % Something special can happen
      if betan == 0    % beta2 = 0
        if alfa == 0   % alfa1 = 0
          flag = 0;    % Ab = 0 and x = 0         ("A" = (A - shift*I))
          break
        else
          flag = -1;   % Ab = alfa1 b,  x = b/alfa1, an eigenvector
          x    = full(b)/alfa;
          break
        end
      end
    end
  else
    r3 = minresxxxM(M,r2);     betan = r2'*r3;
    if betan > 0
      betan = sqrt(betan);
    else
      error('"M" appears to be indefinite or singular.');
    end
  end
  if iter <= 2
    pnorm  = norm([alfa betan]);
  else
    pnorm  = norm([betal alfa betan]);
  end

  if debug
    fprintf('\n\nLanczos iteration %d :\n', iter);
    fprintf('\n  v_%d     = ', iter );  fprintf('%s ', num2str(v(1:min(n,5))' ) );
    fprintf('\n  r1_%d    = ', iter );  fprintf('%s ', num2str(r1(1:min(n,5))') );
    fprintf('\n  r2_%d    = ', iter );  fprintf('%s ', num2str(r2(1:min(n,5))') );
    fprintf('\n  r3_%d    = ', iter );  fprintf('%s ', num2str(r3(1:min(n,5))') );
    fprintf('\n  alpha_%d = %s, beta_%d = %s, beta_%d = %s pnorm_%d = %s ',...
      iter, num2str(alfa), iter, num2str(beta), iter+1, num2str(betan), iter, num2str(pnorm) );
  end


%% Apply previous left reflection Q_{k-1}
  dbar  = dltan;
  dlta  = cs*dbar + sn*alfa;    epln     = eplnn;
  gbar  = sn*dbar - cs*alfa;    eplnn    = sn*betan;
  dltan = -cs*betan;            dlta_QLP = dlta;

  if debug
    fprintf('\n\nApply previous left reflection Q_{%d,%d}:\n', iter-1, iter');
    fprintf('\n  c_%d     = %s, s_%d    = %s', ...
      iter-1, num2str(cs), iter-1, num2str(sn) );
    fprintf('\n  dlta_%d = %s, gbar_%d = %s', ...
      iter, num2str(dlta), iter, num2str(gbar) );
    fprintf('\n  epln_%d = %s, dbar_%d = %s', ...
      iter+1, num2str(eplnn), iter+1, num2str(dltan) );
  end

%% Compute the current left reflection Q_k
  gamal2 = gamal;      gamal    = gama;  % BvdS:  gamal3 not used. 
  [cs,sn,gama] = SymOrtho(gbar, betan);    gama_tmp = gama;
  taul2  = taul;       taul   = tau;       tau      = cs*phi;
  Axnorm = norm([Axnorm tau]);             phi      = sn*phi;

  if debug
    fprintf('\n\nCompute the current left reflection Q_{%d,%d}:\n', iter, iter+1      );
    fprintf('\n  c_%d     = %s, s_%d    = %s ', iter, num2str(cs) , iter, num2str(sn) );
    fprintf('\n  tau_%d   = %s, phi_%d  = %s ', iter, num2str(tau), iter, num2str(phi));
    fprintf('\n  gama_%d = %s ', iter, num2str(gama)  );
  end

%% Apply the previous right reflection P{k-2,k}
  if iter > 2
    veplnl2  = veplnl;     etal2 = etal;     etal = eta;
    dlta_tmp = sr2*vepln - cr2*dlta;
    veplnl   = cr2*vepln + sr2*dlta;
    dlta     = dlta_tmp;   eta = sr2*gama;   gama = -cr2*gama;

  if debug   % BvdS:  apply condtional to debug print
    fprintf('\n\nApply the previous right reflection P_{%d,%d}:\n', iter-2, iter)
    fprintf('\n  cr2_%d   = %s, sr2_%d    = %s', ...
      iter, num2str(cr2), iter, num2str(sr2) );
    fprintf('\n  gama_%d = %s, gama_%d  = %s, gama_%d = %s', ...
      iter-2, num2str(gamal2), iter-1, num2str(gamal), iter, num2str(gama));
    fprintf('\n  dlta_%d = %s, vepln_%d = %s, eta_%d   = %s', ...
      iter, num2str(dlta), iter-1, num2str(veplnl), iter, num2str(eta) );
  end
  end

%% Compute the current right reflection P{k-1,k}, P_12, P_23,...
  if iter > 1
    [cr1, sr1, gamal] = SymOrtho(gamal, dlta);
    vepln =   sr1*gama;
    gama  = - cr1*gama;

    if debug
       fprintf('\n\nCompute the second current right reflections P_{%d,%d}:\n', iter-1, iter');
       fprintf('\n  cr1_%d   = %s, sr1_%d   = %s', ...
                    iter, num2str(cr1), iter, num2str(sr1) );
       fprintf('\n  gama_%d = %s, gama_%d = %s, vepln_%d = %s',...
                    iter-1, num2str(gamal), iter, num2str(gama), iter, num2str(vepln) );
     end
  end

%% Update xnorm
  xnorml = xnorm;     ul4 = ul3;     ul3   = ul2;
  if iter > 2
    ul2 = (taul2 - etal2*ul4 - veplnl2*ul3) / gamal2;
  end
  if iter > 1
    ul = ( taul  - etal *ul3 - veplnl *ul2) / gamal;
  end
  xnorm_tmp = norm([xl2norm ul2 ul]);
  likeLS =  (relresl >= relAresl);
  if abs(gama) > realmin && xnorm_tmp < maxxnorm
    u = (tau - eta*ul2 - vepln*ul) / gama;
    if norm([xnorm_tmp u]) > maxxnorm && likeLS
      u = 0;      flag = 6;
    end
  else
    u = 0;     flag = 9;
  end
  xl2norm = norm([xl2norm ul2]);
  xnorm   = norm([xl2norm ul u]);

%% Update w. Update x except if it will become too big  BvdS: fixed flag check
  wl2 = wl;     wl = w;
  w   = (v - epln*wl2 - dlta_QLP*wl) * (1/gama_tmp);
  if xnorm < maxxnorm
    x = x + tau*w;
  else
    flag = 6;
  end

  if debug
    fprintf('\n\nUpdate w:\n');
    fprintf('\n  w_%d     = ', iter-1 ); fprintf('%s ', num2str(wl(1:min(n,5))') );
    fprintf('\n  w_%d     = ', iter );   fprintf('%s ', num2str(w(1:min(n,5))') );
    fprintf('\n\nUpdate u, x and xnorm:\n');
    fprintf('\n  u_%d     = %s, u_%d     = %s, u_%d     = %s',...
      iter-2, num2str(ul2), iter-1, num2str(ul), iter, num2str(u) );
    fprintf('\n  x_%d     = ', iter ); fprintf('%s ', num2str(x(1:min(n,5))') );
    fprintf('\n  ||x_%d|| = ', iter ); fprintf('%s ', num2str(xnorm) );
  end

%% Compute the next right reflection P{k-1,k+1}
  gamal_tmp = gamal;
  [cr2,sr2,gamal] = SymOrtho(gamal,eplnn);

%% Estimate various norms
  abs_gama = abs(gama);      Anorml = Anorm;
  Anorm = max([Anorm, gamal, abs_gama, pnorm]);

  if iter == 1
    gmin   = gama;    gminl = gmin;
  elseif iter > 1
    gminl2 = gminl;   gminl = gmin;    gmin = min([gminl2, gamal, abs_gama]);
  end
  Acondl   = Acond;     Acond   = Anorm/gmin;
  rnorml   = rnorm;     relresl = relres;
  if (flag ~= 9)
    rnorm = phi;
  end
  relres   = rnorm / (Anorm*xnorm + beta1);
  rootl    = norm([gbar; dltan]);
  Arnorml  = rnorml*rootl;
  relAresl = rootl / Anorm;

%% See if any of the stopping criteria are satisfied.
  epsx = Anorm*xnorm*eps;
  if (flag == flag0) || (flag == 9)
    t1 = 1 + relres;
    t2 = 1 + relAresl;
    if iter     >= maxit   , flag = 8; end  % Too many itns
    if Acond    >= Acondlim, flag = 7; end  % Huge Acond
    if xnorm    >= maxxnorm, flag = 6; end  % xnorm exceeded its limit
    if epsx     >= beta1   , flag = 5; end  % x is an eigenvector
    if t2       <= 1       , flag = 4; end  % Accurate LS solution
    if t1       <= 1       , flag = 3; end  % Accurate Ax=b solution
    if relAresl <= rtol    , flag = 2; end  % Good enough LS solution
    if relres   <= rtol    , flag = 1; end  % Good enough Ax=b solution
  end

  if debug
    fprintf('\n\nUpdate other norms:\n');
    fprintf('\n  gmin_%d   = ', iter );   fprintf('%s ', num2str(gmin));
    fprintf('\n  pnorm_%d  = ', iter );   fprintf('%s ', num2str(pnorm));
    fprintf('\n  rnorm_%d  = ', iter );   fprintf('%s ', num2str(rnorm));
    fprintf('\n  Arnorm_%d = ', iter-1);  fprintf('%s ', num2str(Arnorml));
    fprintf('\n  Acond_%d  = ', iter );   fprintf('%s ', num2str(Acond));
    fprintf('\n\n');
  end


% The "disable" option allowed iterations to continue until xnorm
% became large and x was effectively a nullvector.
% We know that r will become a nullvector much sooner,
% so we now disable the disable option :)

%  if disable && (iter < maxit)
%    flag = 0;
%    if Axnorm < rtol*Anorm*xnorm
%      flag = 10;
%    end
%  end

  if flag == 2 || flag == 4 || (flag == 6 && likeLS) || flag == 7   % Possibly singular
    iter  = iter - 1;
    Acond = Acondl;   rnorm = rnorml;   relres = relresl;
  else
    if ~isempty(resvec)
      resvec(iter+1) = rnorm;
      Aresvec(iter)  = Arnorml;
    end

    if show && mod(iter-1,lines) == 0
      if iter == 101
        lines = 10;     headlines = 20*lines;
      elseif iter == 1001
        lines = 100;    headlines = 20*lines;
      end
      fprintf('%s', ' ')
      fprintf('%7g %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n',...
        iter-1, rnorml, Arnorml, relresl, relAresl, Anorml, Acondl, xnorml)
      if (iter > 1 && mod(iter,headlines) == 1)
        fprintf('\n%s\n', head)
      end
    end
  end
end % while


%% We have exited the main loop.
  fprintf('%s', ' ')

%% Compute final quantities directly.
r1      = b - minresxxxA(A,x) + shift*x;   % r1 is workspace for residual vector
rnorm   = norm(r1);
Arnorm  = norm(minresxxxA(A,r1) - shift*r1);
xnorm   = norm(x);
relres  = rnorm / (Anorm*xnorm + beta1);
relAres = 0;
if rnorm > realmin
  relAres = Arnorm / (Anorm*rnorm);
end

if ~isempty(Aresvec)
  Aresvec(iter+1) = Arnorm;
  Aresvec         = Aresvec(1:iter+1);
end
if ~isempty(resvec)
  resvec          = resvec(1:iter+1);
end

if show
  if rnorm > realmin
    fprintf('%7g %10.2e %10.2eD%10.2e %10.2eD%10.2e %10.2e %10.2e\n\n',...
            iter, rnorm, Arnorm, relres, relAres, Anorm, Acond, xnorm)
  else
    fprintf('%7g %10.2e %10.2eD%10.2e %10.2e             %10.2e %10.2e\n\n',...
            iter, rnorm, Arnorm, relres,          Anorm, Acond, xnorm)
  end
end

fprintf('\n%s flag  =%7g  %s'                                , last, flag , msg(flag+2,:))
fprintf('\n%s iter  =%7g'       , last, iter)
fprintf('\n%s rnorm = %11.4e     rnorm  direct = %11.4e'     , last, rnorm, norm(r1))
fprintf('\n%s                         Arnorm direct = %11.4e', last, Arnorm)
fprintf('\n%s xnorm = %11.4e     xnorm  direct = %11.4e'     , last, xnorm, norm(x))
fprintf('\n%s Anorm = %11.4e     Acond         = %11.4e\n'   , last, Anorm, Acond)


%% Private functions
function p = minresxxxA( A, x )
  if isa(A,'function_handle')
    p = A(x);
  else
    p = A*x;
  end
end % BvdS: added this

function p = minresxxxM( M, x )
  persistent R
  if isa(M,'function_handle')
    p = M(x);
  else  % BvdS:  Split up elseif into two separate lines
      if ~logical(exist('R','var')) || isempty(R)
        R = chol(M);
      end
    p = R'\x;
    p = R\p;
  end
end  % BvdS:  added this

%% end of minres1


function [c, s, r] = SymOrtho(a, b)
 
% SymOrtho: Stable Symmetric Householder reflection 
%  
%  USAGE:
%     [c, s, r] = SymOrtho(a, b)
%
%  INPUTS:
%    a      first element of a two-vector  [a; b]      
%    b      second element of a two-vector [a; b] 
%
%  OUTPUTS:
%    c      cosine(theta), where theta is the implicit angle of rotation
%           (counter-clockwise) in a plane-rotation
%    s      sine(theta)
%    r      two-norm of [a; b]
%  
%  DESCRIPTION:
%     Stable symmetric Householder reflection that gives c and s such that 
%        [ c  s ][a] = [d],
%        [ s -c ][b]   [0]  
%     where d = two-norm of vector [a, b],
%        c = a / sqrt(a^2 + b^2) = a / d, 
%        s = b / sqrt(a^2 + b^2) = b / d.
%     The implementation guards against overlow in computing sqrt(a^2 + b^2).
%
%  EXAMPLE:
%      description
%
%  SEE ALSO:
%     TESTSYMGIVENS.m, 
%     PLANEROT (MATLAB's function) --- 4 divisions while 2 would be enough, 
%     though not too time-consuming on modern machines
%  
%  REFERENCES:
%    Algorithm 4.9, stable *unsymmetric* Givens rotations in
%     Golub and van Loan's book Matrix Computations, 3rd edition. 
%
%  MODIFICATION HISTORY:
%    10/06/2004: Replace d = norm([a,b]) by 
%                        d = a/c if |b| < |a| or b/s otherwise.
%    10/07/2004: First two cases (either b or a == 0) rewritten to make sure 
%                (1) d >= 0 
%                (2) if [a,b] = 0, then c = 1 and s = 0 but not c = s = 0. 
%    09/27/2011: Change filename from SYMGIVENS2 to SYMORTHO.
%    01/16/2012: Change file from SYMORTHO to SYMREFL.
%     
%
%  KNOWN BUGS:
%     MM/DD/2004: description
%
%  AUTHORS: Sou-Cheng (Terrya) Choi, CI, University of Chicago
%          Michael Saunders, SOL, Stanford University
%
%  CREATION DATE: 09/28/2004

absa  = abs(a);
absb  = abs(b);
signa = sign(a);
signb = sign(b);
  
if isreal([a b])  
  %------------------------------
  % Both a and b are real numbers
  %------------------------------
     
  %...........................
  % Special cases: a or b is 0
  %...........................
  if b == 0
    if a == 0
      c = 1; 
    else 
      c = signa;  % NOTE: sign(0) = 0 in MATLAB
    end 
    s = 0;
    r = absa;
    return

  elseif a == 0
    c = 0;
    s = signb;
    r = absb;
    return
  end

  %...........................
  % Both a and b are non-zero  
  %...........................
  if absb > absa
    t = a/b;
    s = signb / sqrt(1 + t^2); 
    c = s*t;
    r = b/s; % computationally better than d = a / c since |c| <= |s|

  else 
    t = b/a; 
    c = signa / sqrt(1 + t^2); 
    s = c*t;
    r = a/c; % computationally better than d = b / s since |s| <= |c|
  end

  return
end

%---------------------------------
% a and/or b are complex numbers
%---------------------------------
  %...........................
  % Special cases: a or b is 0
  %...........................
  if b == 0
    c = 1; 
    s = 0;
    r = a;
    return

  elseif a == 0
    c = 0;
    s = 1;
    r = b;
    return
  end

  %...........................
  % Both a and b are non-zero  
  %...........................
  if absb > absa
    t = absa/absb;
    c = 1/sqrt(1+t^2); % temporary
    s = c*conj(signb/signa); 
    c = c*t;
    r = b/conj(s);  
  else 
    t = absb/absa; 
    c = 1/sqrt(1+t^2); 
    s = c*t*conj(signb/signa);
    r = a/c; 
  end
end %BvdS:  added this.
end %BvdS: added this