function x = solve_with_gepp(A,b)
% x = solve_with_gepp(A,b)
% Solve a linear system A*x=b for x using
% Gaussian elimination with partial pivoting

[n,n] = size(A);
L = zeros(n,n);
p = zeros(n,1);


% factor A = P*L*U.  Overwrite A with U.

for k = 1 : n
  [scrap,pivrowoffset] = max(abs(A(k:n,k)));
  pivrow = k + pivrowoffset - 1;
  A([pivrow,k], k:n) = A([k,pivrow], k:n);
  L([pivrow,k], 1:k-1) = L([k,pivrow], 1:k-1);
  p(k) = pivrow;
  piv = A(k,k);
  if piv == 0.0
    error('Singular matrix!')
  end
  L(k+1:n,k) = A(k+1:n,k) / piv;
  L(k,k) = 1;
  A(k+1:n,k:n) = A(k+1:n,k:n) - L(k+1:n,k) * A(k,k:n);
end

% Apply swaps

for k = 1 : n
  if k ~= p(k)
    b([k,p(k)]) = b([p(k),k]);
  end
end

% Forward sub with L

w = zeros(n,1);

for j = 1 : n
  w(j) = (b(j) - L(j,1:j-1) * w(1:j-1)) / L(j,j);
end

% Backsub with U.

x = zeros(n,1);

for j = n : -1 : 1
  x(j) = (w(j) - A(j,j+1:n) * x(j+1:n)) / A(j,j);
end



