function[L,U,P] = decomp_lu(A)

[n,m] = size(A);
L = zeros(n);
P = eye(n);

for r = 1:n-1

    [val,pos] = max(abs(A(r:n,r)));
    linha = pos + (r - 1);

    if linha ~= r
    aux = A(r,:);
    A(r,:) = A(linha,:);
    A(linha,:) = aux;

    aux = P(r,:);
    P(r,:) = P(linha,:);
    P(linha,:) = aux;

    aux = L(r,:);
    L(r,:) = L(linha,:);
    L(linha,:) = aux;
    end

  L(r,r) = 1;
  for i = r + 1 : n
    mir = - A(i,r)/A(r,r);
    L(i,r) = -mir;
    for j = r:m
      A(i,j) = A(i,j) + mir * A(r,j);
    end
  end

end
L(n,n) = 1;
U = A;

end