use LinBox as native matrix representation for dense matrices over GF(p)

[malb]

Copying to and from LinBox uses up precious RAM and the point of fast linear algebra is to deal with large matrices. We should consider switching to LinBox as the native representation of matrices over GF(p)

Without Patch

sage: A = random_matrix(GF(97),2000,2000)
sage: %time A*A
CPU times: user 9.66 s, sys: 0.12 s, total: 9.77 s
Wall time: 9.82 s

With Patch

sage: A = random_matrix(GF(97),2000,2000)
sage: %time A*A
CPU times: user 1.32 s, sys: 0.00 s, total: 1.32 s
Wall time: 1.35 s

Magma

> A:=RandomMatrix(GF(97),2000,2000);
> time C:=A*A;                      
Time: 1.560

• Install http://sage.math.washington.edu/home/malb/spkgs/linbox-1.1.6.p5.spkg
• Apply attachment: trac_4260-linbox_default.patch
• Apply attachment: trac_4260-dense_ctypes_template.patch
• Apply attachment: trac_4260-matrix-modn-docs.patch
• Apply attachment: trac_4260_more_doctests.patch
• Apply attachment: trac_4260_echelonformdomain.patch
• Apply attachment: trac_4260-minor_fixes.patch
• Apply attachment: trac_4260_bugfix.patch

CC: @simon-king-jena @rbeezer @sagetrac-drkirkby

Component: linear algebra

Keywords: linbox, sd32, sd34

Author: Burcin Erocal, Martin Albrecht, Rob Beezer

Reviewer: Burcin Erocal, Simon King, Martin Albrecht, Jeroen Demeyer

Merged: sage-4.8.alpha3

Issue created by migration from https://trac.sagemath.org/ticket/4260

[ClementPernet]

comment:1

I will work on it as a coding sprint at SD10.

[burcin]

Author: Burcin Erocal

[burcin]

comment:3

I finally rebased the patch from SD16. The template class in the patch contains the updates made to the modn_dense class since then, like changes to the sig_* functions. Apparently the modn_dense class representation now allows permuting the rows by changing pointers in the _matrix array. We can't allow that if we want to pass the _entries to linbox, so I skipped those changes.

Sage builds with the attached patches, and you can construct matrices. However, there are lots of bugs, some linbox wrappers are still stubs, etc. Expect crashes and wrong results.

With the patch applied, I get a crash with the following:

sage: a = matrix(GF(97),3,4,range(12))
sage: a.echelonize()
*** glibc detected *** python: free(): invalid next size (fast): 0x000000000270b370 ***
======= Backtrace: =========
<snip>

AFAICT, the new cython code is an exact copy of the wrapper function in linbox-sage.C. Here is what valgrind says:

==3026== Invalid write of size 8
==3026==    at 0x39E49EF1: T.4552 (ffpack_ludivine.inl:420)
==3026==    by 0x39E49AA0: T.4552 (ffpack_ludivine.inl:486)
==3026==    by 0x39E4ABBF: __pyx_pf_4sage_6matrix_24matrix_modn_dense_double_26M
atrix_modn_dense_template_20_echelonize_linbox(_object*, _object*) (ffpack.h:113
2)
==3026==    by 0x4E74082: PyObject_Call (abstract.c:2492)
==3026==    by 0x39E2CA8A: __pyx_pf_4sage_6matrix_24matrix_modn_dense_double_26M
atrix_modn_dense_template_19echelonize(_object*, _object*, _object*) (matrix_mod
n_dense_double.cpp:8738)
==3026==    by 0x4F17FF9: PyEval_EvalFrameEx (ceval.c:3706)
==3026==    by 0x4F19CDC: PyEval_EvalCodeEx (ceval.c:2968)
==3026==    by 0x4F19DB1: PyEval_EvalCode (ceval.c:522)
==3026==    by 0x4F19083: PyEval_EvalFrameEx (ceval.c:4401)
==3026==    by 0x4F19CDC: PyEval_EvalCodeEx (ceval.c:2968)
<snip lots more Py* lines>
==3026==  Address 0x6ca16e8 is 0 bytes after a block of size 24 alloc'd
==3026==    at 0x4C267CE: malloc (vg_replace_malloc.c:236)
==3026==    by 0x39E4AA5A: __pyx_pf_4sage_6matrix_24matrix_modn_dense_double_26Matrix_modn_dense_template_20_echelonize_linbox(_object*, _object*) (memory.h:32)
==3026==    by 0x4E74082: PyObject_Call (abstract.c:2492)
==3026==    by 0x39E2CA8A: __pyx_pf_4sage_6matrix_24matrix_modn_dense_double_26Matrix_modn_dense_template_19echelonize(_object*, _object*, _object*) (matrix_modn_dense_double.cpp:8738)
==3026==    by 0x4F17FF9: PyEval_EvalFrameEx (ceval.c:3706)
==3026==    by 0x4F19CDC: PyEval_EvalCodeEx (ceval.c:2968)
==3026==    by 0x4F19DB1: PyEval_EvalCode (ceval.c:522)
==3026==    by 0x4F19083: PyEval_EvalFrameEx (ceval.c:4401)
==3026==    by 0x4F19CDC: PyEval_EvalCodeEx (ceval.c:2968)
==3026==    by 0x4F18074: PyEval_EvalFrameEx (ceval.c:3802)
<snip lots of Py* lines>

I'd appreciate any pointers about the problem above, though I don't know if I'll have the time to come back to this before the bug days in August (when I presume Martin will take over?).

[rbeezer]

comment:4

These are the files in sage/matrix with failures:

        sage -t  devel/sage-main/sage/matrix/matrix_cyclo_dense.pyx # 22 doctests failed
        sage -t  devel/sage-main/sage/matrix/strassen.pyx # 2 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix0.pyx # 2 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix_integer_dense.pyx # 5 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix_space.py # 1 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix_window_modn_dense.pyx # 1 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix_modn_sparse.pyx # 1 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix_integer_dense_saturation.py # 0 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix_rational_dense.pyx # 44 doctests failed
        sage -t  devel/sage-main/sage/matrix/matrix2.pyx # Time out
        sage -t  devel/sage-main/sage/matrix/matrix_modn_dense.pyx # Time out
        sage -t  devel/sage-main/sage/matrix/matrix_modn_dense_template.pxi # Time out

[malb]

Attachment: trac_4260-linbox_default.patch.gz

make matrix space constructor use the new classes

[malb]

Description changed:

--- 
+++ 
@@ -1 +1,25 @@
 Copying to and from LinBox uses up precious RAM and the point of fast linear algebra is to deal with large matrices. We should consider switching to LinBox as the native representation of matrices over GF(p)
+
+**Without Patch**
+
+```
+!python
+sage: A = random_matrix(GF(97),2000,2000)
+sage: %time A*A
+CPU times: user 9.66 s, sys: 0.12 s, total: 9.77 s
+Wall time: 9.82 s
+```
+
+**With Patch**
+
+```
+!python
+sage: A = random_matrix(GF(97),2000,2000)
+sage: %time A*A
+CPU times: user 1.32 s, sys: 0.00 s, total: 1.32 s
+Wall time: 1.35 s
+```
+
+* **Install** http://sage.math.washington.edu/home/malb/spkgs/linbox-1.1.6.p4.spkg
+* **Apply** [attachment: trac_4260-linbox_default.patch](https://github.com/sagemath/sage-prod/files/10642185/trac_4260-linbox_default.patch.gz)
+* **Apply** [attachment: trac_4260-dense_ctypes_template.patch](https://github.com/sagemath/sage-prod/files/10642186/trac_4260-dense_ctypes_template.patch.gz)

[malb]

comment:5

I fixed a few issues and segfaults but the thing is far from done. However, one can probably do higher level stuff now, i.e. it shouldn't crash that much any more.

We need a new LinBox SPKG because Modular<float> didn't have a NonZeroRandIter which is needed by the charpoly code. LinBox 1.1.7 fixes this issue but I tried unsuccessfully to upgrade to 1.1.7 for like 10 hours (cf. #11718).