Trac #16803: Reimplement matrix_integer_dense using FLINT

In this ticket we reimplement all of matrix_integer_dense using FLINT
fmpz_mat_t.

The speed-up is substantial. With the new code:

{{{
sage: A = Matrix(ZZ,1000,1000,range(10^6))
sage: %time B = A*A
CPU times: user 883 ms, sys: 3.33 ms, total: 887 ms
Wall time: 888 ms
}}}

This takes more than 1 minute in the same computer using Sage 6.3.

URL: http://trac.sagemath.org/16803
Reported by: mmasdeu
Ticket author(s): Marc Masdeu
Reviewer(s): William Stein, Jeroen Demeyer

src/module_list.py

@@ -129,6 +129,7 @@ def uname_specific(name, value, alternative):
     Extension('sage.algebras.quatalg.quaternion_algebra_cython',
                sources = ['sage/algebras/quatalg/quaternion_algebra_cython.pyx'],
                language='c++',
+               depends = flint_depends,
                libraries = ["flint", "gmp", "gmpxx", "m", "stdc++", "ntl"]),
 
     ################################
@@ -699,7 +700,7 @@ def uname_specific(name, value, alternative):
               language="c++",
               include_dirs = [SAGE_INC + '/fplll'],
               extra_compile_args=["-DFPLLL_V3_COMPAT"],
-              depends = [SAGE_INC + "/fplll/fplll.h"]),
+              depends = [SAGE_INC + "/fplll/fplll.h"] + flint_depends),
 
     Extension('sage.libs.linbox.linbox',
               sources = ['sage/libs/linbox/linbox.pyx'],
@@ -750,7 +751,11 @@ def uname_specific(name, value, alternative):
 
     Extension('sage.libs.pari.pari_instance',
               sources = ["sage/libs/pari/pari_instance.pyx"],
-              libraries = ['pari', 'gmp']),
+              language="c",
+              extra_compile_args = ["-std=c99",  "-D_XPG6"],
+              libraries = [ 'pari',  'gmp','flint'],
+              depends =  flint_depends),
+
 
     Extension('sage.libs.ppl',
               sources = ['sage/libs/ppl.pyx', 'sage/libs/ppl_shim.cc'],
@@ -1052,13 +1057,12 @@ def uname_specific(name, value, alternative):
               sources = ['sage/matrix/matrix_integer_2x2.pyx'],
               libraries = ['gmp']),
 
-    # TODO -- change to use BLAS at some point.
     Extension('sage.matrix.matrix_integer_dense',
               sources = ['sage/matrix/matrix_integer_dense.pyx'],
-              extra_compile_args = m4ri_extra_compile_args,
-              depends = [SAGE_INC + '/m4ri/m4ri.h'],
+              extra_compile_args = ['-std=c99'] + m4ri_extra_compile_args,
               # order matters for cygwin!!
-              libraries = ['iml', 'pari', 'gmp', 'm', BLAS, BLAS2]),
+              libraries = ['iml', 'pari', 'ntl', 'gmp', 'm', 'flint', BLAS, BLAS2],
+              depends = [SAGE_INC + '/m4ri/m4ri.h']+ flint_depends),
 
     Extension('sage.matrix.matrix_integer_sparse',
               sources = ['sage/matrix/matrix_integer_sparse.pyx'],
@@ -1080,7 +1084,11 @@ def uname_specific(name, value, alternative):
 
     Extension('sage.matrix.matrix_modn_dense',
               sources = ['sage/matrix/matrix_modn_dense.pyx'],
-              libraries = ['gmp']),
+              language="c++",
+              extra_compile_args = ["-D_XPG6"]+ m4ri_extra_compile_args,
+              # order matters for cygwin!!
+              libraries = ['iml', 'pari', 'ntl', 'gmp', 'm', 'flint', BLAS, BLAS2],
+              depends = [SAGE_INC + '/m4ri/m4ri.h']+ flint_depends),
 
     Extension('sage.matrix.matrix_modn_dense_float',
               sources = ['sage/matrix/matrix_modn_dense_float.pyx'],
@@ -1092,7 +1100,7 @@ def uname_specific(name, value, alternative):
               sources = ['sage/matrix/matrix_modn_dense_double.pyx'],
               language="c++",
               libraries = ['linbox', 'givaro', 'mpfr', 'gmpxx', 'gmp', BLAS, BLAS2],
-              extra_compile_args = ['-DDISABLE_COMMENTATOR'] + givaro_extra_compile_args),
+              extra_compile_args = ["-D_XPG6" '-DDISABLE_COMMENTATOR']+ m4ri_extra_compile_args + givaro_extra_compile_args),
 
     Extension('sage.matrix.matrix_modn_sparse',
               sources = ['sage/matrix/matrix_modn_sparse.pyx'],
@@ -1107,7 +1115,11 @@ def uname_specific(name, value, alternative):
 
     Extension('sage.matrix.matrix_rational_dense',
               sources = ['sage/matrix/matrix_rational_dense.pyx'],
-              libraries = ['pari', 'gmp']),
+              language="c",
+              extra_compile_args = ["-std=c99",  "-D_XPG6"]+ m4ri_extra_compile_args,
+              # order matters for cygwin!!
+              libraries = ['iml', 'pari', 'ntl', 'gmp', 'm', 'flint', BLAS, BLAS2],
+              depends = [SAGE_INC + '/m4ri/m4ri.h']+ flint_depends),
 
     Extension('sage.matrix.matrix_rational_sparse',
               sources = ['sage/matrix/matrix_rational_sparse.pyx'],
@@ -1918,7 +1930,6 @@ def uname_specific(name, value, alternative):
                          flint_depends,
               libraries = ['flint', 'gmp', 'ratpoints']),
 
-
     Extension('sage.schemes.elliptic_curves.period_lattice_region',
               sources = ['sage/schemes/elliptic_curves/period_lattice_region.pyx'],
               include_dirs = numpy_include_dirs),
@@ -2095,13 +2106,13 @@ def uname_specific(name, value, alternative):
 
     Extension('sage.sat.solvers.satsolver',
               sources = ['sage/sat/solvers/satsolver.pyx']),
-              
+
     ################################
-    ## 
+    ##
     ## sage.schemes
     ##
     ################################
-              
+
     Extension('sage.schemes.projective.projective_morphism_helper',
               sources = ['sage/schemes/projective/projective_morphism_helper.pyx']),
     ]

src/sage/algebras/quatalg/quaternion_algebra_cython.pyx

@@ -88,21 +88,32 @@ def integral_matrix_and_denom_from_rational_quaternions(v, reverse=False):
     # Now fill in each row x of A, multiplying it by q = d/denom(x)
     cdef mpz_t q
     cdef mpz_t* row
+    cdef mpz_t tmp
     mpz_init(q)
+    mpz_init(tmp)
     for i in range(n):
         x = v[i]
         mpz_fdiv_q(q, d.value, x.d)
         if reverse:
-            mpz_mul(A._matrix[n-i-1][3], q, x.x)
-            mpz_mul(A._matrix[n-i-1][2], q, x.y)
-            mpz_mul(A._matrix[n-i-1][1], q, x.z)
-            mpz_mul(A._matrix[n-i-1][0], q, x.w)
+            mpz_mul(tmp, q, x.x)
+            A.set_unsafe_mpz(n-i-1,3,tmp)
+            mpz_mul(tmp, q, x.y)
+            A.set_unsafe_mpz(n-i-1,2,tmp)
+            mpz_mul(tmp, q, x.z)
+            A.set_unsafe_mpz(n-i-1,1,tmp)
+            mpz_mul(tmp, q, x.w)
+            A.set_unsafe_mpz(n-i-1,0,tmp)
         else:
-            mpz_mul(A._matrix[i][0], q, x.x)
-            mpz_mul(A._matrix[i][1], q, x.y)
-            mpz_mul(A._matrix[i][2], q, x.z)
-            mpz_mul(A._matrix[i][3], q, x.w)
+            mpz_mul(tmp, q, x.x)
+            A.set_unsafe_mpz(i,0,tmp)
+            mpz_mul(tmp, q, x.y)
+            A.set_unsafe_mpz(i,1,tmp)
+            mpz_mul(tmp, q, x.z)
+            A.set_unsafe_mpz(i,2,tmp)
+            mpz_mul(tmp, q, x.w)
+            A.set_unsafe_mpz(i,3,tmp)
     mpz_clear(q)
+    mpz_clear(tmp)
     return A, d
 
 def rational_matrix_from_rational_quaternions(v, reverse=False):
@@ -195,6 +206,8 @@ def rational_quaternions_from_integral_matrix_and_denom(A, Matrix_integer_dense
     a = Integer(A.invariants()[0])
     b = Integer(A.invariants()[1])
     cdef Py_ssize_t i, j
+    cdef mpz_t tmp
+    mpz_init(tmp)
 
     if reverse:
         rng = range(H.nrows()-1,-1,-1)
@@ -207,19 +220,28 @@ def rational_quaternions_from_integral_matrix_and_denom(A, Matrix_integer_dense
         mpz_set(x.a, a.value)
         mpz_set(x.b, b.value)
         if reverse:
-            mpz_init_set(x.x, H._matrix[i][3])
-            mpz_init_set(x.y, H._matrix[i][2])
-            mpz_init_set(x.z, H._matrix[i][1])
-            mpz_init_set(x.w, H._matrix[i][0])
+            H.get_unsafe_mpz(i,3,tmp)
+            mpz_init_set(x.x, tmp)
+            H.get_unsafe_mpz(i,2,tmp)
+            mpz_init_set(x.y, tmp)
+            H.get_unsafe_mpz(i,1,tmp)
+            mpz_init_set(x.z, tmp)
+            H.get_unsafe_mpz(i,0,tmp)
+            mpz_init_set(x.w, tmp)
         else:
-            mpz_init_set(x.x, H._matrix[i][0])
-            mpz_init_set(x.y, H._matrix[i][1])
-            mpz_init_set(x.z, H._matrix[i][2])
-            mpz_init_set(x.w, H._matrix[i][3])
+            H.get_unsafe_mpz(i,0,tmp)
+            mpz_init_set(x.x, tmp)
+            H.get_unsafe_mpz(i,1,tmp)
+            mpz_init_set(x.y, tmp)
+            H.get_unsafe_mpz(i,2,tmp)
+            mpz_init_set(x.z, tmp)
+            H.get_unsafe_mpz(i,3,tmp)
+            mpz_init_set(x.w, tmp)
         mpz_init_set(x.d, d.value)
         # WARNING -- we do *not* canonicalize the entries in the quaternion.  This is
         # I think _not_ needed for quaternion_element.pyx
         v.append(x)
+    mpz_clear(tmp)
     return v
 
 

src/sage/libs/flint/fmpz_mat.pxi

@@ -0,0 +1,32 @@
+include "sage/libs/flint/fmpz.pxi"
+include "sage/libs/flint/fmpz_poly.pxi"
+
+cdef extern from "flint/fmpz_mat.h":
+    ctypedef void* fmpz_mat_t
+    void fmpz_mat_init(fmpz_mat_t mat, unsigned long rows, unsigned long cols)
+    void fmpz_mat_init_set(fmpz_mat_t mat, const fmpz_mat_t src)
+    void fmpz_mat_set(fmpz_mat_t result, fmpz_mat_t mat)
+    void fmpz_mat_clear(fmpz_mat_t mat)
+    int fmpz_mat_print_pretty(const fmpz_mat_t mat)
+    fmpz_t fmpz_mat_entry(fmpz_mat_t mat, long i, long j)
+    void fmpz_mat_one(fmpz_mat_t mat)
+    void fmpz_mat_zero(fmpz_mat_t mat)
+    void fmpz_mat_neg(fmpz_mat_t f, fmpz_mat_t g)
+    void fmpz_mat_scalar_mul_si(fmpz_mat_t B, const fmpz_mat_t A, long c)
+    void fmpz_mat_scalar_mul_fmpz(fmpz_mat_t B, const fmpz_mat_t A, const fmpz_t c)
+    void fmpz_mat_mul(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B)
+    void fmpz_mat_sqr(fmpz_mat_t B, const fmpz_mat_t A)
+    void fmpz_mat_add(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B)
+    void fmpz_mat_sub(fmpz_mat_t C, const fmpz_mat_t A, const fmpz_mat_t B)
+    void fmpz_mat_pow(fmpz_mat_t C, const fmpz_mat_t A, unsigned long n)
+    int fmpz_mat_is_zero(const fmpz_mat_t mat)
+    void fmpz_mat_charpoly(fmpz_poly_t cp, const fmpz_mat_t mat)
+    void fmpz_mat_det(fmpz_t det, const fmpz_mat_t A)
+    int fmpz_mat_inv(fmpz_mat_t Ainv, fmpz_t den, const fmpz_mat_t A)
+    void fmpz_mat_transpose(fmpz_mat_t B, const fmpz_mat_t A)
+    long fmpz_mat_rref(fmpz_mat_t B, fmpz_t den, const fmpz_mat_t A)
+    long fmpz_mat_fflu(fmpz_mat_t B, fmpz_poly_t den, long *perm, const fmpz_mat_t A, int rank_check)
+    long fmpz_mat_rref_fraction_free(long * perm, fmpz_mat_t B, fmpz_t den, const fmpz_mat_t A)
+    long fmpz_mat_rank(const fmpz_mat_t A)
+    int fmpz_mat_solve(fmpz_mat_t X, fmpz_t den, const fmpz_mat_t A, const fmpz_mat_t B)
+    long fmpz_mat_nullspace(fmpz_mat_t B, const fmpz_mat_t A)

src/sage/libs/fplll/fplll.pxd

@@ -4,3 +4,12 @@ include "fplll.pxi"
 cdef class FP_LLL:
     cdef object fp_map
     cdef ZZ_mat[mpz_t] *_lattice
+
+cdef extern from "flint/fmpz.h":
+    ctypedef void* fmpz_t
+    void fmpz_get_mpz(mpz_t x, const fmpz_t f)
+    void fmpz_set_mpz(fmpz_t f, const mpz_t val)
+
+cdef extern from "flint/fmpz_mat.h":
+    ctypedef void* fmpz_mat_t
+    fmpz_t fmpz_mat_entry(fmpz_mat_t mat ,long i ,long j)

src/sage/libs/fplll/fplll.pyx

@@ -139,19 +139,20 @@ cdef class FP_LLL:
             ValueError: fpLLL cannot handle matrices with zero rows.
 
         """
-        cdef int i,j
-
-        cdef Z_NR[mpz_t] t
+        cdef unsigned long i,j
 
         if A._nrows == 0:
             raise ValueError('fpLLL cannot handle matrices with zero rows.')
 
         self._lattice = new ZZ_mat[mpz_t](A._nrows,A._ncols)
-
-        for i in range(A._nrows):
-            for j in range(A._ncols):
-                t.set(A._matrix[i][j])
-                self._lattice[0][i][j] = t
+        cdef mpz_t tmp
+        mpz_init(tmp)
+        for i from 0 <= i < A._nrows:
+            for j from 0 <= j < A._ncols:
+                A.get_unsafe_mpz(i,j,tmp)
+                # mpz_set(self._lattice[0][i][j],tmp)
+                self._lattice[0][i][j].set(tmp)
+        mpz_clear(tmp)
 
     def __dealloc__(self):
         """
@@ -1461,5 +1462,5 @@ cdef to_sage(ZZ_mat[mpz_t] *A):
 
     for i from 0 <= i < A.getRows():
         for j from 0 <= j < A.getCols():
-            mpz_set(B._matrix[i][j], A[0][i][j].getData())
+            B.set_unsafe_mpz(i,j,A[0][i][j].getData())
     return B

src/sage/libs/linkages/padics/mpz.pxi

@@ -77,7 +77,7 @@ cdef inline int ccmp(mpz_t a, mpz_t b, long prec, bint reduce_a, bint reduce_b,
     cdef int ans
     if reduce_a or reduce_b:
         mpz_sub(holder.value, a, b)
-        mpz_mod(holder.value, holder.value, prime_pow.pow_mpz_t_tmp(prec))
+        mpz_mod(holder.value, holder.value, &prime_pow.pow_mpz_t_tmp(prec)[0])
         return mpz_sgn(holder.value)
     else:
         ans = mpz_cmp(a,b)

src/sage/libs/pari/pari_instance.pxd

@@ -1,4 +1,6 @@
 include 'decl.pxi'
+include 'sage/libs/flint/fmpz.pxi'
+include 'sage/libs/flint/fmpz_mat.pxi'
 
 cimport sage.structure.parent_base
 cimport cython
@@ -27,8 +29,8 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase):
     cdef gen new_ref(self, GEN g, gen parent)
     cdef gen _empty_vector(self, long n)
     cdef long get_var(self, v)
-    cdef GEN _new_GEN_from_mpz_t_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc)
-    cdef GEN _new_GEN_from_mpz_t_matrix_rotate90(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc)
-    cdef gen integer_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf)
+    cdef GEN _new_GEN_from_fmpz_mat_t(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc)
+    cdef GEN _new_GEN_from_fmpz_mat_t_rotate90(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc)
+    cdef gen integer_matrix(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf)
     cdef GEN _new_GEN_from_mpq_t_matrix(self, mpq_t** B, Py_ssize_t nr, Py_ssize_t nc)
     cdef gen rational_matrix(self, mpq_t** B, Py_ssize_t nr, Py_ssize_t nc)

src/sage/libs/pari/pari_instance.pyx

@@ -889,7 +889,7 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase):
         """
         return objtogen(s)
 
-    cdef GEN _new_GEN_from_mpz_t_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc):
+    cdef GEN _new_GEN_from_fmpz_mat_t(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc):
         r"""
         Create a new PARI ``t_MAT`` with ``nr`` rows and ``nc`` columns
         from a ``mpz_t**``.
@@ -900,13 +900,16 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase):
         cdef GEN x
         cdef GEN A = zeromatcopy(nr, nc)
         cdef Py_ssize_t i, j
+        cdef mpz_t tmp
+        mpz_init(tmp)
         for i in range(nr):
             for j in range(nc):
-                x = self._new_GEN_from_mpz_t(B[i][j])
+                fmpz_get_mpz(tmp,fmpz_mat_entry(B,i,j))
+                x = self._new_GEN_from_mpz_t(tmp)
                 set_gcoeff(A, i+1, j+1, x)  # A[i+1, j+1] = x (using 1-based indexing)
         return A
 
-    cdef GEN _new_GEN_from_mpz_t_matrix_rotate90(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc):
+    cdef GEN _new_GEN_from_fmpz_mat_t_rotate90(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc):
         r"""
         Create a new PARI ``t_MAT`` with ``nr`` rows and ``nc`` columns
         from a ``mpz_t**`` and rotate the matrix 90 degrees
@@ -920,13 +923,16 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase):
         cdef GEN x
         cdef GEN A = zeromatcopy(nc, nr)
         cdef Py_ssize_t i, j
+        cdef mpz_t tmp
+        mpz_init(tmp)
         for i in range(nr):
             for j in range(nc):
-                x = self._new_GEN_from_mpz_t(B[i][nc-j-1])
+                fmpz_get_mpz(tmp,fmpz_mat_entry(B,i,nc-j-1))
+                x = self._new_GEN_from_mpz_t(tmp)
                 set_gcoeff(A, j+1, i+1, x)  # A[j+1, i+1] = x (using 1-based indexing)
         return A
 
-    cdef gen integer_matrix(self, mpz_t** B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf):
+    cdef gen integer_matrix(self, fmpz_mat_t B, Py_ssize_t nr, Py_ssize_t nc, bint permute_for_hnf):
         """
         EXAMPLES::
 
@@ -936,9 +942,9 @@ cdef class PariInstance(sage.structure.parent_base.ParentWithBase):
         pari_catch_sig_on()
         cdef GEN g
         if permute_for_hnf:
-            g = self._new_GEN_from_mpz_t_matrix_rotate90(B, nr, nc)
+            g = self._new_GEN_from_fmpz_mat_t_rotate90(B, nr, nc)
         else:
-            g = self._new_GEN_from_mpz_t_matrix(B, nr, nc)
+            g = self._new_GEN_from_fmpz_mat_t(B, nr, nc)
         return self.new_gen(g)
 
     cdef GEN _new_GEN_from_mpq_t_matrix(self, mpq_t** B, Py_ssize_t nr, Py_ssize_t nc):

src/sage/matrix/change_ring.pyx

@@ -39,7 +39,7 @@ def integer_to_real_double_dense(Matrix_integer_dense A):
                                          S, None, None, None)
     for i from 0 <= i < A._nrows:
         for j from 0 <= j < A._ncols:
-            M.set_unsafe_double(i, j, mpz_get_d(A._matrix[i][j]))
+            M.set_unsafe_double(i,j,A.get_unsafe_double(i,j))
     return M