% solve Ax=b for x using Gaussian elimination with partial pivoting (GEPP)
% function [x,PG] = geppsolve(A,b)
% Inputs:
%   A - matrix
%   b - vector
% Outputs:
%   x - solution of A*x=b, computed usign GENP
%   PG - pivot growth, defined as the
%        largest entry in U divided by the largest entry in A
%        (note that the largest entry in L is 1, because of pivoting)
%        Should be at least 1. Expect an inaccurate answer if >> 1.
function [x,PG] = geppsolve(A,b)
%
% Use Gaussian elimination to compute L, U and P such that P*A=L*U and
%   P is a permutation, i.e. a matrix whose rows are the rows of
%     the identity in a permuted order. 
%     This means P*A is the same as A with its rows in permuted order.
%   L is lower triangular with ones on its diagonal
%   U is upper triangular 
n=max(size(A));  % dimension of A
[L,U,P] = lu(A);
%
% Compute pivot growth
PG = max(max(abs(L)))*max(max(abs(U)))/max(max(abs(A)));
%
% Solve L*y=P*b for y using forward substitution
bb = P*b;
y(1) = bb(1);
for i=2:n,
   y(i) = (bb(i) - L(i,1:i-1)*y(1:i-1)')/L(i,i);
end
%
% Solve U*x=y
x(n) = y(n)/U(n,n);
for i=n-1:-1:1,
   x(i) = (y(i) - U(i,i+1:n)*x(i+1:n)')/U(i,i);
end
%
% Make x a column vector
x=x';