function [L,U, P, A] = MyLUpp(A)
%
%  LU factorization with partial pivoting
%  call as [L, U , P, B] = MyLUpp(A)
%  then    L*U = P*A
%  9/20/2006
%

[m,n] = size(A);
msg = 0 ; % for warning message
Amax = max(abs(A(:)));
tol = sqrt(eps)*Amax; % tolerance for small pivots
p = 1:m;  % tracks permutations

for j = 1:m-1
    [piv, k] = max( abs( A(j:m,j) ) ) ;   
    if piv == 0
        fprintf('   zero pivot encountered in step %8.0f \n',j)
    else
        if abs(piv) <= tol, msg = 1; end;
        if k ~= 1  % row changes
            temp = A(j,:);
            A(j,:) = A(j+k-1,:);
            A(j+k-1,:) = temp;
            temp = p(j); p(j) = p(j+k-1); p(j+k-1)=temp;
        end
        
        for k = j+1:m
            alpha = A(k,j)/A(j,j);
            A(k,j) = alpha;
            for i = j+1:n
                A(k,i) = A(k,i)-alpha*A(j,i);
            end
        end
    end
end

if msg, fprintf('  Warning: Small pivots were encountered  \n'), end

%  Build matrices for the output

U = zeros(size(A));
for i=1:m
    for j = i:n
        U(i,j)=A(i,j);
    end
end

L = eye(m);
for i=2:m
    for j= 1:(i-1)
        L(i,j) = A(i,j);
    end
end

P = zeros(m,m);
for j = 1:m
    P(j,p(j))=1;
end