'Compute a fraction-free LU decomposition.
Returns 4 matrices P, L, D, U such that PA = L D**-1 U.
If the elements of the matrix belong to some integral domain I, then all
elements of L, D and U are guaranteed to belong to I.
**Reference**
- W. Zhou & D.J. Jeffrey, "Fraction-free matrix factors: new forms
for LU and QR factors". Frontiers in Computer Science in China,
Vol 2, no. 1, pp. 67-80, 2008.
See Also
LUdecomposition
LUdecomposition_Simple
LUsolve'

def LUdecompositionFF(self):

    from sympy.matrices import SparseMatrix
    zeros = SparseMatrix.zeros
    eye = SparseMatrix.eye
    (n, m) = (self.rows, self.cols)
    (U, L, P) = (self.as_mutable(), eye(n), eye(n))
    DD = zeros(n, n)
    oldpivot = 1
    for k in range((n - 1)):
        if (U[(k, k)] == 0):
            for kpivot in range((k + 1), n):
                if U[(kpivot, k)]:
                    break
            else:
                raise ValueError('Matrix    is    not    full    rank')
            (U[k, k:], U[kpivot, k:]) = (U[kpivot, k:], U[k, k:])
            (L[k, :k], L[kpivot, :k]) = (L[kpivot, :k], L[k, :k])
            (P[k, :], P[kpivot, :]) = (P[kpivot, :], P[k, :])
        L[(k, k)] = Ukk = U[(k, k)]
        DD[(k, k)] = (oldpivot * Ukk)
        for i in range((k + 1), n):
            L[(i, k)] = Uik = U[(i, k)]
            for j in range((k + 1), m):
                U[(i, j)] = (((Ukk * U[(i, j)]) - (U[(k, j)] * Uik)) / oldpivot)
            U[(i, k)] = 0
        oldpivot = Ukk
    DD[((n - 1), (n - 1))] = oldpivot
    return (P, L, DD, U)
