function [ x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm ] = ...
           minres( n, b, Aname, Mname, iw, rw,  ...
                   precon, shift, show, check, itnlim, rtol )

%        [ x, istop, itn, rnorm, Arnorm, Anorm, Acond, ynorm ] = ...
%          minres( n, b, A, M, iw, rw,  ...
%                  precon, shift, show, check, itnlim, rtol )
%
%  minres solves the system of linear equations Ax = b
%  or the least squares problem           min ||Ax - b||_2^2,
%  where A is a symmetric matrix (possibly indefinite or singular)
%  and b is a given vector.
%  
%  "A" may be a dense or sparse matrix (preferably sparse!)
%  or the name of a function such that
%               y = A( n, x, iw, rw )
%  returns the product y = Ax for any given vector x.
%
%  If precon is true,
%  "M" defines a positive-definite preconditioner M = C C'.
%  "M" may be a dense or sparse matrix (preferably sparse!)
%  or the name of a function such that
%               y = M( n, x, iw, rw )
%  solves the system My = x for any given vector x.
%
%  WARNING:   The files containing the functions "A" and "M"
%             must not be called Aname.m or Mname.m !!!!
%
%  If shift ~= 0, minres really solves (A - shift*I)x = b
%  or the corresponding least squares problem if shift is an
%  eigenvalue of A.
%

%  Code authors:Michael Saunders, SOL, Stanford University
%               Sou Cheng Choi,  SCCM, Stanford University
%
%  02 Sep 2003: Date of Fortran 77 version, based on 
%               C. C. Paige and M. A. Saunders (1975),
%               Solution of sparse indefinite systems of linear equations,
%               SIAM J. Numer. Anal. 12(4), pp. 617-629.
%
%  02 Sep 2003: ||Ar|| now estimated as Arnorm.
%  17 Oct 2003: f77 version converted to MATLAB.
%  ------------------------------------------------------------------

%  Initialize                               

first = 'Enter minres.   ';
last  = 'Exit  minres.   ';
space = ' ';
true  = 1;
false = 0;
msg   =[' beta2 = 0.  If M = I, b and x are eigenvectors    '   % -1
        ' beta1 = 0.  The exact solution is  x = 0          '   %  0
        ' A solution to Ax = b was found, given rtol        '   %  1
        ' A least-squares solution was found, given rtol    '   %  2
        ' Reasonable accuracy achieved, given eps           '   %  3
        ' x has converged to an eigenvector                 '   %  4
        ' acond has exceeded 0.1/eps                        '   %  5
        ' The iteration limit was reached                   '   %  6
        ' Aname  does not define a symmetric matrix         '   %  7
        ' Mname  does not define a symmetric matrix         '   %  8
        ' Mname  does not define a pos-def preconditioner   ']; %  9
 
if show
   disp( space )
   disp( [ first 'Solution of symmetric Ax = b' ] )
   disp( sprintf( 'n      =  %3g     precon =  %4g           shift  =  %23.14e',...
                   n,                precon,                 shift ) )
   disp( sprintf( 'itnlim =  %3g     rtol   =  %11.2e\n' ,...
                   itnlim,           rtol ) )
end

istop = 0;   itn   = 0;   Anorm = 0;    Acond = 0;
rnorm = 0;   ynorm = 0;   done  = false;
x     = zeros(n,1);

%---------------------------------------------------------------------
% Decode Aname.
%---------------------------------------------------------------------
operator  =  ischar(Aname) | isa(Aname,'function_handle');
explicitA = ~operator;

if explicitA            % assume Aname is an explicit matrix A.
  nnzA   = nnz(Aname);
  if  issparse(Aname)
    fprintf('The matrix A is an explicit sparse matrix')
    fprintf('\nnnz(A) =%9g', nnzA)
  else
    fprintf('The matrix A is an explicit dense matrix' )
  end
else
  if ischar(Aname)
    fname = Aname;
  else  % assume Aname is a function handle
    fname = func2str(@Aname);	 
  end
  disp(['The matrix A is an operator defined by ' fname])
end

%---------------------------------------------------------------------
% Decode Mname.
%---------------------------------------------------------------------
if precon
  operator  =  ischar(Mname) | isa(Mname,'function_handle');
  explicitM = ~operator;

  if explicitM            % assume Mname is an explicit matrix M.
    nnzM   = nnz(Mname);
    if  issparse(Mname)
      fprintf('\nThe matrix M is an explicit sparse matrix')
      fprintf('\nnnz(M) =%9g', nnzM)
    else
      fprintf('\nThe matrix M is an explicit dense matrix' )
    end
  else
    if ischar(Mname)
      fname = Mname;
    else  % assume Mname is a function handle
      fname = func2str(@Mname);	 
    end
    disp(['The matrix M is an operator defined by ' fname])
  end
end % if precon

%------------------------------------------------------------------
% Set up y and v for the first Lanczos vector v1.
% y  =  beta1 P' v1,  where  P = C**(-1).
% v is really P' v1.
%------------------------------------------------------------------
y     = b;
r1    = b;
if precon, y = minresxxxM( Mname, n, b, iw, rw ); end
beta1 = b'*y;
                                                
%  Test for an indefinite preconditioner.
%  If b = 0 exactly, stop with x = 0.

if beta1< 0, istop = 8;  show = true;  done = true; end
if beta1==0,             show = true;  done = true; end

if beta1> 0
  beta1  = sqrt( beta1 );       % Normalize y to get v1 later.

  % See if M is symmetric.

  if check & precon
    r2     = minresxxxM( Mname, n, y, iw, rw );
    s      = y' *y;
    t      = r1'*r2;
    z      = abs( s - t );
    epsa   = (s + eps) * eps^(1/3);
    if z > epsa, istop = 7;  show = true;  done = true; end
  end

  % See if A is symmetric.

  if check
     w    = minresxxxA( Aname, n, y, iw, rw );
     r2   = minresxxxA( Aname, n, w, iw, rw );
     s    = w' *w;
     t    = y' *r2;
     z    = abs( s - t );
     epsa = (s + eps) * eps^(1/3);
     if z > epsa, istop = 6;  done  = true;  show = true;  end
  end
end

%------------------------------------------------------------------
% Initialize other quantities.
% ------------------------------------------------------------------
oldb   = 0;       beta   = beta1;   dbar   = 0;       epsln  = 0;
qrnorm = beta1;   phibar = beta1;   rhs1   = beta1;
rhs2   = 0;       tnorm2 = 0;       ynorm2 = 0;
cs     = -1;      sn     = 0;
w      = zeros(n,1);
w2     = zeros(n,1);
r2     = r1;

if show
  disp(' ')
  disp(' ')
  head1 = '   Itn     x(1)     Compatible    LS';
  head2 = '       norm(A)  cond(A)';
  head2 = [head2 ' gbar/|A|'];   %%%%%% Check gbar
  disp( [head1 head2] )
end

%---------------------------------------------------------------------
% Main iteration loop.
% --------------------------------------------------------------------
if ~done                              % k = itn = 1 first time through
  while itn < itnlim
    itn    = itn  +  1;

    %-----------------------------------------------------------------
    % Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
    % The general iteration is similar to the case k = 1 with v0 = 0:
    %
    %   p1      = Operator * v1  -  beta1 * v0,
    %   alpha1  = v1'p1,
    %   q2      = p2  -  alpha1 * v1,
    %   beta2^2 = q2'q2,
    %   v2      = (1/beta2) q2.
    %
    % Again, y = betak P vk,  where  P = C**(-1).
    % .... more description needed.
    %-----------------------------------------------------------------
    s = 1/beta;                 % Normalize previous vector (in y).
    v = s*y;                    % v = vk if P = I

    y      = minresxxxA( Aname, n, v, iw, rw );
    y      = (- shift)*v + y;
    if itn >= 2
      y = y - (beta/oldb)*r1;
    end

    alfa   = v'*y;              % alphak
    y      = (- alfa/beta)*r2 + y;
    r1     = r2;
    r2     = y;
    if precon,  y = minresxxxM( Mname, n, r2, iw, rw );  end
    oldb   = beta;              % oldb = betak
    beta   = r2'*y;             % beta = betak+1^2
    if beta < 0, istop = 6;  break;  end
    beta   = sqrt( beta );
    tnorm2 = tnorm2 + alfa^2 + oldb^2 + beta^2;

    if itn==1                   % Initialize a few things.
      if beta/beta1 <= 10*eps   % beta2 = 0 or ~ 0.
	istop = -1;             % Terminate later.
      end
      %tnorm2 = alfa**2  ??
      gmax   = abs( alfa );     % alpha1
      gmin   = gmax;            % alpha1
    end

    % Apply previous rotation Qk-1 to get
    %   [deltak epslnk+1] = [cs  sn][dbark    0   ]
    %   [gbar k dbar k+1]   [sn -cs][alfak betak+1].

    oldeps = epsln;
    delta  = cs * dbar  +  sn * alfa; % delta1 = 0         deltak
    gbar   = sn * dbar  -  cs * alfa; % gbar 1 = alfa1     gbar k
    epsln  =               sn * beta; % epsln2 = 0         epslnk+1
    dbar   =            -  cs * beta; % dbar 2 = beta2     dbar k+1
    root   = norm([gbar dbar]);
    Arnorm = phibar * root;

    % Compute the next plane rotation Qk

    gamma  = norm([gbar beta]);       % gammak
    gamma  = max([gamma eps]);
    cs     = gbar / gamma;            % ck
    sn     = beta / gamma;            % sk
    phi    = cs * phibar ;            % phik
    phibar = sn * phibar ;            % phibark+1

    % Update  x.

    denom = 1/gamma;
    w1    = w2;
    w2    = w;
    w     = (v - oldeps*w1 - delta*w2) * denom;
    x     = x  +  phi*w;

    % Go round again.

    gmax   = max([gmax gamma]);
    gmin   = min([gmin gamma]);
    z      = rhs1 / gamma;
    ynorm2 = z^2  +  ynorm2;
    rhs1   = rhs2 -  delta*z;
    rhs2   =      -  epsln*z;

    % Estimate various norms and test for convergence.

    Anorm  = sqrt( tnorm2 );
    ynorm  = sqrt( ynorm2 );
    epsa   = Anorm * eps;
    epsx   = Anorm * ynorm * eps;
    epsr   = Anorm * ynorm * rtol;
    diag   = gbar;
    if diag==0, diag = epsa; end

    qrnorm = phibar;
    rnorm  = qrnorm;
    test1  = rnorm / (Anorm*ynorm);    %  ||r|| / (||A|| ||x||)
    test2  = root  /  Anorm;           % ||Ar|| / (||A|| ||r||)

    % Estimate  cond(A).
    % In this version we look at the diagonals of  R  in the
    % factorization of the lower Hessenberg matrix,  Q * H = R,
    % where H is the tridiagonal matrix from Lanczos with one
    % extra row, beta(k+1) e_k^T.

    Acond  = gmax/gmin;

    % See if any of the stopping criteria are satisfied.
    % In rare cases, istop is already -1 from above (Abar = const*I).

    if istop==0
      t1 = 1 + test1;      % These tests work if rtol < eps
      t2 = 1 + test2;
      if t2    <= 1      , istop = 2; end
      if t1    <= 1      , istop = 1; end
      
      if itn   >= itnlim , istop = 6; end
      if Acond >= 0.1/eps, istop = 4; end
      if epsx  >= beta1  , istop = 3; end
     %if rnorm <= epsx   , istop = 2; end
     %if rnorm <= epsr   , istop = 1; end
      if test2 <= rtol   , istop = 2; end
      if test1 <= rtol   , istop = 1; end
    end


    % See if it is time to print something.

    prnt   = false;
    if n      <= 40       , prnt = true; end
    if itn    <= 10       , prnt = true; end
    if itn    >= itnlim-10, prnt = true; end
    if mod(itn,10)==0     , prnt = true; end
    if qrnorm <= 10*epsx  , prnt = true; end
    if qrnorm <= 10*epsr  , prnt = true; end
    if Acond  <= 1e-2/eps , prnt = true; end
    if istop  ~=  0       , prnt = true; end

    if show & prnt
      str1 = sprintf ( '%6g %12.5e %10.3e', itn, x(1), test1 );
      str2 = sprintf ( ' %10.3e',           test2 );
      str3 = sprintf ( ' %8.1e %8.1e',      Anorm, Acond );
      str3 = [str3 sprintf( ' %8.1e', gbar/Anorm)];
      disp( [str1 str2 str3] )

      if istop > 0, break; end

      if mod(itn,10)==0, disp(' '); end

%     if (debug) then  % Print true Arnorm
%        call Aprod ( n, x, v )               % v = A*x
%        do i = 1, n
%           v(i) = b(i) - v(i) + shift*x(i)   % v = b - (A - shift*I)*x
%        end do
%        call Aprod ( n, v, ww )             %% ww = A*v = Ar
%        t2true = dnrm2 ( n, ww, 1 ) / (Anorm * rnorm)
%        write(nout, '(25x, a, 1p, e10.2)') 'True test2 =',t2true
%      end if

    end % show & prnt
  end % main loop
end % if ~done early

% Display final status.

if show
  disp( space )
  disp( [ last sprintf( ' istop   =  %3g               itn   =%5g',...
			istop,                      itn  ) ] )
  disp( [ last sprintf( ' Anorm   =  %12.4e      Acond =  %12.4e' ,...
			  Anorm,                 Acond  ) ] )
  disp( [ last sprintf( ' rnorm   =  %12.4e      ynorm =  %12.4e' ,...
			  rnorm,                 ynorm  ) ] )
  disp( [ last sprintf( ' Arnorm  =  %12.4e', Arnorm ) ] )
  disp( [ last msg(istop+2,:) ])
end
%-----------------------------------------------------------------------
% End function pdco.m
%-----------------------------------------------------------------------


function y = minresxxxA( Aname, n, x, iw, rw )

%        y = minresxxxA( Aname, n, x, iw, rw )
%    computes y = Ax for a matrix A defined by input parameter Aname.

%-----------------------------------------------------------------------
% 17 Oct 2003: Default A*x function for minres.m.
%              User parameters Aname, n, iw, rw are passed thru to here.
%-----------------------------------------------------------------------

  if ischar(Aname)
    y = feval( Aname, mode, n, x, iw, rw );
  else
    y = Aname*x;
  end

%-----------------------------------------------------------------------
% End private function minresxxxA
%-----------------------------------------------------------------------


function y = minresxxxM( Mname, n, x, iw, rw )

%        y = minresxxxM( Mname, n, x, iw, rw )
%    solves My = x for a matrix M defined by input parameter Mname.

%-----------------------------------------------------------------------
% 17 Oct 2003: Default M\x for minres.m.
%              User parameters Mname, n, iw, rw are passed thru to here.
%-----------------------------------------------------------------------

  if ischar(Mname)
    y = feval( Mname, n, x, iw, rw );
  else
    y = Mname\x;
  end

%-----------------------------------------------------------------------
% End private function minresxxxM
%-----------------------------------------------------------------------