function [L,U,p] = mylu2(A)
% LU factorization with partial pivoting

n = length(A);
p = 1:n;
for k=1:n-1
    % find max abs value
    amax = abs(A(k,k));
    imax = k;
    for i=k+1:n
        t = abs(A(i,k));
        if t > amax
            amax = t;
            imax = i;
        end
    end
    if amax == 0, error('zero division'), end
    % swap k-th row for imax-th row.
    if k ~= imax
        kt = p(k);
        p(k) = p(imax);
        p(imax) = kt;
        for j=1:n
            t = A(k,j);
            A(k,j) = A(imax,j);
            A(imax,j) = t;
        end
    end
    % part of L
    for i=k+1:n
        A(i,k) = A(i,k)/A(k,k);
    end
    % part of U
    for j=k+1:n
        for i=k+1:n
            A(i,j) = A(i,j) - A(i,k)*A(k,j);
        end
    end
end
L = eye(n) + tril(A,-1);
U = triu(A);