Add LinBox algorithm for right kernel of sparse matrix over the rationals

src/sage/libs/linbox/conversion.pxd

@@ -26,18 +26,22 @@ in the module ``linbox_flint_interface``.
 
 
 from libcpp.vector cimport vector as cppvector
+from libcpp.utility cimport move
 
 from sage.libs.gmp.mpz cimport mpz_set
+from sage.libs.gmp.mpq cimport mpq_get_num, mpq_get_den
 
-from sage.libs.linbox.givaro cimport Modular_uint64, ZRing, Integer
-from sage.libs.linbox.linbox cimport SparseMatrix_Modular_uint64, SparseMatrix_integer, DenseVector_integer
+from sage.libs.linbox.givaro cimport Modular_uint64, ZRing, Integer, QField, Rational
+from sage.libs.linbox.linbox cimport SparseMatrix_Modular_uint64, SparseMatrix_integer, SparseMatrix_rational, DenseVector_integer
 
 from sage.matrix.matrix_modn_sparse cimport Matrix_modn_sparse
 from sage.matrix.matrix_integer_sparse cimport Matrix_integer_sparse
+from sage.matrix.matrix_rational_sparse cimport Matrix_rational_sparse
 
 from sage.modules.vector_modn_sparse cimport c_vector_modint
 from sage.modules.vector_integer_dense cimport Vector_integer_dense
-from sage.modules.vector_integer_sparse cimport mpz_vector,  mpz_vector_get_entry, mpz_vector_set_entry
+from sage.modules.vector_integer_sparse cimport mpz_vector
+from sage.modules.vector_rational_sparse cimport mpq_vector
 
 ########################################
 # algorithm for solving linear systems #
@@ -135,7 +139,7 @@ cdef inline SparseMatrix_integer * new_linbox_matrix_integer_sparse(ZRing &ZZ, M
     r"""
     Return a new LinBox matrix from a Sage matrix.
 
-    Suc matrix has to be deallocated with a "del" statement.
+    Such a matrix has to be deallocated with a "del" statement.
 
     INPUT:
 
@@ -145,6 +149,44 @@ cdef inline SparseMatrix_integer * new_linbox_matrix_integer_sparse(ZRing &ZZ, M
     set_linbox_matrix_integer_sparse(A[0], m)
     return A
 
+##########################
+# matrix rational sparse #
+##########################
+
+cdef inline void set_linbox_matrix_rational_sparse(SparseMatrix_rational& A, Matrix_rational_sparse m) noexcept:
+    r"""
+    Set the entries of a LinBox matrix from a Sage matrix.
+
+    INPUT:
+
+    - ``A`` -- LinBox matrix
+    - ``m`` -- Sage matrix
+    """
+    cdef size_t i, j, k
+    cdef mpq_vector * v
+    cdef Integer num, denom
+    for i in range(<size_t> m._nrows):
+        v = m._matrix + i
+        for k in range(<size_t> v.num_nonzero):
+            j = v.positions[k]
+            mpq_get_num(num.get_mpz(), v.entries[k])
+            mpq_get_den(denom.get_mpz(), v.entries[k])
+            A.setEntry(i, j, Rational(move(num), move(denom), 0))
+
+cdef inline SparseMatrix_rational * new_linbox_matrix_rational_sparse(QField &QQ, Matrix_rational_sparse m) noexcept:
+    r"""
+    Return a new LinBox matrix from a Sage matrix.
+
+    Such a matrix has to be deallocated with a "del" statement.
+
+    INPUT:
+
+    - ``m`` -- Sage matrix
+    """
+    cdef SparseMatrix_rational * A = new SparseMatrix_rational(QQ, <size_t> m._nrows, <size_t> m._ncols)
+    set_linbox_matrix_rational_sparse(A[0], m)
+    return A
+
 ########################
 # vector integer dense #
 ########################

src/sage/libs/linbox/givaro.pxd

@@ -25,6 +25,23 @@ cdef extern from "gmp++/gmp++.h" namespace "Givaro":
         mpz_ptr get_mpz()
         mpz_srcptr get_mpz_const()
 
+cdef extern from "givaro/givrational.h" namespace "Givaro":
+    cdef cppclass Rational:
+        Rational()
+        Rational(int32_t)
+        Rational(int64_t)
+        Rational(uint32_t)
+        Rational(uint64_t)
+        Rational(Integer&)
+        Rational(int32_t, int32_t)
+        Rational(int64_t, int64_t)
+        Rational(uint32_t, uint32_t)
+        Rational(uint64_t, uint64_t)
+        Rational(Integer&, Integer&, int)
+
+        Integer nume()
+        Integer deno()
+
 cdef extern from "givaro/givcategory.h" namespace "Givaro":
     cdef cppclass Sporadic:
         pass
@@ -41,6 +58,14 @@ cdef extern from "givaro/zring.h":
         Element one
         Element mone
 
+cdef extern from "givaro/qfield.h":
+    ## template<class RatElement> class QField
+    cdef cppclass QField "Givaro::QField<Givaro::Rational>":
+        ctypedef Rational Element
+        Element one
+        Element mOne
+        Element zero
+
 cdef extern from "givaro/modular-integral.h":
     cdef cppclass Modular_uint64 "Givaro::Modular<uint64_t>":
         ctypedef uint64_t Element

src/sage/libs/linbox/linbox.pxd

@@ -7,6 +7,7 @@
 
 from libc.stdint cimport uint32_t, uint64_t
 from libcpp.vector cimport vector as cppvector
+from libcpp.pair cimport pair
 
 from sage.libs.linbox.givaro cimport *
 
@@ -46,6 +47,12 @@ cdef extern from "linbox/matrix/dense-matrix.h":
 
         ostream& write(ostream&)
 
+cdef extern from "linbox/vector/sparse.h":
+    cdef cppclass Sparse_Vector_rational "LinBox::Sparse_Vector<Givaro::Rational>":
+        ctypedef pair[unsigned int, Rational] Element
+        size_t size()
+        Element& operator[](size_t i)
+
 cdef extern from "linbox/matrix/sparse-matrix.h":
     ## template<class _Field, class _Storage = SparseMatrixFormat::SparseSeq >
     ## class SparseMatrix ;
@@ -74,6 +81,19 @@ cdef extern from "linbox/matrix/sparse-matrix.h":
 
         ostream& write(ostream&)
 
+    cdef cppclass SparseMatrix_rational "LinBox::SparseMatrix<Givaro::QField<Givaro::Rational>>":
+        ctypedef QField Field
+        ctypedef Rational Element
+        SparseMatrix_rational(Field &F, size_t m, size_t n)
+        size_t rowdim()
+        size_t coldim()
+        void setEntry(size_t i, size_t j, Element &a)
+        Element &getEntry(size_t i, size_t j)
+        Field& field()
+        Sparse_Vector_rational& getRow(size_t i)
+
+        ostream& write(ostream&)
+
 cdef extern from "linbox/polynomial/dense-polynomial.h":
     ## template<class Field>
     ## class DensePolynomial : public Givaro::Poly1FactorDom<Field, Givaro::Dense>::Element
@@ -158,6 +178,12 @@ cdef extern from "linbox/algorithms/gauss.h":
                                              unsigned long Ni,
                                              unsigned long Nj)
 
+    cdef cppclass GaussDomain_rational "LinBox::GaussDomain<Givaro::QField<Givaro::Rational>>":
+        ctypedef QField Field
+        ctypedef Rational Element
+        GaussDomain_rational(Field &)
+        SparseMatrix_rational& nullspacebasisin(SparseMatrix_rational &X, SparseMatrix_rational &A)
+
 cdef extern from "linbox/solutions/echelon.h" namespace "LinBox":
     size_t rowEchelon (DenseMatrix_Modular_float&, const DenseMatrix_Modular_float&)
     size_t rowEchelonize (DenseMatrix_Modular_float&)

src/sage/matrix/matrix2.pyx

@@ -4296,9 +4296,11 @@ cdef class Matrix(Matrix1):
           - ``'generic'`` -- naive algorithm usable for matrices over any field
           - ``'flint'`` -- FLINT library code for matrices over the rationals
             or the integers
+          - ``'linbox'`` -- LinBox library code for sparse matrices over the
+            rationals
           - ``'pari'`` -- PARI library code for matrices over number fields
             or the integers
-          - ``'padic'`` -- padic algorithm from IML library for matrices
+          - ``'padic'`` -- `p`-adic algorithm from the IML library for matrices
             over the rationals and integers
           - ``'pluq'`` -- PLUQ matrix factorization for matrices mod 2
 
@@ -4348,12 +4350,17 @@ cdef class Matrix(Matrix1):
 
         Over the Rational Numbers:
 
-        Kernels are computed by the IML library in
-        :meth:`~sage.matrix.matrix_rational_dense.Matrix_rational_dense._right_kernel_matrix`.
-        Setting the `algorithm` keyword to 'default', 'padic' or unspecified
-        will yield the same result, as there is no optional behavior.
-        The 'computed' format of the basis vectors are exactly the negatives
-        of the vectors in the 'pivot' format. ::
+        The ``algorithm`` keyword can be set to ``'default'`` or ``'padic'`` to
+        use the `p`-adic algorithm from the IML library (for dense or sparse
+        matrices), or to ``'linbox'`` to use the LinBox library (only for
+        sparse matrices).
+        Kernels are computed by the IML library for dense matrices in
+        :meth:`~sage.matrix.matrix_rational_dense.Matrix_rational_dense._right_kernel_matrix`,
+        and sparse matrices are first converted to dense matrices.
+        Kernels are computed by the LinBox library for sparse matrices in
+        :meth:`~sage.matrix.matrix_rational_sparse.Matrix_rational_sparse._right_kernel_matrix_linbox`.
+        The 'computed' format of the basis vectors returned by IML are exactly
+        the negatives of the vectors in the 'pivot' format. ::
 
             sage: A = matrix(QQ, [[1, 0, 1, -3, 1],
             ....:                 [-5, 1, 0, 7, -3],
@@ -4418,6 +4425,11 @@ cdef class Matrix(Matrix1):
             sage: set_verbose(0)
             sage: D == S
             True
+            sage: K = B.right_kernel_matrix(algorithm='linbox', basis='computed'); K
+            [ 1 -2  2  1  0]
+            [-1 -2  0  0  1]
+            sage: B*K.transpose() == zero_matrix(QQ, 4, 2, sparse=True)
+            True
 
         Over Number Fields:
 
@@ -4839,14 +4851,19 @@ cdef class Matrix(Matrix1):
         algorithm = kwds.pop('algorithm', None)
         if algorithm is None:
             algorithm = 'default'
-        elif algorithm not in ['default', 'generic', 'flint', 'pari', 'padic', 'pluq']:
+        elif algorithm not in ['default', 'generic', 'flint', 'linbox', 'pari', 'padic', 'pluq']:
             raise ValueError("matrix kernel algorithm '%s' not recognized" % algorithm)
         elif algorithm == 'padic' and not isinstance(R, (IntegerRing_class,
                                                          RationalField)):
             raise ValueError("'padic' matrix kernel algorithm only available over the rationals and the integers, not over %s" % R)
         elif algorithm == 'flint' and not isinstance(R, (IntegerRing_class,
                                                          RationalField)):
             raise ValueError("'flint' matrix kernel algorithm only available over the rationals and the integers, not over %s" % R)
+        elif algorithm == 'linbox' and not isinstance(self, sage.matrix.matrix_rational_sparse.Matrix_rational_sparse):
+            if isinstance(R, RationalField):
+                raise ValueError("'linbox' matrix kernel algorithm only available for sparse matrices over the rationals")
+            else:
+                raise ValueError("'linbox' matrix kernel algorithm only available over the rationals, not over %s" % R)
         elif algorithm == 'pari' and not (isinstance(R, (IntegerRing_class, NumberField)) and not isinstance(R, RationalField)):
             raise ValueError("'pari' matrix kernel algorithm only available over non-trivial number fields and the integers, not over %s" % R)
         elif algorithm == 'generic' and R not in _Fields:
@@ -4996,9 +5013,11 @@ cdef class Matrix(Matrix1):
           - ``'generic'`` -- naive algorithm usable for matrices over any field
           - ``'flint'`` -- FLINT library code for matrices over the rationals
             or the integers
+          - ``'linbox'`` -- LinBox library code for sparse matrices over the
+            rationals
           - ``'pari'`` -- PARI library code for matrices over number fields
             or the integers
-          - ``'padic'`` -- padic algorithm from IML library for matrices
+          - ``'padic'`` -- `p`-adic algorithm from the IML library for matrices
             over the rationals and integers
           - ``'pluq'`` -- PLUQ matrix factorization for matrices mod 2
 
@@ -5233,7 +5252,7 @@ cdef class Matrix(Matrix1):
             Quaternion Algebra (-1, -1) with base ring Rational Field
 
         Sparse matrices, over the rationals and the integers,
-        use the same routines as the dense versions. ::
+        use the same routines as the dense versions by default. ::
 
             sage: A = matrix(ZZ, [[0, -1, 1, 1, 2],
             ....:                 [1, -2, 0, 1, 3],
@@ -5365,9 +5384,11 @@ cdef class Matrix(Matrix1):
           - ``'generic'`` -- naive algorithm usable for matrices over any field
           - ``'flint'`` -- FLINT library code for matrices over the rationals
             or the integers
+          - ``'linbox'`` -- LinBox library code for sparse matrices over the
+            rationals
           - ``'pari'`` -- PARI library code for matrices over number fields
             or the integers
-          - ``'padic'`` -- padic algorithm from IML library for matrices
+          - ``'padic'`` -- `p`-adic algorithm from the IML library for matrices
             over the rationals and integers
           - ``'pluq'`` -- PLUQ matrix factorization for matrices mod 2