parameters: (self,*vars)
return_statement: return P, L, DD, U
docstring: 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.
docstring_summary: Compute a fraction-free LU decomposition.
docstring_tokens: Compute a fraction - free LU decomposition .
function: def LUdecompositionFF(self,*vars):
          """
          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.
          """
          n, m = self.rows, self.cols
          U, L, P = self[:,:], eye(n), eye(n)
          DD = zeros(n) # store it smarter since it's just diagonal
          oldpivot = S.One
  
          for k in range(n-1):
              log.info('row=%d', k)
              if U[k,k] == 0:
                  for kpivot in range(k+1, n):
                      if U[kpivot, k] != 0:
                          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] = simplify((Ukk * U[i,j] - U[k,j]*Uik) / oldpivot)
                      D = Ukk * U[i,j] - U[k,j]*Uik
                      if len(vars) > 0 and D != S.Zero and not oldpivot.is_number:
                          try:
                              D,r = div(Poly(D,*vars),oldpivot)
                          except UnificationFailed:
                              log.warn('unification failed, trying direct division')
                              D /= oldpivot
                      else:
                          D /= oldpivot
                      # save
                      if D.is_Atom:
                          U[i,j] = D.as_expr()
                      else:
                          if len(vars) > 0:
                              U[i,j] = D.as_expr()
                          else:
                              U[i,j] = D.cancel()
                  U[i,k] = 0
              oldpivot = Ukk
          DD[n-1,n-1] = oldpivot
          return P, L, DD, U
function_tokens: def LUdecompositionFF ( self , * vars ) : n , m = self . rows , self . cols U , L , P = self [ : , : ] , eye ( n ) , eye ( n ) DD = zeros ( n ) # store it smarter since it's just diagonal oldpivot = S . One for k in range ( n - 1 ) : log . info ( 'row=%d' , k ) if U [ k , k ] == 0 : for kpivot in range ( k + 1 , n ) : if U [ kpivot , k ] != 0 : 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] = simplify((Ukk * U[i,j] - U[k,j]*Uik) / oldpivot) D = Ukk * U [ i , j ] - U [ k , j ] * Uik if len ( vars ) > 0 and D != S . Zero and not oldpivot . is_number : try : D , r = div ( Poly ( D , * vars ) , oldpivot ) except UnificationFailed : log . warn ( 'unification failed, trying direct division' ) D /= oldpivot else : D /= oldpivot # save if D . is_Atom : U [ i , j ] = D . as_expr ( ) else : if len ( vars ) > 0 : U [ i , j ] = D . as_expr ( ) else : U [ i , j ] = D . cancel ( ) U [ i , k ] = 0 oldpivot = Ukk DD [ n - 1 , n - 1 ] = oldpivot return P , L , DD , U