Ticket #11195: trac_11195-LU-decomposition-exact.patch

File trac_11195-LU-decomposition-exact.patch, 3.1 KB (added by rbeezer, 10 years ago)
  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1302757190 25200
    # Node ID d32a6ce9ab25a66cc097a2af2f81e082b6fd754b
    # Parent  461868191947707e8c6a764078c3bec7274088a6
    imported patch lu-exact
    
    diff -r 461868191947 -r d32a6ce9ab25 sage/matrix/matrix2.pyx
    a b  
    65646564            self.cache('cholesky', L)
    65656565        return L
    65666566
     6567    def LU(self, pivot='partial', format='plu'):
     6568        r"""
     6569        """
     6570        cdef Py_ssize_t m, n, d, i, j, k, p, max_location
     6571        cdef Matrix M
     6572        R = self.base_ring()
     6573        if not R.is_field():
     6574            try:
     6575                R = R.fraction_field()
     6576                M = self.change_ring(R)
     6577            except:
     6578                raise TypeError('matrix entries need to be in a field, not {0}'.format(R))
     6579        else:
     6580            M = self.__copy__()
     6581        if pivot == 'partial':
     6582            try:
     6583                abs(R.an_element())
     6584                ordered = True
     6585            except:
     6586                raise TypeError("cannot take absolute value of matrix entries, try 'pivot=None'")
     6587        elif pivot == None:
     6588            ordered = False
     6589        else:
     6590            raise ValueError("pivot strategy must be 'partial' or None, not {0}".format(pivot))
     6591        m, n = self._nrows, self._ncols
     6592        d = min(m, n)
     6593        perm = range(m)
     6594        zero = R(0)
     6595        for k in range(d):
     6596            max_location = -1
     6597            max_entry = zero
     6598            if ordered:
     6599                for i in range(k,m):
     6600                    entry = abs(M.get_unsafe(i,k))
     6601                    if entry > max_entry:
     6602                        max_location = i
     6603                        max_entry = entry
     6604            else:
     6605                for i in range(k,m):
     6606                    if M.get_unsafe(i,k) != zero:
     6607                        max_location = i
     6608                        break
     6609            if max_location != -1:
     6610                perm[k], perm[max_location] = perm[max_location], perm[k]
     6611                M.swap_rows(k, max_location)
     6612                for j in range(k+1, m):
     6613                    scale = -M.get_unsafe(j,k)/M.get_unsafe(k,k)
     6614                    M.set_unsafe(j,k, -scale)
     6615                    for p in range(k+1,n):
     6616                        M.set_unsafe(j,p, M.get_unsafe(j,p) + scale*M.get_unsafe(k,p))
     6617
     6618        if format == 'compact':
     6619            return perm, M
     6620        elif format == 'plu':
     6621            import sage.matrix.constructor
     6622            import sage.combinat.permutation
     6623            perm = [perm[i]+1 for i in range(m)]
     6624            P = sage.combinat.permutation.Permutation(perm).to_matrix()
     6625            L = sage.matrix.constructor.identity_matrix(R, m)
     6626            for i in range(1, m):
     6627                for k in range(min(i,d)):
     6628                    L[i,k] = M[i,k]
     6629                    M[i,k] = zero
     6630            return P, L, M
     6631        else:
     6632            raise ValueError("format must be 'plu' or 'compact', not {0}".format(format))
     6633
     6634
     6635
     6636
     6637
     6638
     6639
     6640
     6641
     6642
     6643
     6644
    65676645    def hadamard_bound(self):
    65686646        r"""
    65696647        Return an int n such that the absolute value of the determinant of