clear ;

rand('seed', 1234); % Initialize uniform random number generator

n = 8; % Matrix dimension

A = rand(n,n); % Create random matrix
xsol = rand(n,1);
b = A*xsol;

P = eye(n);
for k=1:n-1
  % Choose best pivot
  [piv, m] = max(abs(A(k:n,k)));
  m = m+k-1;

  % Apply permutation of rows in matrix
  A([k,m],:) = A([m,k],:);

  % Save permutation 
  P([k,m],:) = P([m,k],:);

  % Eliminate the column
  for i=k+1:n
    A(i,k) = A(i,k)/A(k,k);
    A(i,k+1:n) = A(i,k+1:n) - A(i,k)*A(k,k+1:n);
  end
end

% Verifications

U = triu(A); % Extract the upper triangular factor
L = tril(A); % Extract lower triangular part
D = diag(diag(A)); % Create diagonal matrix
L = L - D + eye(n); % Construct the actual lower triangular factor

%Solving the system
y = L \ (P*b);
x = U \ y;
norm(x - xsol)/norm(xsol)