enhance sparse integer matrix with linbox (+ cleaning)

[vdelecroix]

LinBox has algorithms to work with sparse matrices. This ticket provides the necessary interface for:

• solving sparse systems over integers
• computing determinant and rank
• computing minimal and characteristic polynomials

Some cleaning had to be done as a first step (used to be #26178).

LinBox issues

• ​Issue 118 (determinant hangs)
• ​Issue 138 (Compilation failure on OSX)
• ​Issue 144 (row index: unsigned long versus size_t)

(see also #22872 for the linbox task ticket)

[kedlaya]

Just to put in a plug for this, it would be very useful for various aspects of number theory, particularly dealing with Hecke operators which tend to be represented by sparse integer matrices.

[vdelecroix]

And combinatorics and geomtry and ... just waiting for upstream ​github issue 56.

[roed]

That issue has been closed. Are we waiting on a release?

[vdelecroix]

Indeed, Clément wrote on the github issue

Note that we're also planning on releasing all 3 libraries
soon (within a few weeks), and I'll create a sage ticket
for updating the spkg.

[cpernet]

FYI LinBox upgrade ticket is #24214. Current ticket already depends on it by transitivity with #23158.

[vdelecroix]

New commits:

​eb54e4b	24544: cleaning linbox headers
​9d787d3	24544: fix matrix_modn_dense
​99ed944	24544: fix matrix_modn_sparse
​13f4de4	24544: fix documentation
​a1e4933	24544: do not compile the sage extension support anymore
​956610a	24544: remove linbox GFq from benchmarks
​7e92d8e	23214: sparse integer matrices and LinBox
​2712020	23214: empty conversion.pyx to make cython happy

[vdelecroix]

This hangs forever

sage: m = diagonal_matrix(ZZ, 68, [2]*66 + [1,1])
sage: type(m)
<type 'sage.matrix.matrix_integer_sparse.Matrix_integer_sparse'>
sage: m.det()

(from matrix/matrix_integer_dense_hnf.py)

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​74ccd06	23214: sparse integer matrices and LinBox
​1cd8037	23214: don't use LinBox for determinant

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​5b7b84a	23214: sparse integer matrices and LinBox
​8423615	23214: don't use LinBox for determinant

[jdemeyer]

sage: Matrix(IntegerModRing(16), 3, 3, lambda i,j: i+j, sparse=True).determinant()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-2d0c9b09dcc5> in <module>()
----> 1 Matrix(IntegerModRing(Integer(16)), Integer(3), Integer(3), lambda i,j: i+j, sparse=True).determinant()

/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/matrix/matrix_modn_sparse.pyx in sage.matrix.matrix_modn_sparse.Matrix_modn_sparse.determinant (build/cythonized/sage/matrix/matrix_modn_sparse.cpp:10523)()
    924 
    925         if algorithm is None or algorithm == "linbox":
--> 926             r, d = self._rank_det_linbox()
    927             self.cache('rank', r)
    928             self.cache('det', d)

/usr/local/src/sage-config/local/lib/python2.7/site-packages/sage/matrix/matrix_modn_sparse.pyx in sage.matrix.matrix_modn_sparse.Matrix_modn_sparse._rank_det_linbox (build/cythonized/sage/matrix/matrix_modn_sparse.cpp:9553)()
    766 
    767         if not is_prime(self.p):
--> 768             raise TypeError("only GF(p) supported via LinBox")
    769 
    770         cdef Modular_int64 * F = new Modular_int64(self.p)

TypeError: only GF(p) supported via LinBox

[vdelecroix]

update milestone 8.3 -> 8.4

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​2c157ff	24544: cleaning linbox headers
​b037c87	24544: fix matrix_modn_dense
​a12f0a1	24554: fix matrix modn sparse
​52d8a34	24544: fix documentation
​8acfbc9	24544: remove linbox GFq from benchmarks
​43babde	23214: use linbox in sparse integer matrices

[vdelecroix]

patchbot is happy :-)

[alexjbest]

Oops, I got this confused with the other needs review linbox ticket where compiling on OSX was an issue, I'll try them both out together now I guess.

[vdelecroix]

@alexjbest: this ticket strictly contains the commit of the other one #24544.

[tscrim]

Small typo (missing leading colon): meth:`sage.matrix.matrix_sparse.Matrix_sparse.charpoly` Otherwise I am happy with the current state (and will assume there are no failures on OSX).

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1d8d266

23214: small typo

[tscrim]

LGTM. Thanks. Alex, if you end up hitting an error during your testing, please revert.

[vbraun]

Doesn't build on OSX:

[sagelib-8.4.beta3] In file included from build/cythonized/sage/matrix/matrix_integer_sparse.cpp:708:
[sagelib-8.4.beta3] In file included from /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/matrix/dense-matrix.h:79:
[sagelib-8.4.beta3] In file included from /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/matrix/densematrix/blas-matrix.h:44:
[sagelib-8.4.beta3] /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/vector/blas-vector.h:306:11: error: member reference base type 'const long' is not a structure or union
[sagelib-8.4.beta3]                         _size(V.size()),_1stride(1),_rep(V.size(), F.zero),_ptr(&_rep[0]),_field(&F)
[sagelib-8.4.beta3]                               ~^~~~~
[sagelib-8.4.beta3] build/cythonized/sage/matrix/matrix_integer_sparse.cpp:10932:21: note: in instantiation of function template specialization 'LinBox::BlasVector<Givaro::ZRing<Givaro::Integer>, std::__1::vector<Givaro::Integer, std::__1::allocator<Givaro::Integer> > >::BlasVector<long>' requested here
[sagelib-8.4.beta3]   __pyx_v_res = new LinBox::DenseVector<Givaro::ZRing<Givaro::Integer>>(__pyx_v_givZZ, __pyx_v_self->__pyx_base.__pyx_base.__pyx_base.__pyx_base.__pyx_base._ncols);
[sagelib-8.4.beta3]                     ^
[sagelib-8.4.beta3] In file included from build/cythonized/sage/matrix/matrix_integer_sparse.cpp:708:
[sagelib-8.4.beta3] In file included from /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/matrix/dense-matrix.h:79:
[sagelib-8.4.beta3] In file included from /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/matrix/densematrix/blas-matrix.h:44:
[sagelib-8.4.beta3] /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/vector/blas-vector.h:306:38: error: member reference base type 'const long' is not a structure or union
[sagelib-8.4.beta3]                         _size(V.size()),_1stride(1),_rep(V.size(), F.zero),_ptr(&_rep[0]),_field(&F)
[sagelib-8.4.beta3]                                                          ~^~~~~
[sagelib-8.4.beta3] /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/vector/blas-vector.h:306:61: warning: field '_rep' is uninitialized when used here [-Wuninitialized]
[sagelib-8.4.beta3]                         _size(V.size()),_1stride(1),_rep(V.size(), F.zero),_ptr(&_rep[0]),_field(&F)
[sagelib-8.4.beta3]                                                                                  ^
[sagelib-8.4.beta3] /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/vector/blas-vector.h:129:13: error: type 'long' cannot be used prior to '::' because it has no members
[sagelib-8.4.beta3]                         typename OtherVector::const_iterator jt = V.begin();
[sagelib-8.4.beta3]                                  ^
[sagelib-8.4.beta3] /Users/buildslave-sage/slave/sage_git/build/local/include/linbox/vector/blas-vector.h:311:4: note: in instantiation of function template specialization 'LinBox::BlasVector<Givaro::ZRing<Givaro::Integer>, std::__1::vector<Givaro::Integer, std::__1::allocator<Givaro::Integer> > >::createBlasVector<long>' requested here
[sagelib-8.4.beta3]                         createBlasVector(V);
[sagelib-8.4.beta3]                         ^
[sagelib-8.4.beta3] build/cythonized/sage/matrix/matrix_integer_sparse.cpp:10932:21: note: in instantiation of function template specialization 'LinBox::BlasVector<Givaro::ZRing<Givaro::Integer>, std::__1::vector<Givaro::Integer, std::__1::allocator<Givaro::Integer> > >::BlasVector<long>' requested here
[sagelib-8.4.beta3]   __pyx_v_res = new LinBox::DenseVector<Givaro::ZRing<Givaro::Integer>>(__pyx_v_givZZ, __pyx_v_self->__pyx_base.__pyx_base.__pyx_base.__pyx_base.__pyx_base._ncols);
[sagelib-8.4.beta3]                     ^
[sagelib-8.4.beta3] 8 warnings and 3 errors generated.
[sagelib-8.4.beta3] error: command 'gcc' failed with exit status 1

[vdelecroix]

With the relevant lines from matrix_integer_sparse.cpp:

  /* "sage/matrix/matrix_integer_sparse.pyx":1074
 *         cdef linbox.SparseMatrix_integer * A = new_linbox_matrix_integer_sparse(givZZ, self)
 *         cdef linbox.DenseVector_integer * b = new_linbox_vector_integer_dense(givZZ, v)
 *         cdef linbox.DenseVector_integer * res = new linbox.DenseVector_integer(givZZ, self._ncols)             # <<<<<<<<<<<<<<
 *         cdef givaro.Integer D
 * 
 */
  __pyx_v_res = new LinBox::DenseVector<Givaro::ZRing<Givaro::Integer>>(__pyx_v_givZZ, __pyx_v_self->__pyx_base.__pyx_base.__pyx_base.__pyx_base.__pyx_base._ncols);

[vklein]

This compile failure can be reproduced outside of sage. Issue opened on upstream ​linbox/issues/138.

[chapoton]

what is the status of this ticket ?

[vdelecroix]

Replying to chapoton:

what is the status of this ticket ?

Waiting for upstream compilation fix, see 28.

[cpernet]

As mentionned in ​linbox/issues/138, passing a long instead of a size_t to BlasVector? constructor is causing an weird type resolution problem, which I can still not fully understand. Meanwhile 2 options can fix the problem:

• have sage cython code send size_t parameters to LinBox, as required by the interface. This means doing explicit casts (inelegant) or changing the type of ncols and nrows in matrix classes (heavy change)
• removing the templated constructor in LinBox::BlasVector that is in the way. This change is not harmfull for the library, but it feels not right to remove features for some weird type resolution issue on OSX due to sage not respecting the interface. This second option would also delay the fixing of this issue to the integration of the next release of LinBox or patching sage's LinBox.

[vdelecroix]

Replying to cpernet:

As mentionned in ​linbox/issues/138, passing a long instead of a size_t to BlasVector? constructor is causing an weird type resolution problem, which I can still not fully understand. Meanwhile 2 options can fix the problem:

• have sage cython code send size_t parameters to LinBox, as required by the interface. This means doing explicit casts (inelegant) or changing the type of ncols and nrows in matrix classes (heavy change)

One possible obstruction for the change Py_ssize_t -> size_t is that matrices can be accessed with negative indices. Perhaps it is not so dramatic. However, if this change happend we want to do it coherently with vectors, etc.

For the other proposition, calling properly LinBox with casting where necessary do not look too ugly for me. We can arrange so that it only appears in the code at sage/libs/linbox/.

• removing the templated constructor in LinBox::BlasVector that is in the way. This change is not harmfull for the library, but it feels not right to remove features for some weird type resolution issue on OSX due to sage not respecting the interface. This second option would also delay the fixing of this issue to the integration of the next release of LinBox or patching sage's LinBox.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​a9377c3	23214: use linbox in sparse integer matrices
​4ec1fa5	23214: disable sage extension from linbox

[vdelecroix]

Please Mac OSX testing needed!

[vklein]

Replying to vdelecroix:

Please Mac OSX testing needed!

make ptestlong is running.

[vklein]

We have almost the same error but for matrix_modn_sparse.

[sagelib-8.5.beta1] In file included from build/cythonized/sage/matrix/matrix_modn_sparse.cpp:725:
[sagelib-8.5.beta1] In file included from /Users/doctorant/sage/local/include/linbox/matrix/dense-matrix.h:79:
[sagelib-8.5.beta1] In file included from /Users/doctorant/sage/local/include/linbox/matrix/densematrix/blas-matrix.h:44:
[sagelib-8.5.beta1] /Users/doctorant/sage/local/include/linbox/vector/blas-vector.h:306:11: error: member reference base type 'const long' is not a structure or union
[sagelib-8.5.beta1]                         _size(V.size()),_1stride(1),_rep(V.size(), F.zero),_ptr(&_rep[0]),_field(&F)
[sagelib-8.5.beta1]                               ~^~~~~
[sagelib-8.5.beta1] build/cythonized/sage/matrix/matrix_modn_sparse.cpp:12191:19: note: in instantiation of function template specialization 'LinBox::BlasVector<Givaro::ZRing<Givaro::Integer>, std::__1::vector<Givaro::Integer, std::__1::allocator<Givaro::Integer> > >::BlasVector<long>' requested here
[sagelib-8.5.beta1]   __pyx_v_b = new LinBox::DenseVector<Givaro::ZRing<Givaro::Integer>>(__pyx_v_givZZ, __pyx_v_self->__pyx_base.__pyx_base.__pyx_base.__pyx_base.__pyx_base._nrows);
[sagelib-8.5.beta1]                   ^
[sagelib-8.5.beta1] In file included from build/cythonized/sage/matrix/matrix_modn_sparse.cpp:725:
[sagelib-8.5.beta1] In file included from /Users/doctorant/sage/local/include/linbox/matrix/dense-matrix.h:79:
[sagelib-8.5.beta1] In file included from /Users/doctorant/sage/local/include/linbox/matrix/densematrix/blas-matrix.h:44:
[sagelib-8.5.beta1] /Users/doctorant/sage/local/include/linbox/vector/blas-vector.h:306:38: error: member reference base type 'const long' is not a structure or union
[sagelib-8.5.beta1]                         _size(V.size()),_1stride(1),_rep(V.size(), F.zero),_ptr(&_rep[0]),_field(&F)
[sagelib-8.5.beta1]                                                          ~^~~~~
[sagelib-8.5.beta1] /Users/doctorant/sage/local/include/linbox/vector/blas-vector.h:306:61: warning: field '_rep' is uninitialized when used here [-Wuninitialized]
[sagelib-8.5.beta1]                         _size(V.size()),_1stride(1),_rep(V.size(), F.zero),_ptr(&_rep[0]),_field(&F)
[sagelib-8.5.beta1]                                                                                  ^
[sagelib-8.5.beta1] /Users/doctorant/sage/local/include/linbox/vector/blas-vector.h:129:13: error: type 'long' cannot be used prior to '::' because it has no members
[sagelib-8.5.beta1]                         typename OtherVector::const_iterator jt = V.begin();
[sagelib-8.5.beta1]                                  ^
[sagelib-8.5.beta1] /Users/doctorant/sage/local/include/linbox/vector/blas-vector.h:311:4: note: in instantiation of function template specialization 'LinBox::BlasVector<Givaro::ZRing<Givaro::Integer>, std::__1::vector<Givaro::Integer, std::__1::allocator<Givaro::Integer> > >::createBlasVector<long>' requested here
[sagelib-8.5.beta1]                         createBlasVector(V);
[sagelib-8.5.beta1]                         ^
[sagelib-8.5.beta1] build/cythonized/sage/matrix/matrix_modn_sparse.cpp:12191:19: note: in instantiation of function template specialization 'LinBox::BlasVector<Givaro::ZRing<Givaro::Integer>, std::__1::vector<Givaro::Integer, std::__1::allocator<Givaro::Integer> > >::BlasVector<long>' requested here
[sagelib-8.5.beta1]   __pyx_v_b = new LinBox::DenseVector<Givaro::ZRing<Givaro::Integer>>(__pyx_v_givZZ, __pyx_v_self->__pyx_base.__pyx_base.__pyx_base.__pyx_base.__pyx_base._nrows);
[sagelib-8.5.beta1]                   ^
[sagelib-8.5.beta1] 8 warnings and 3 errors generated.
[sagelib-8.5.beta1] error: command 'gcc' failed with exit status 1
[sagelib-8.5.beta1] [ 3/86] make[4]: *** [sage] Error 1

[cpernet]

There was indeed some casts missing. I pushed a fixed and checked that all casts have been done. @vklein, could you check again please?

---

New commits:

​d0f5a99

fix missing cast to size_t

[vklein]

Sure.

[chapoton]

I got 3 failing doctests, see patchbot.

[vklein]

Tests are still running on OSX. That mean off course that all compilation problems are solved.

[vklein]

I get 5 errors on OSX:

sage -t --long src/sage/combinat/designs/ext_rep.py  # Killed due to abort
sage -t --long src/sage/doctest/external.py  # Killed due to abort
sage -t --long src/sage/matrix/matrix_modn_sparse.pyx  # 3 doctests failed

matrix_modn_sparse.pyx errors are the sames as the patchbot ones and the two others are not ticket dependent (see Justin C. Walker message in ​Sage 8.5.beta1 released ).

[vdelecroix]

My bad. Error should be fixed now.

[chapoton]

forgot to push ?

[vdelecroix]

Not quite.

---

New commits:

​ea33422

23214: fix doctests

[tscrim]

There seems to be a test failure on one of the patchbots:

File "src/sage/matrix/matrix_space.py", line 1982, in sage.matrix.matrix_space.MatrixSpace._an_element_
Failed example:
    M.density()
Expected:
    99/1000000
Got:
    49/500000
**********************************************************************
1 item had failures:
   1 of   7 in sage.matrix.matrix_space.MatrixSpace._an_element_
    [383 tests, 1 failure, 6.40 s]
----------------------------------------------------------------------
sage -t --long --warn-long 54.6 src/sage/matrix/matrix_space.py  # 1 doctest failed

Can you reproduce this locally (sorry, I cannot test this right now on my laptop)?

[vdelecroix]

This concerns petitbonum and, on the other hand, everything is fine on pc78-math.

I can not reproduce either.

[vdelecroix]

This one on ​tmonteil-debian-stretch-32 (32 bits) is more annoying

[sagelib-8.5.beta2] build/cythonized/sage/matrix/matrix_modn_sparse.cpp:8931:132: error: no matching function for call to 'LinBox::GaussDomain<Givaro::Modular<long long unsigned int, long long unsigned int> >::InPlaceLinearPivoting(size_t&, uint64_t&, LinBox::SparseMatrix<Givaro::Modular<long long unsigned int, long long unsigned int>, LinBox::SparseMatrixFormat::SparseSeq>&, size_t, size_t)'
[sagelib-8.5.beta2]    (void)(__pyx_v_dom->InPlaceLinearPivoting(__pyx_v_A_rank, __pyx_v_A_det, (__pyx_v_A[0]), __pyx_v_A->rowdim(), __pyx_v_A->coldim()));
[sagelib-8.5.beta2]                                                                                                                                     ^

@Clément: any clue of what might be wrong on 32 bit architecture?

[vdelecroix]

@Clement: note that in LinBox the function is declared as (linbox/algorithms/gauss.h line 227)

// Sparsest method
//   erases elements while computing rank/det.
template <class _Matrix>                                      
unsigned long& InPlaceLinearPivoting(unsigned long &rank,
                                     Element& determinant,
                                     _Matrix        &A,
                                     unsigned long Ni,
                                     unsigned long Nj) const;

Here I don't understand why rank, Ni and Nj are not size_t as they should rather be.

[vdelecroix]

I filled the ​Issue 144 on linbox website.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​98ac984

23214: follow LinBox declarations

[vdelecroix]

Good news: tests pass on 32 bits.

@vinklein: can you test one more time on OSX?

If everything went fine on OSX, I think we should get this ticket merged as is. I will change unsigned long to size_t as soon as LinBox gets updated.

[tscrim]

This is also good on my machine and on a 32 and 64 bit patchbot.

[vklein]

Replying to vdelecroix:

@vinklein: can you test one more time on OSX?

Sure.

[vklein]

Build & tests pass on OSX (98ac984)

[vdelecroix]

Merci Vincent!

[vklein]

De rien !

[dimpase]

I strongly suspect that this ticket caused #28586 (breakage in solving sparse linear systems over prime fields), as it was the most substantial related change that got in in Sage 8.5.beta3. (and I bisected #28586 to 8.5.beta3).