matlup(M) =
{
  my (L = matid(#M), U = M, P = L);

  for (i = 1, #M-1,     \\ pivoting
    p = M[z=i,i];
    for (k = i, #M, if (M[k,i] > p, p = M[z=k,i]));

    if (i != z,         \\ swap rows
      k = U[i,]; U[i,] = U[z,]; U[z,] = k;
      k = P[i,]; P[i,] = P[z,]; P[z,] = k;
    );
  );

  for (i = 1, #M-1,     \\ decompose
    for (k = i+1, #M,
      L[k,i] = U[k,i] / U[i,i];
      for (j = i, #M, U[k,j] -= L[k,i] * U[i,j])
    )
  );

  [L,U,P]       \\ return L,U,P triple matrix
}