 function [L, U, p] = lufactP(A)
%----------------------------------------------
% function [L, U, p] = lufactP (A)
% LU factorization with partial pivoting
% on return A(p,:) = L*U
%---------------------------------------------- 
 n = size(A,1) ;
%-------------------- initialize permutation
 p = [1:n];    
 for k=1:n-1 
    fprintf('*** Step k =  %d  ***\n',k)
    p
    A
    pause 
%%
    [t, ip] = max(abs(A(k:n,k))); 
    ip = ip+k-1 ;
    fprintf('swap rows k=%d and i=%d\n',k,ip);
%%---------- swap rows 
    temp = A(k,:) ;
    A(k,:) =  A(ip,:); 
    A(ip,:) = temp;
%%----------update permutation (swap .. )
    i = p(k); p(k) = p(ip); p(ip) = i;
    p
    A
    pause 
%%---------- elimination
    for i=k+1:n
      piv = A(i,k) / A(k,k) ;
      A(i,k) = piv ;
      A(i,k+1:n) = A(i,k+1:n) - piv*A(k,k+1:n);
    end
 end 
%%
 p
 A
 pause 
%%----------- get L and U matrices
 L = eye(n) + tril(A,-1);
 U = triu(A);

