function [x, aa, P,f] = gepp(a,b)
% ge  Solve linear system a*x=b using Gaussian Elimination with partial pivoting
%
% Output variables:
% x = a^-1*b, the solution
% aa is an nxn matrix, the upper (including the diagonal) is the
% matrix U of the P*A=L*U decomposition of A, the lower (not including
% the diagonal, is the L matrix without the diagonal)
% P is the permutation matrix
% f=1, if A is invertible, f=0 if not


% To recover L and U
% l=eye(n);l=l+tril(aa,-1);
% u=triu(aa,0);



%function [pr,l,u,err]=gepp(a)

f=1;
n=max(size(a));
P=eye(n);
aa=a;
for i=1:n-1,
%  find largest absolute entry in 
%  the first column of the   
%  (i:n,i:n) submatrix
   [aux,I]=max(abs(aa(i:n,i)));
   if aux==0
       error('non invertible matrix')
   end
   imax=I(1)+i-1;
%  swap rows
   if (imax ~= i),
      temp = P(i,:);
      P(i,:) = P(imax,:);
      P(imax,:) = temp;
      temp = aa(i,:);
      aa(i,:) = aa(imax,:);
      aa(imax,:) = temp;
   end
%  eliminate
   aa(i+1:n,i) = aa(i+1:n,i)/aa(i,i);
   aa(i+1:n,i+1:n) = aa(i+1:n,i+1:n) - aa(i+1:n,i)*aa(i,i+1:n);
end
l=eye(n);l=l+tril(aa,-1);
u=triu(aa,0);




y = forbacksubs(l,P*b,0);
x = forbacksubs(u,y,1);
