def lup(matrix):
    assert matrix.shape[0] == matrix.shape[1]
    n = matrix.shape[0]
    L = np.eye(n, dtype=matrix.dtype)
    U = np.copy(matrix)
    P = np.eye(n, dtype=matrix.dtype)
    for i in range(n):
        if np.abs(U[i,i]) < tol:
            for k in range(i + 1, n):
                if np.abs(U[k, i]) > tol:
                    U[[i,k]] = U[[k, i]]
                    L[[i,k]] = L[[k, i]]
                    L[:, [i,k]] = L[:, [k,i]]
                    P[[i,k]] = P[[k, i]]
                    break

        for j in range(i + 1, n):
            coef = U[j,i]/U[i,i]
            U[j, :] -= coef * U[i, :]
            L[j, i] = coef
    return L, U, P