function [x,P,L,U] = LU_rowPivot(A,b,n)
% LU factorization of an n by n matrix A
% using Gauss elimination with row pivoting.
% b is an n by 1 vector.
% x is the solution to the linear system Ax = b.
% P is the permutation matrix.
% L is the unit lower triangular matrix.
% U is upper triangular and is stored in the upper triangular
% part of A; the remaining part has no special meaning
% and must be zeroed out to get U (not done here since
% back-substitution does not read those elements any way).
% Note the matrices must be such that: PA = LU.
% K. Ming Leung, 02/23/04

L = eye(n); P = eye(n); 
A = [A b];  % form augmented system 
for k = 1:n-1
    [pivot, p] = max(abs(A(k:n,k)));
    p = k+p-1;
    if (p ~= k)
        A = swapRow(A,k,p);
        L(:,1:k-1) = swapRow(L(:,1:k-1),k,p);
        P = swapRow(P,k,p);
    end
    
    if (A(k,k) ~= 0) 
        L(k+1:n,k) = A(k+1:n,k)/A(k,k);
        i = k+1:n; j = k+1:n+1;
        A(i,j) = A(i,j) - L(i,k)*A(k,j);
    end
end

% 
b = A(:,n+1); U = triu(A(:,1:n)); 
x = backSubstitution(U,b,n);


     