Ticket #4260: trac_4260_echelonformdomain.patch

File trac_4260_echelonformdomain.patch, 7.1 KB (added by malb, 9 years ago)
  • new file sage/libs/linbox/echelonform.pxd

    # HG changeset patch
    # User Martin Albrecht <martinralbrecht@googlemail.com>
    # Date 1317130575 -3600
    # Node ID d2e28530d98ce4c4bccd40682cec7f9f81d1a90f
    # Parent  69ac0934111197aa81be803516516d71d8077f3e
    #4260: use EchelonFormDomain, needs more memory but is twice as fast
    
    diff -r 69ac09341111 -r d2e28530d98c sage/libs/linbox/echelonform.pxd
    - +  
     1cdef extern from "linbox/matrix/blas-matrix.h" namespace "LinBox":
     2    cdef cppclass BlasMatrix[T]:
     3        BlasMatrix(size_t nrows, size_t ncols)
     4        void setEntry(size_t i, size_t j, T)
     5        T &getEntry(size_t i, size_t j)
     6
     7from sage.libs.linbox.modular cimport ModDoubleField, ModDoubleFieldElement, ModFloatField, ModFloatFieldElement
     8
     9cdef extern from "linbox/algorithms/echelon-form.h":
     10    cdef cppclass EchelonFormDomainDouble "LinBox::EchelonFormDomain<LinBox::Modular<double> >":
     11        EchelonFormDomainDouble(ModDoubleField)
     12        int rowReducedEchelon(BlasMatrix[ModDoubleFieldElement], BlasMatrix[ModDoubleFieldElement])
     13
     14    cdef cppclass EchelonFormDomainFloat "LinBox::EchelonFormDomain<LinBox::Modular<float> >":
     15        EchelonFormDomainFloat(ModFloatField)
     16        int rowReducedEchelon(BlasMatrix[ModFloatFieldElement], BlasMatrix[ModFloatFieldElement])
     17 
  • sage/libs/linbox/fflas.pxd

    diff -r 69ac09341111 -r d2e28530d98c sage/libs/linbox/fflas.pxd
    a b  
    11from modular cimport ModDoubleField, ModFloatField, ModDoubleFieldElement, ModFloatFieldElement
    22
    3 ctypedef unsigned long mod_int
    4 
    53cdef extern from "linbox/fflas/fflas.h" namespace "std":
    64    cdef cppclass vector[T]:
    75        cppclass iterator:
  • sage/matrix/matrix_modn_dense_double.pyx

    diff -r 69ac09341111 -r d2e28530d98c sage/matrix/matrix_modn_dense_double.pyx
    a b  
    2020# randstate in template needs this
    2121include '../ext/random.pxi'
    2222
     23from sage.libs.linbox.echelonform cimport  BlasMatrix
     24
    2325from sage.libs.linbox.modular cimport ModDoubleField as ModField, ModDoubleFieldElement as ModFieldElement
     26from sage.libs.linbox.echelonform cimport EchelonFormDomainDouble as EchelonFormDomain
    2427
    2528from sage.libs.linbox.fflas cimport ModDouble_fgemm as Mod_fgemm, ModDouble_fgemv as Mod_fgemv, \
    2629    ModDoubleDet as ModDet, \
  • sage/matrix/matrix_modn_dense_float.pyx

    diff -r 69ac09341111 -r d2e28530d98c sage/matrix/matrix_modn_dense_float.pyx
    a b  
    2020# randstate in template needs this
    2121include '../ext/random.pxi'
    2222
     23from sage.libs.linbox.echelonform cimport BlasMatrix
     24
    2325from sage.libs.linbox.modular cimport ModFloatField as ModField, ModFloatFieldElement as ModFieldElement
     26from sage.libs.linbox.echelonform cimport EchelonFormDomainFloat as EchelonFormDomain
    2427
    2528from sage.libs.linbox.fflas cimport ModFloat_fgemm as Mod_fgemm, ModFloat_fgemv as Mod_fgemv, \
    2629        ModFloatDet as ModDet, \
  • sage/matrix/matrix_modn_dense_template.pxi

    diff -r 69ac09341111 -r d2e28530d98c sage/matrix/matrix_modn_dense_template.pxi
    a b  
    156156    del F
    157157    return r, pivots
    158158
     159cdef inline linbox_echelonize_efd(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols):
     160    cdef ModField *F = new ModField(<long>modulus)
     161    cdef EchelonFormDomain *EF = new EchelonFormDomain(F[0])
     162    cdef BlasMatrix[ModFieldElement] *A = new BlasMatrix[ModFieldElement](nrows, ncols)
     163    cdef BlasMatrix[ModFieldElement] *E = new BlasMatrix[ModFieldElement](nrows, ncols)
     164
     165    cdef Py_ssize_t i,j
     166
     167    # TODO: can we avoid this copy?
     168    for i in range(nrows):
     169        for j in range(ncols):
     170            A.setEntry(i, j, <ModFieldElement>entries[i*ncols+j])
     171
     172    cdef int r = EF.rowReducedEchelon(E[0], A[0])
     173    for i in range(nrows):
     174        for j in range(ncols):
     175            entries[i*ncols+j] = <celement>E.getEntry(i,j)
     176
     177    cdef Py_ssize_t ii = 0
     178    cdef list pivots = []
     179    for i in range(r):
     180        for j in range(ii,ncols):
     181            if entries[i*ncols+j] == 1:
     182                pivots.append(j)
     183                ii = j+1
     184                break
     185
     186    # TODO: recover pivots
     187    del F, A, E, EF
     188    return r, pivots
     189
    159190cdef inline celement *linbox_copy(celement modulus, celement *entries,  Py_ssize_t nrows, Py_ssize_t ncols) except NULL:
    160191    """
    161192    Create a copy of the entries array.
     
    15741605
    15751606        - ``algorithm``
    15761607
    1577           - ``linbox`` - uses the LinBox library
     1608          - ``linbox`` - uses the LinBox library (``EchelonFormDomain`` implementation, default)
     1609
     1610          - ``linbox_noefd`` - uses the LinBox library (FFPACK directly, less memory but slower)
    15781611
    15791612          - ``gauss`` - uses a custom slower `O(n^3)` Gauss
    15801613            elimination implemented in Sage.
     
    17611794            raise NotImplementedError("Echelon form not implemented over '%s'."%self.base_ring())
    17621795
    17631796        if algorithm == 'linbox':
    1764             self._echelonize_linbox()
    1765 
     1797            self._echelonize_linbox(efd=True)
     1798        elif algorithm == 'linbox_noefd':
     1799            self._echelonize_linbox(efd=False)
    17661800        elif algorithm == 'gauss':
    17671801            self._echelon_in_place_classical()
    17681802
    17691803        elif algorithm == 'all':
    17701804            A = self.__copy__()
    1771             self._echelonize_linbox()
     1805            B = self.__copy__()
     1806            self._echelonize_linbox(efd=True)
    17721807            A._echelon_in_place_classical()
    1773             if A != self:
     1808            B._echelonize_linbox(efd=False)
     1809            if A != self or A != B:
    17741810                raise ArithmeticError("Bug in echelon form.")
    17751811        else:
    17761812            raise ValueError("Algorithm '%s' not known"%algorithm)
    17771813
    1778     def _echelonize_linbox(self):
     1814    def _echelonize_linbox(self, efd=True):
    17791815        """
    17801816        Puts ``self`` in row echelon form using LinBox.
    17811817
    17821818        This function is called by echelonize if
    17831819        ``algorithm='linbox'``.
    17841820
     1821        INPUT:
     1822
     1823        - ``efd`` - if ``True`` LinBox's ``EchelonFormDomain``
     1824          implementation is used, which is faster than the direct
     1825          ``LinBox::FFPACK`` implementation, since the latter also
     1826          computes the transformation matrix (which we
     1827          ignore). However, ``efd=True`` uses more memory than FFLAS
     1828          directly (default=``True``)
     1829
    17851830        EXAMPLES::
    17861831
    17871832            sage: A = random_matrix(GF(7), 10, 20); A
     
    18121857        self.clear_cache()
    18131858
    18141859        t = verbose('Calling echelonize mod %d.'%self.p)
    1815         r, pivots = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols)
     1860        if efd:
     1861            r, pivots = linbox_echelonize_efd(self.p, self._entries, self._nrows, self._ncols)
     1862        else:
     1863            r, pivots = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols)
    18161864        verbose('done with echelonize',t)
    18171865
    18181866        self.cache('in_echelon_form',True)