function [x,flag,relres,iter,resvec] = mycgs(A,b,tol,tolmethod,maxit,maxtime,x0,callback,varargin)
%MYCGS   Conjugate Gradients Squared Method.
%   See CGS for the documentation.  This version behaves exactly the same
%   except that it eliminates some unnecessary code for options we never
%   use, adds a limit on the allowed wallclock time, and prints out an
%   indication of progress.

% mycgs assumed that before it was called, a message was written to the
% command window that was not terminated by a newline.  Hence neednewline
% is initially true.
neednewline = true;
havecallback = nargin >= 8;

if (nargin < 2)
   error('MATLAB:cgs:NotEnoughInputs', 'Not enough input arguments.');
end

% Check matrix and right hand side vector inputs have appropriate sizes
[m,n] = size(A);
if (m ~= n)
  error('MATLAB:mycgs:NonSquareMatrix', 'Matrix must be square.');
end
if ~isequal(size(b),[m,1])
  error('MATLAB:mycgs:RSHsizeMatchCoeffMatrix', ...
     ['Right hand side must be a column vector of' ...
     ' length %d to match the coefficient matrix.'],m);
end

% Assign default values to unspecified parameters
if (nargin < 3) || isempty(tol)
   tol = 1e-6;
end
if (nargin < 4) || isempty(maxit)
   maxit = min(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(mfilename(),tol,maxit,0,flag,iter,NaN);
   end
   if neednewline
       fprintf( 1, '\n' );
       neednewline = false;
   end
   fprintf( 1, '%s returning: right hand side is zero.\n', mfilename() );
   return
end

if ((nargin >= 7) && ~isempty(x0))
   if ~isequal(size(x0),[n,1])
      warning('MATLAB:cgs:WrongInitGuessSize', ...
         ['Initial guess has length %d, problem size is %d.'],prod(size(x0)),n);
      x = zeros(n,1);
   else
      x = x0(:);
   end
else
   x = zeros(n,1);
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
r = b - A*x;
normr = norm(r);                   % Norm of residual

NORMERROR = strcmp(tolmethod, 'norm');
if NORMERROR
    tolb = tol * n2b;                  % Relative tolerance
    endcondition = normr <= tolb;
else
    tolb = tol * max(abs(b));                  % Relative tolerance
    endcondition = max(abs(r)) <= tolb;
end

if endcondition                 % Initial guess is a good enough solution
   flag = 0;
   relres = normr / n2b;
   iter = 0;
   resvec = normr;
   if (nargout < 2)
      itermsg(mfilename(),tol,maxit,0,flag,iter,relres);
   end
   fprintf( 1, '%s returning: initial guess was close enough: normr %f, tolb %f.\n', ...
       mfilename(), normr, tolb );
   return
end

rt = r;                            % Shadow residual
resvec = zeros(maxit+1,1);         % Preallocate vector for norms of residuals
resvec(1) = normr;                 % resvec(1) = norm(b-A*x0)
normrmin = normr;                  % Norm of residual from xmin
rho = 1;
stag = 0;                          % stagnation of the method

% loop over maxit iterations (unless convergence or failure)

startTime = clock();
elapsedTime = 0;
ETIMEFLAG = 20;

itersPerDot = 20;
dotsPerLine = 50;
itersPerLine = itersPerDot*dotsPerLine;
for i = 1 : maxit
   rho1 = rho;
   rho = rt' * r;
   if (rho == 0) || isinf(rho)
      flag = 4;
      break
   end
   if i == 1
      u = r;
      p = u;
   else
      beta = rho / rho1;
      if (beta == 0) || isinf(beta)
         flag = 4;
         break
      end
      u = r + beta * q;
      p = u + beta * (q + beta * p);
   end
   ph1 = p;
   ph = ph1;
   vh = A*ph;
   rtvh = rt' * vh;
   if rtvh == 0
      flag = 4;
      break
   else
      alpha = rho / rtvh;
   end
   if isinf(alpha)
      flag = 4;
      break
   end
   if alpha == 0                    % stagnation of the method
      stag = 1;
   end
   q = u - alpha * vh;
   uh1 = u+q;
   uh = uh1;

   % Check for stagnation of the method
   if stag == 0
      stagtest = zeros(n,1);
      ind = (x ~= 0);
      stagtest(ind) = uh(ind) ./ x(ind);
      stagtest(~ind & uh ~= 0) = Inf;
      if abs(alpha)*norm(stagtest,inf) < eps
         stag = 1;
      end
   end

   x = x + alpha * uh;              % form the new iterate
   errvec = b - A*x;
   normr = norm(errvec);
   
   resvec(i+1) = normr;

    if NORMERROR
        endcondition = normr <= tolb;
    else
        endcondition = max(abs(errvec)) <= tolb;
    end

   if endcondition                 % 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
   qh = A*uh;
   r = r - alpha * qh;

   if (havecallback && (mod(i,10)==0) && callback(varargin{:}))
      if neednewline
          fprintf( 1, '\n' );
      end
      fprintf( 1, 'mycgs: detected stop on iteration %d\n', i );
      flag = 8;
      break;
   end
   
   if mod(i,itersPerLine)==1
       fprintf( 1, '%6.3g (%g/%g) ', normr / n2b, normr, n2b );
       neednewline = true;
       if maxtime > 0
           elapsedTime = etime( clock(), startTime );
           if elapsedTime > maxtime
               flag = ETIMEFLAG;
               break;
           end
       end
   end
   if mod(i,itersPerDot)==0
       fprintf( 1, '.' );
       neednewline = mod(i,itersPerLine) ~= 0;
       if ~neednewline % really
           fprintf( 1, '\n' );
       end
   end
end                                % for i = 1 : maxit
if neednewline
    fprintf( 1, ' %.3g\n', normr / n2b );
else
    fprintf( 1, '%6.3g\n', normr / n2b );
end

% 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(mfilename(),tol,maxit,i,flag,iter,relres);
end
