% solve Ax=b for x using Gaussian elimination with complete pivoting (GECP)
% function [x,PG] = gecpsolve(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] = gecpsolve(A,b)
%
% Use Gaussian elimination to compute L, U, Pr and Pc such that Pr'*A*Pc=L*U 
%   or equivalently A = Pr*L*U*Pc'.
%   Pr and Pc are permutations, i.e. a matrix whose rows are the rows of
%     the identity in a permuted order. 
%     This means Pr'*A*Pc is the same as A with its rows and columns 
%     in permuted order.
%   L is lower triangular with ones on its diagonal
%   U is upper triangular 
[Pr,Pc,L,U] = gecp(A);
%
% Compute pivot growth
PG = max(max(abs(U)))/max(max(abs(A)));
%
% Solve L*y=Pr*b for y using forward substitution
bb = Pr'*b;
y(1) = bb(1);
n=max(size(A));  % dimension of A
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';    
%
% Permute x
x = Pc*x;