function [P,L,U]=lu_fac_pivot(A)
% PA = LU factorization for matrix A
% L is lower triangular
% U is upper triangular
% P is permutation matrix that permutes A into stable form
  
  n=length(A);U=A; P=eye(n);
  for k=1:n-1 % Loop over columns
    [y,kk]=max(abs(A(k:n,k))) ; kk = kk + k-1;% Find coordinate kk of max element
    if kk ~= k % Interchage rows kk and k of matrix and permutation matrix
      temp = A(k,:) ; A(k,:) = A(kk,:) ; A(kk,:) = temp;
      temp = P(k,:) ; P(k,:) = P(kk,:) ; P(kk,:)  = temp;
    end
    if A(k,k)==0 
      continue % Skip current column if it's already zero
    end
    A(k+1:n,k)=A(k+1:n,k)/A(k,k); % Compute multipliers for column k
    for j = k+1:n % Apply elimination matrix
	A(k+1:n,j) = A(k+1:n,j) - A(k+1:n,k)*A(k,j);		
    end
  end
  U=triu(A);           % Upper trianguler matrix
  L=tril(A,-1)+eye(n); % Lower triangular matrix