Ticket #3684: mzd_kernel.patch

File mzd_kernel.patch, 6.1 KB (added by malb, 12 years ago)
  • sage/matrix/matrix_mod2_dense.pxd

    # HG changeset patch
    # User Martin Albrecht <malb@informatik.uni-bremen.de>
    # Date 1260898714 0
    # Node ID 7e790bf432473ff2e87ca3b13ce0f5723ae3da0e
    # Parent  964c2f4ce74db0417a771de0b0cfc951b1fab73c
    optimize computation of kernel for dense matrices over GF(2) (fixes #3684)
    
    diff -r 964c2f4ce74d -r 7e790bf43247 sage/matrix/matrix_mod2_dense.pxd
    a b  
    164164    # reduced row echelon form using PLUQ factorization
    165165    cdef long mzd_echelonize_pluq(mzd_t *A, int full)
    166166
     167    # reduced row echelon form using PLUQ factorization
     168    cdef mzd_t *mzd_kernel_left_pluq(mzd_t *A, int cutoff)
     169
    167170    ##################################
    168171    # Internal Functions for Debugging
    169172    ##################################
  • sage/matrix/matrix_mod2_dense.pyx

    diff -r 964c2f4ce74d -r 7e790bf43247 sage/matrix/matrix_mod2_dense.pyx
    a b  
    107107
    108108from sage.misc.misc import verbose, get_verbose, cputime
    109109
     110from sage.modules.free_module import VectorSpace
     111
    110112cdef extern from "gd.h":
    111113    ctypedef struct gdImagePtr "gdImagePtr":
    112114        pass
     
    15291531        else:
    15301532            return matrix_dense.Matrix_dense.density(self)
    15311533
    1532     def rank(self):
     1534    def rank(self, algorithm='lqup'):
    15331535        """
    15341536        Return the rank of this matrix.
    15351537
     1538        On average 'lqup' should be faster than 'm4ri' and hence it is
     1539        the default choice. However, for small - i.e. quite few
     1540        thousand rows & columns - and sparse matrices 'm4ri' might be
     1541        a better choice.
     1542
     1543        INPUT:
     1544
     1545        - ``algorithm`` - either "lqup" or "m4ri"
     1546
    15361547        EXAMPLE:
    15371548            sage: A = random_matrix(GF(2), 1000, 1000)
    15381549            sage: A.rank()
     
    15481559        if self._nrows == 0 or self._ncols == 0:
    15491560            return 0
    15501561        cdef mzd_t *A = mzd_copy(NULL, self._entries)
    1551         # for now we always call m4ri but in the future we will call a
    1552         # yet to be written mzd_rank(0) function.
    1553         r = mzd_echelonize_m4ri(A, 0, 0)
     1562        cdef mzp_t *P, *Q
     1563
     1564        if algorithm == 'lqup':
     1565            P = mzp_init(self._entries.nrows)
     1566            Q = mzp_init(self._entries.ncols)
     1567            r = mzd_lqup(A, P, Q, 0)
     1568            mzp_free(P)
     1569            mzp_free(Q)
     1570        elif algorithm == 'm4ri':
     1571            r = mzd_echelonize_m4ri(A, 0, 0)
     1572        else:
     1573            raise ValueError("Algorithm '%s' unknown."%algorithm)
    15541574        mzd_free(A)
    15551575        self.cache('rank', r)
    15561576        return r
    15571577
     1578    def right_kernel(self, algorithm='pluq'):
     1579        """
     1580        Return the right kernel of this matrix, as a vector
     1581        space. This is the space of vectors x such that ``self*x=0``.
     1582        A left kernel can be found with :meth:`left_kernel()` or just
     1583        :meth:`kernel`.
     1584           
     1585        INPUT:
     1586       
     1587        - ``algorithm`` - either "pluq" or "generic"
     1588       
     1589        By convention if self has 0 columns, the kernel is of dimension 0,
     1590        whereas the kernel is whole domain if self has 0 rows.
     1591       
     1592        .. note::
     1593
     1594           Preference is given to left kernels in that the generic method
     1595           name :meth:`kernel` returns a left kernel.  However most computations
     1596           of kernels are implemented as right kernels.
     1597       
     1598        EXAMPLES:
     1599       
     1600        A trivial right kernel::
     1601       
     1602            sage: A = MatrixSpace(GF(2), 2)([1,0,0,1])
     1603            sage: A.right_kernel()
     1604            Vector space of degree 2 and dimension 0 over Finite Field of size 2
     1605            Basis matrix:
     1606            []
     1607       
     1608        Right kernel of a zero matrix::
     1609       
     1610            sage: A = MatrixSpace(GF(2), 2)(0)
     1611            sage: A.right_kernel()
     1612            Vector space of degree 2 and dimension 2 over Finite Field of size 2
     1613            Basis matrix:
     1614            [1 0]
     1615            [0 1]
     1616       
     1617        Right kernel of a non-square matrix::
     1618       
     1619            sage: A = MatrixSpace(GF(2),2,3)(range(6))
     1620            sage: A.right_kernel()
     1621            Vector space of degree 3 and dimension 1 over Finite Field of size 2
     1622            Basis matrix:
     1623            [1 0 1]
     1624       
     1625        A non-trivial kernel computation::
     1626
     1627            sage: A = random_matrix(GF(2),1000,1010)
     1628            sage: A.right_kernel()
     1629            Vector space of degree 1010 and dimension 10 over Finite Field of size 2
     1630            Basis matrix:
     1631            10 x 1010 dense matrix over Finite Field of size 2
     1632        """
     1633        if algorithm == 'generic':
     1634            return matrix_dense.Matrix_dense.right_kernel(self)
     1635        if algorithm != 'pluq':
     1636            raise ValueError("Algorithm '%s' is unknown."%algorithm)
     1637
     1638        cdef Matrix_mod2_dense M
     1639        K = self.fetch('right_kernel')
     1640        if not K is None:
     1641            return K
     1642
     1643        R = self._base_ring
     1644        if self._ncols == 0:    # from a degree-0 space
     1645            V = VectorSpace(R, self._ncols)
     1646            Z = V.zero_subspace()
     1647            self.cache('right_kernel', Z)
     1648            return Z
     1649        elif self._nrows == 0:  # to a degree-0 space
     1650            Z = VectorSpace(R, self._ncols)
     1651            self.cache('right_kernel', Z)
     1652            return Z
     1653
     1654        cdef mzd_t *A = mzd_copy(NULL, self._entries)
     1655        cdef mzd_t *k = mzd_kernel_left_pluq(A, 0) # well, we don't
     1656                                                   # agree on the name
     1657                                                   # of this thing
     1658        mzd_free(A)
     1659        if k != NULL:
     1660            M = self.new_matrix(nrows = k.ncols, ncols = k.nrows)
     1661            mzd_transpose(M._entries, k)
     1662            mzd_free(k)
     1663            M = M.echelon_form()
     1664            basis = M.rows()
     1665        else:
     1666            basis = []
     1667
     1668        V = R**self._ncols
     1669        W = V.submodule(basis, check=False, already_echelonized=True)
     1670        self.cache('right_kernel', W)
     1671        return W
     1672
    15581673# Used for hashing
    15591674cdef int i, k
    15601675cdef unsigned long parity_table[256]