% solve Ax=b for x using Gaussian elimination with no pivoting (GENP)
% function [x,PG] = genpsolve(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 L times the
%        largest entry in U divided by the largest entry in A
%        Will be at least 1. Expect an inaccurate answer if >> 1.
function [x,PG] = genpsolve(A,b)
%
% Use Gaussian elimination to compute L and U such that A=L*U and
%   L is lower triangular with ones on its diagonal
%   U is upper triangular 
n=max(size(A));  % dimension of A
normA = max(max(abs(A)));
for i=1:n-1,
   A(i+1:n,i) = A(i+1:n,i)/A(i,i);
   A(i+1:n,i+1:n) = A(i+1:n,i+1:n) - A(i+1:n,i)*A(i,i+1:n);
end
L = eye(n)+tril(A,-1);
U = triu(A,0);
%
% Compute pivot growth
PG = max(max(abs(L)))*max(max(abs(U)))/normA;
%
% Solve L*y=b for y using forward substitution
y(1) = b(1);
for i=2:n,
   y(i) = (b(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';