function  [z0,sigma_sequence,iterations] = gmres_x(A,b,L,U,P,tol,varargin)
% Input:  A,b, L,U,P, tol, [z0 = initial estimate of soln, which is optional]
% Output: z0,sigma_sequence,iterations
% Function gmres_x.m is a version of GMRES organized as suggested by
% R.J. Hanson and D.R. Kincaid,"Notes on GMRES Algorithm Organization",
% Technical Report TR-05-05, March 2005,
% Computer Sciences Department, University of Texas at Austin.
% Date = 18 April 2005.

% Prepare to enter busy loop that computes GMRES(itmax) approximation.
% Inline form of GMRES with left and right pre-conditioners.

while 1
    itmax = input(' Enter number ( > 0) of inner GMRES iterations (itmax) = ');
    if(itmax > 0)
        break
    end
end
[n,n]=size(A);
V = zeros(n,itmax); sigma = zeros(itmax+1,1);
H = zeros(itmax*(itmax+1)/2,1); %B = zeros(itmax,itmax);
n_in=nargin;
if(n_in == 6)
    % No initial soln estimate is provided.  Start with randomness.
    z0 = rand(n,1);
elseif(n_in == 7)
    % Use provided initial soln estimate.
    z0=varargin{1};
    [nn,kk]=size(z0);
    if(nn ~= n)
        disp('Error - Initial estimate must have size that matches size of A')
        return
    end
end
iterations = 0;
% Overall GMRES(itmax) iteration loop ---------------------------------------
while 1
  x = z0;            % Start the iteration at an approximate solution.
  sigma = 0*sigma;
  H = 0*H;
  % Inner GMRES iteration phase starts here.---------------------------------
  for j=0:itmax
    x = A*x;   % One place that operator gets applied
    if(j==0)   % An organizational difference from the
      x = b-x;   % other published GMRES algorithms is here.
    end
    x = L\(P*x);   % One place that left pre-conditioner (M_L^-1) gets applied.
    % MGS loop is done here.  Note that the loop is skipped when j==0.
    for i=1:j
      t = x'*V(:,i);
      % Upper triangular matrix is stored in H, with row major order.
      % During the solution phase this allows for unit strides in applying the
      % plane rotations to the matrix rows.
      x = x-t*V(:,i); H(itmax*(i-1)-(i*(i-1))/2+j) = t;
    end
    sigma(j+1) = norm(x);
    % Stop iteration if residual is small enough or have come to maximum number
    % of interations.  If the initial value z_0 is a solution, this will exit
    % immediately with j==0.
    if(sigma(j+1)<=tol || j==itmax)
      break
    else
      x = x*(1/sigma(j+1));
      V(:,j+1) = x;
    end
    x = U\x;      % One place that right pre-conditioner (M_R^_1) gets applied.
  end
  % Inner GMRES iteration phase ends here.-----------------------------------
  % Solution phase starts here. ---------------------------------------------
  if(j >= 0)
    y = zeros(j+1,1); y(1) = sigma(1);
  end
  % Solve least squares system with upper Hessenberg matrix.  Note that the
  % sub-diagonal terms of H are stored in the array sigma(:), offset by 1.
  k = 1;
  for i=1:j
    % Note the unit stride when applying the plane rotation.
    m = k+itmax-i+1; in = j-i;
    % The lower diagonal terms are kept in array sigma(2:j+1).
    [H(k),c,s] = rotg(H(k),sigma(i+1));
    T =            c*H(k+1:k+in)+s*H(m:m+in-1);
    H(m:m+in-1) = -s*H(k+1:k+in)+c*H(m:m+in-1);
    H(k+1:k+in) = T;
    % Apply rotation to the right hand side.  Use fact that y(i+1) = 0.
    t  =  c*y(i);
    y(i+1) = -s*y(i);
    y(i) = t;
    % This test leaves k at the right value for the back solve step.
    if(i==j)
      break
    end
    k = m;
  end

  if(j > 0)
    % System is now upper triangular.
    % Back solve using dot products since these involve a unit stride.
    % The current value of k is correct because of the triangularization step.
    y(j) = y(j)/H(k);
    for i=j-1:-1:1
      k = k-(itmax-i+1);
      y(i) = (y(i)-H(k+1:k+j-i)'*y(i+1:j))/H(k);
    end
    x = V(:,1:j)*y(1:j);
    x = U\x;       % Account for right pre-conditioner coordinate change here.
    z0 = z0+x;     % Update solution.
  end
  % Solution phase ends here.------------------------------------------------
  iterations = iterations + 1;

  sigma_sequence(iterations) = sigma(1);
  % Decreasing values of sigma(1) imply convergence.  Note if
  % stagnation has onset.
  if(sigma(1) <= tol || abs(sigma(1)-abs(y(j+1))) <= sqrt(eps)*sigma_sequence(1))
    sigma_sequence(iterations+1) = abs(y(j+1));
    break
  end
end