Opened 5 years ago

Closed 5 years ago

Last modified 5 years ago

#16803 closed enhancement (fixed)

Reimplement matrix_integer_dense using FLINT

Reported by: mmasdeu Owned by:
Priority: major Milestone: sage-6.4
Component: linear algebra Keywords: flint, matrix
Cc: pbruin Merged in:
Authors: Marc Masdeu Reviewers: William Stein, Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: 812a509 (Commits) Commit:
Dependencies: Stopgaps:

Description

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.

Change History (98)

comment:1 Changed 5 years ago by mmasdeu

  • Branch set to u/mmasdeu/16803-matrix_integer_dense_flint
  • Commit set to a51a2fb4f4576685c13f9ec01bd252b04e6dd503

comment:2 Changed 5 years ago by mmasdeu

By the way, after rebasing it does not compile. I will fix it tomorrow hopefully.

comment:3 Changed 5 years ago by fredrik.johansson

A word of caution about algorithm defaults: FLINT should generally be faster for small input but IML/Linbox etc. probably have much better code for things like nullspace, charpoly... as soon as the matrices start to get large.

Replacing _rational_kernel_iml with _rational_kernel_flint, for example, is strange. Surely both algorithms should be available.

comment:4 Changed 5 years ago by git

  • Commit changed from a51a2fb4f4576685c13f9ec01bd252b04e6dd503 to d180ee24dc9547212762499ff271b5a119dadfa9

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

f1862b3Trac 16803: Fixed compilation problems.
d180ee2Trac 16803: Finished re-adding legacy functionality.

comment:5 Changed 5 years ago by mmasdeu

I rewrote the code to have all the previous functionality available. The code compiles and passes all the doctests in the matrix/ folder. I will have it run all the doctests during the weekend.

comment:6 Changed 5 years ago by git

  • Commit changed from d180ee24dc9547212762499ff271b5a119dadfa9 to e26c3472968016bc9dcb82cca70803bbd946b886

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

e26c347Trac 16803: Modified to pass doctests from outside matrix folder.

comment:7 Changed 5 years ago by git

  • Commit changed from e26c3472968016bc9dcb82cca70803bbd946b886 to f2ac03c5052e46df91afa80420ed2b08d43d2823

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

a6e9a3dMoved ntl dependency inside the pyx file.
f2ac03cTrac 16803: Added sig_on/off, changed some defaults after benchmarking.

comment:8 Changed 5 years ago by git

  • Commit changed from f2ac03c5052e46df91afa80420ed2b08d43d2823 to 166bae537c52aa4eee968275dcd2f2aef9f52a82

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

888451bTrac 16803: Fixed a missing sig_off() and doctests.
166bae5All doctests pass, documentation builds well also.

comment:9 follow-up: Changed 5 years ago by was

Hi,

A first remark playing around with timings is that in sage-6.3 somebody or something has definitely *completely massively screwed* up Linbox, or how Linbox is built, or how we use linbox, which was the default representation and library for arithmetic. Consider this:

sage: n=200; set_random_seed(0); a = random_matrix(ZZ,n); b = random_matrix(ZZ,n)
sage: %time c=a._multiply_linbox(b)
CPU times: user 326 ms, sys: 3.63 ms, total: 330 ms
Wall time: 347 ms
sage: %time c=a._multiply_multi_modular(b)
CPU times: user 51.5 ms, sys: 1.6 ms, total: 53.1 ms
Wall time: 56.1 ms
sage: n=400; set_random_seed(0); a = random_matrix(ZZ,n); b = random_matrix(ZZ,n)
sage: %time c=a._multiply_linbox(b)
CPU times: user 5.28 s, sys: 22.8 ms, total: 5.3 s
Wall time: 5.37 s
sage: %time c=a._multiply_multi_modular(b)
CPU times: user 230 ms, sys: 6.08 ms, total: 236 ms
Wall time: 248 ms

Linbox should easily beat the multiply_multi_modular approach -- that's just basically a naive version of the same thing linbox would be doing if it were properly built and linked. Instead, Linbox is massively slower. In fact, almost exactly as much slower as I would expect if the wrong BLAS were being linked. Note that I tested against the new flint backend and now a._multiply_multi_modular(b) is asymptotically similar to flint, though asymptotically (for n=2000 already), they are pretty close to the same in speed (8.55s versus 11.4s). In any case, somebody should really look into why linbox is so screwed up -- we switched to it initially since it was much better than doing a._multiply_multi_modular(b) and also our own strassen implementation.

comment:10 Changed 5 years ago by was

OK, I've been skeptically looking at and testing this patch much of today. I made a few tiny changes, but frankly I think this is overall AWESOME. FLINT's linear algebra seems quite good.

The following three methods are now still gone. Add them back so I can run some randomized consistency tests more easily...

     a._multiply_classical      a._multiply_linbox         a._multiply_multi_modular

When these are added back, and I run tests of correctness by comparing with them, then I'll give this an enthusiastic positive review.

comment:11 Changed 5 years ago by was

  • Status changed from new to needs_review

comment:12 Changed 5 years ago by was

  • Status changed from needs_review to needs_work

comment:13 Changed 5 years ago by git

  • Commit changed from 166bae537c52aa4eee968275dcd2f2aef9f52a82 to 5f21e3dbfda7a126d4ba388681c75f7c559a1e24

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

7f4f903Merge branch 'develop' into matflint
5f21e3dMinor changes during review. Put back old multiplication functionality.

comment:14 Changed 5 years ago by git

  • Commit changed from 5f21e3dbfda7a126d4ba388681c75f7c559a1e24 to a13b79c9decdba47551a70f2ce1822a4e7e082bd

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

a13b79cFixed errors, all doctests pass now.

comment:15 Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

I have added the three functions back, and run the long tests in sage/matrix/.

comment:16 Changed 5 years ago by was

MAJOR PROBLEM! Now that multiply_multi_modular is back, I wrote code to do some random testing of correctness. I did not find situations where things are wrong... however, I quickly ran into major trouble because this ticket introduces a huge memory leak. I just defined two matrices and multiplied them a few hundred times, and had major ram problems. In particular, in sage-6.3 (without this branch):

sage: get_memory_usage()
1047.1875
sage: for i in range(50):
....:         a=random_matrix(ZZ,300); del a
....: 
sage: get_memory_usage()
1047.1875

With this ticket:

sage: get_memory_usage()
2132.09765625
sage: for i in range(50):
    a=random_matrix(ZZ,300); del a
....: 
sage: get_memory_usage()
2231.3203125

Conclusion: making nontrivial use of the code on this ticket will very quickly use up all RAM on your computer.

So... needs work until that is fixed.

comment:17 Changed 5 years ago by git

  • Commit changed from a13b79c9decdba47551a70f2ce1822a4e7e082bd to 0d37b701fbd845ee36fdaa6dcf62bae71bf4d6af

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

0d37b70Fixed memory leak.

comment:18 Changed 5 years ago by was

I don't think the changeset with name "Fixed memory leak" fixes the memory leak, at least not for matrix multiplication. I'm not sure if you're done yet. With your fix, if you run the functions test1, test2 or test3 in the script 2014-08-29-180550-stupid-test.sage below, then the memory usage (print out second) grows quite rapidly and doesn't go down. Doing the same thing in sage-6.3 results in almost no leakage of the heap. In particular, with your patch after running those 3 tests I ended up with over 2GB extra memory used (oer startup), even after doing import gc; gc.collect(). Doing the identical thing with sage-6.3 uses only a few MB's.

def t0(n):
    a = random_matrix(ZZ,n)
    b = random_matrix(ZZ,n)
    assert a*b == a*b

def t(n):
    a = random_matrix(ZZ,n)
    b = random_matrix(ZZ,n)
    assert a*b == a._multiply_multi_modular(b)

def test1():
    for n in range(1,50):
        print n, get_memory_usage()
        set_random_seed(0)
        for i in range(100):
            t(n)

def test2():
    for n in range(51,150):
        print n, get_memory_usage()
        set_random_seed(0)
        for i in range(10):
            t(n)


def t2(n):
    a = random_matrix(ZZ,n,2*n)
    b = random_matrix(ZZ,2*n,n)
    at = a.transpose()
    bt = b.transpose()
    assert a*b == a._multiply_multi_modular(b)
    assert bt*at == bt._multiply_multi_modular(at)

def test3():
    for n in range(1,50):
        print n, get_memory_usage()
        set_random_seed(0)
        for i in range(50):
            t2(n)

comment:19 Changed 5 years ago by fredrik.johansson

Just from a quick glance at the code, it looks like _lift_crt doesn't mpz_clear out tmp.

comment:20 Changed 5 years ago by git

  • Commit changed from 0d37b701fbd845ee36fdaa6dcf62bae71bf4d6af to f6ac760964f32d4e1105a541333e4df38a22aecd

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

2f2fa14Merge branch 'develop' into matflint
f6ac760Fixed other memory leaks. Now all doctests pass again as well.

comment:21 follow-up: Changed 5 years ago by mmasdeu

Thanks Fredrik, I fixed that (although it was not the only problem, by far).

The previous commit was not finished yet (sorry for the confusing commit message). I believe that now the memory leaks are gone (at least those that were detected by William's tests). The main problem was that the get_unsafe_mpz function was initializing the value that it modifies in place, and all functions calling it were also mpz_init'ing it as well. All should be good now.

comment:22 in reply to: ↑ 9 ; follow-up: Changed 5 years ago by jdemeyer

Replying to was:

A first remark playing around with timings is that in sage-6.3 somebody or something has definitely *completely massively screwed* up Linbox, or how Linbox is built, or how we use linbox, which was the default representation and library for arithmetic.

Are you sure this isn't an imaginary slow-down? Just for fun, I installed some old versions of Sage on boxen and executed the command

$ time ./sage -c '_ = Matrix(ZZ,1000,1000,range(10^6))^5'

I tested Sage versions 2.9(*), 2.11(*), 3.0.6, 3.2.3, 4.1, 4.5, 4.7, 5.0, 5.12, 6.1.1, 6.3 and the above command always took about 12 minutes. Note that boxen is not an idle machine, so it's difficult to get accurate timings, but there is certainly no massive slow-down.

(*) For these versions, ATLAS didn't compile so I used the ATLAS package from Sage 3.0.6.

comment:23 in reply to: ↑ 22 Changed 5 years ago by was

Replying to jdemeyer:

Replying to was:

A first remark playing around with timings is that in sage-6.3 somebody or something has definitely *completely massively screwed* up Linbox, or how Linbox is built, or how we use linbox, which was the default representation and library for arithmetic.

Are you sure this isn't an imaginary slow-down? Just for fun, I installed some old versions of Sage on boxen and executed the command

$ time ./sage -c '_ = Matrix(ZZ,1000,1000,range(10^6))^5'

I tested Sage versions 2.9(*), 2.11(*), 3.0.6, 3.2.3, 4.1, 4.5, 4.7, 5.0, 5.12, 6.1.1, 6.3 and the above command always took about 12 minutes. Note that boxen is not an idle machine, so it's difficult to get accurate timings, but there is certainly no massive slow-down.

(*) For these versions, ATLAS didn't compile so I used the ATLAS package from Sage 3.0.6.

Thanks for testing. I wouldn't be at all surprised if the problem was introduced quite a long time ago -- I haven't benchmarked sage linear algebra at all since around the time when Clement was my postdoc (about 2007-2008), and at that time definitely a properly built linbox provided the fastest matrix multiply (and many other things). It might still today -- I don't know.

comment:24 in reply to: ↑ 21 Changed 5 years ago by was

Replying to mmasdeu:

Thanks Fredrik, I fixed that (although it was not the only problem, by far).

The previous commit was not finished yet (sorry for the confusing commit message). I believe that now the memory leaks are gone (at least those that were detected by William's tests). The main problem was that the get_unsafe_mpz function was initializing the value that it modifies in place, and all functions calling it were also mpz_init'ing it as well. All should be good now.

Marc -- can you update the code in our SD61 project in sagetest-marc so I can test this out. I tried what was there, and it leaked like a sieve.

comment:25 Changed 5 years ago by mmasdeu

It just finished compiling now. I got tired of bug-fixing in the cloud (some things are better done in private), and worked from my local copy in the last bit.

comment:26 Changed 5 years ago by was

  • Reviewers set to William Stein
  • Status changed from needs_review to positive_review

OK, I'm now officially thrilled/very happy by this code :-)

comment:27 Changed 5 years ago by vbraun

Merge conflict, probably from #16583

comment:28 Changed 5 years ago by vbraun

  • Status changed from positive_review to needs_work

comment:29 Changed 5 years ago by jdemeyer

  • Branch changed from u/mmasdeu/16803-matrix_integer_dense_flint to u/jdemeyer/ticket/16803
  • Created changed from 08/12/14 16:39:14 to 08/12/14 16:39:14
  • Modified changed from 09/05/14 11:03:05 to 09/05/14 11:03:05

comment:30 Changed 5 years ago by jdemeyer

  • Commit changed from f6ac760964f32d4e1105a541333e4df38a22aecd to 2e66cee05777a1565f8aa012730df777f798cb0f

Merged with

git merge ticket/16583 -X theirs

Compiles fine, not tested again.


New commits:

3bda21ccleanup of fmpz_t and fmpz_poly_t to avoid Cython warnings about referencing before assignment
dfc8ce1silence Cython warning about complicated declaration
8bb6b48merge
14e730eMerge with 6.4.beta1
18a4842trac #16583, oops, correcting wrong merge
c416feeRemove stuff which is not needed
b3d6aa2Further cleanup
20c767cFix include of fmpq_poly.pxi in Cython doctest
2e66ceeMerge branch 'ticket/16583' into ticket/16803

comment:31 Changed 5 years ago by jdemeyer

  • Status changed from needs_work to positive_review

comment:32 Changed 5 years ago by mmasdeu

  • Branch changed from u/jdemeyer/ticket/16803 to u/mmasdeu/ticket/16803
  • Commit changed from 2e66cee05777a1565f8aa012730df777f798cb0f to 22e6896aadf0d6f99aeebc29041d6fd247569738
  • Status changed from positive_review to needs_review

I have just realized that (in a bad design choice) I was initializing allocating _rows and _entries for every single matrix created, which is not very efficient. Now it is only done when needed (i.e. when calling _init_linbox()). Changes are minimal, and I ran the tests (long version) in /matrix/ for safety, but I will resubmit for review.


New commits:

22e6896Trac 16803: Initialize _rows and _entries only when needed (to save memory).

comment:33 Changed 5 years ago by was

Be more careful about freeing. Right now, in all cases, you do

sage_free(self._rows)
sage_free(self._entries)

You're depending on self._rows being automatically initialized to NULL (?), and sage_free doing the right thing when passed NULL. Maybe that works, but I would be more comfortable if it was explicit, e.g., put self._rows=NULL and only free self._rows if not NULL (or check the source code and see whether it already works that way).

comment:34 Changed 5 years ago by git

  • Commit changed from 22e6896aadf0d6f99aeebc29041d6fd247569738 to 70afcd6d8a4092515da48e2e33d7024a80c8ec8e

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

70afcd6Trac 16803: free _rows and _entries only when needed.

comment:35 Changed 5 years ago by mmasdeu

Done.


New commits:

70afcd6Trac 16803: free _rows and _entries only when needed.

comment:36 follow-up: Changed 5 years ago by vbraun

Its safe to call sage_free with NULL, but an unintialized pointer can have a non-null value. You should explicitly set all pointers to NULL in __cinit__.

Last edited 5 years ago by vbraun (previous) (diff)

comment:37 in reply to: ↑ 36 Changed 5 years ago by was

  • Status changed from needs_review to positive_review

Replying to vbraun:

Its safe to call sage_free with NULL, but an unintialized pointer can have a non-null value. You should explicitly set all pointers to NULL in __cinit__.

Thanks for clarifying this. The change he made does explicitly set pointers to NULL now, but it also checks for NULL before calling sage_free. He should probably change that part back.

Anyway, positive review.

comment:38 Changed 5 years ago by git

  • Commit changed from 70afcd6d8a4092515da48e2e33d7024a80c8ec8e to 5dd0b79e71c51367c1602cb70e0d92291e12e931
  • Status changed from positive_review to needs_review

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

5dd0b79Trac 16803: No need to check for uninitialized pointers, so reverting

comment:39 Changed 5 years ago by mmasdeu

  • Status changed from needs_review to positive_review

I cleaned that up.

comment:40 Changed 5 years ago by vbraun

The commit is ok but the commit message is wrong ;-) You initialize it to NULL, thats the whole point.

comment:41 Changed 5 years ago by mmasdeu

Argh! By uninitialized I meant unallocated. Sorry...

comment:42 Changed 5 years ago by vbraun

  • Status changed from positive_review to needs_work

Conflicts with #15946, please fix.

comment:43 Changed 5 years ago by jdemeyer

  • Branch changed from u/mmasdeu/ticket/16803 to u/jdemeyer/ticket/16803
  • Modified changed from 09/06/14 11:04:18 to 09/06/14 11:04:18

comment:44 Changed 5 years ago by jdemeyer

  • Commit changed from 5dd0b79e71c51367c1602cb70e0d92291e12e931 to 6d819238634f4ca885827df0ce217344899d415b

I merged with #15946, but that revealed some genuine bugs with this ticket, hang on...


New commits:

6fc03e3Get rid of mpz_t_offset hacks
4a5d820Merge branch 'ticket/16910' into HEAD
d76f8e7Fix mpz_t, mpq_t, mpf_t declarations
a52cff2Use mpz_srcptr for padics pow_computer
6d81923Merge branch 'ticket/15946' into ticket/16803

comment:45 Changed 5 years ago by git

  • Commit changed from 6d819238634f4ca885827df0ce217344899d415b to 0c75c2ba9475a2eb02807387375b88140d1137c3

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

0c75c2bRemove unused _lift_crt_rr_with_lcm() function

comment:46 Changed 5 years ago by jdemeyer

Still needs work because of

sage -t src/sage/combinat/q_analogues.py
**********************************************************************
File "src/sage/combinat/q_analogues.py", line 200, in sage.combinat.q_analogues.q_binomial
Failed example:
    q_binomial(4,2,matrix([[2,1],[-1,3]]))
Exception raised:
    Traceback (most recent call last):
      File "/usr/local/src/sage-git/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 480, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/usr/local/src/sage-git/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 843, in compile_and_execute
        exec compiled in globs
      File "<doctest sage.combinat.q_analogues.q_binomial[14]>", line 1, in <module>
        q_binomial(Integer(4),Integer(2),matrix([[Integer(2),Integer(1)],[-Integer(1),Integer(3)]]))
      File "/usr/local/src/sage-git/local/lib/python2.7/site-packages/sage/combinat/q_analogues.py", line 308, in q_binomial
        return numerat/denomin
      File "element.pyx", line 2743, in sage.structure.element.Matrix.__div__ (build/cythonized/sage/structure/element.c:20137)
      File "matrix_integer_dense.pyx", line 923, in sage.matrix.matrix_integer_dense.Matrix_integer_dense._matrix_times_matrix_ (build/cythonized/sage/matrix/matrix_integer_dense.c:11570)
      File "c_lib.pyx", line 97, in sage.ext.c_lib.sig_raise_exception (build/cythonized/sage/ext/c_lib.c:1208)
    SignalError: Segmentation fault
**********************************************************************

Given that this test also fails (not with a segfault though) in the patchbot, I guess the merging isn't to blame.

Please fix this...

comment:47 Changed 5 years ago by mmasdeu

  • Branch changed from u/jdemeyer/ticket/16803 to u/mmasdeu/16803-2
  • Commit changed from 0c75c2ba9475a2eb02807387375b88140d1137c3 to fcb36b1e57fc645c0372617711b2a9ec0423b77f
  • Status changed from needs_work to needs_review

I have fixed this by changing the code of sage.structure.element.Matrix.div, which was assuming that the _matrix_times_matrix method of matrix.matrix_integer_dense would know how to deal with matrices over the rationals. This is a low-level private function which doesn't do type checking, which I think is the right thing to do.

What the current patch does is, when writing something like A/B for matrices A and B, to take the denominator out of ~B, multiply and then put it back. I have run make ptest and all passed except for some that have nothing to do with matrices, for example this one:

sage -t src/sage/rings/integer.pyx
**********************************************************************
File "src/sage/rings/integer.pyx", line 1954, in sage.rings.integer.Integer.__pow__
Failed example:
    2^(2^63-2)
Expected:
    Traceback (most recent call last):
    ...
    MemoryError: failed to allocate 1152921504606847008 bytes
Got:
    Fatal Python error: Insufficient memory
    Traceback (most recent call last):
      File "/usr/local/sage/sage-6.3.beta6/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 480, in _run
        self.execute(example, compiled, test.globs)
      File "/usr/local/sage/sage-6.3.beta6/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 839, in execute
        exec compiled in globs
      File "<doctest sage.rings.integer.Integer.__pow__[27]>", line 1, in <module>
        Integer(2)**(Integer(2)**Integer(63)-Integer(2))
      File "integer.pyx", line 1999, in sage.rings.integer.Integer.__pow__ (build/cythonized/sage/rings/integer.c:14710)
      File "c_lib.pyx", line 85, in sage.ext.c_lib.sig_raise_exception (build/cythonized/sage/ext/c_lib.c:1040)
    RuntimeError: Aborted
**********************************************************************

comment:48 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

The patchbot reports a segmentation fault in matrix_integer_dense.pyx, which is likely due to this ticket.

comment:49 Changed 5 years ago by mmasdeu

I cannot reproduce it. And I don't know where the error comes from, other than the coercion model not checking for types before calling the low-level functions...

~/sagetest-marc$ ./sage -t --long src/sage/modular/modform/ambient_R.py
Running doctests with ID 2014-09-09-08-48-58-b66cd1ea.
Doctesting 1 file.
sage -t --long --warn-long 26.1 src/sage/modular/modform/ambient_R.py
    [26 tests, 29.09 s]
----------------------------------------------------------------------
All tests passed!
----------------------------------------------------------------------
Total time for all tests: 29.2 seconds
    cpu time: 27.7 seconds
    cumulative wall time: 29.1 seconds

comment:50 Changed 5 years ago by jdemeyer

Your last fix is almost certainly going to break over other rings. In any case, this change requires extra doctests over other rings.

I would advise the following: the matrix division coercion model bug is independent of this ticket. I recommend you to create a new ticket for that and make this ticket depend on that.

comment:51 follow-up: Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

I propose the following, which I think is a good solution (especially since the plan is to move matrices over rationals also to FLINT).

def __div__(left,right):
    return left * (~right)

I have tried this and it passes the tests. I will push the code in a minute.

comment:52 Changed 5 years ago by mmasdeu

  • Branch changed from u/mmasdeu/16803-2 to u/mmasdeu/16803-matflint-3
  • Commit changed from fcb36b1e57fc645c0372617711b2a9ec0423b77f to 4cf706d0974ffc39af15e4b24b6e77b8bc984d53

Last 10 new commits:

7069aecFurther cleanup
f1756fbFix include of fmpq_poly.pxi in Cython doctest
78d9bfcTrac 16803: Initialize _rows and _entries only when needed (to save memory).
42243b2Trac 16803: free _rows and _entries only when needed.
902fcc5Trac 16803: No need to check for uninitialized pointers, so reverting
1df2b4dGet rid of mpz_t_offset hacks
dc9ef8cFix mpz_t, mpq_t, mpf_t declarations
c2d5546Use mpz_srcptr for padics pow_computer
d1635d4Remove unused _lift_crt_rr_with_lcm() function
4cf706dTrac 16803: fixed __div__ function.

comment:53 in reply to: ↑ 51 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

Replying to mmasdeu:

I propose the following, which I think is a good solution (especially since the plan is to move matrices over rationals also to FLINT).

def __div__(left,right):
    return left * (~right)

I have tried this and it passes the tests.

You need to add more doctests to this __div__ function.

comment:54 Changed 5 years ago by git

  • Commit changed from 4cf706d0974ffc39af15e4b24b6e77b8bc984d53 to 59d2a7b11c45338c59b5ed16c09167c2dec6809d

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

59d2a7bTrac 16803: added doctests to __div__ function in structure/element.pyx

comment:55 Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

Done.

comment:56 Changed 5 years ago by jdemeyer

  • Reviewers changed from William Stein to William Stein, Jeroen Demeyer
  • Status changed from needs_review to needs_work

I really think more doctests are needed, also over rings which are not ZZ or QQ (see also comment 50). Perhaps GF(3) and RDF and ZZ['t'] to have some wider coverage?

And you should use the :trac:`16803` syntax to refer to Trac tickets (but given that this ticket has nothing to do with division of matrices, perhaps it makes more sense to just remove that comment).

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:57 Changed 5 years ago by jdemeyer

Also, doctest division of a matrix by a scalar...

comment:58 Changed 5 years ago by git

  • Commit changed from 59d2a7b11c45338c59b5ed16c09167c2dec6809d to 6095b103e1bf20752a58926a63a492bbf81cb004

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

6095b10Added more tests.

comment:59 Changed 5 years ago by git

  • Commit changed from 6095b103e1bf20752a58926a63a492bbf81cb004 to f9974ba608c96861627c3f01cc58c5030536e9bd

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

f9974baAnd some more.

comment:60 Changed 5 years ago by mmasdeu

Done!

comment:61 Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

comment:62 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

In the last patches, you messed up the indentation from 4 to 3 spaces.

comment:63 Changed 5 years ago by git

  • Commit changed from f9974ba608c96861627c3f01cc58c5030536e9bd to 17bb2b784da9b2cecc40e8593fb32047ed6ae9cb

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

17bb2b7Fixed wrong indentation.

comment:64 Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

Sorry. Fixed that.

comment:65 Changed 5 years ago by jdemeyer

  • Status changed from needs_review to needs_work

I think this still leaks memory.

The following command works without this patch, but fails with this patch:

( ulimit -v 2000000; ./sage -t --long src/sage/combinat/root_system/branching_rules.py )

comment:66 Changed 5 years ago by git

  • Commit changed from 17bb2b784da9b2cecc40e8593fb32047ed6ae9cb to cd815217e562f84118f53458b18cfd9b68c975a2

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

83bec69Fixed typo in polynomial.
cd81521Fixed bug introduced in outside package.

comment:67 Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

The bug reported in comment:65 may be related to an issue with infinite recursion. I traced it down to the modification of the div function. With the current patch the command in comment:65 works, and all the tests in sage/structure and sage/matrix pass. I am now running make ptest and will report if something arises.

comment:68 Changed 5 years ago by mmasdeu

  • Branch changed from u/mmasdeu/16803-matflint-3 to u/mmasdeu/16803-matflint-4
  • Commit changed from cd815217e562f84118f53458b18cfd9b68c975a2 to 6359d72733219492dccba27caf10a04b909eda97

I have rebased it again. It looks like some other code is a moving target and it is hard to keep pace...


Last 10 new commits:

2b87370Trac 16803: No need to check for uninitialized pointers, so reverting
50cb2aeFix mpz_t, mpq_t, mpf_t declarations
e52c1ffRemove unused _lift_crt_rr_with_lcm() function
fa5a4afTrac 16803: fixed __div__ function.
1b86766Trac 16803: added doctests to __div__ function in structure/element.pyx
3c6f808Added more tests.
d87e348And some more.
5b8ee2dFixed wrong indentation.
eab8a1eFixed bug introduced in outside package.
6359d72Fixed compilation errors.

comment:69 Changed 5 years ago by was

  • Status changed from needs_review to positive_review

Since test pass, and the issues Jereon found were addressed, I'm switching this back to positive review.

comment:70 Changed 5 years ago by vbraun

  • Status changed from positive_review to needs_work

This still fails, I'm pretty sure just or your branch as well:

sage -t --long --warn-long 70.7 src/sage/structure/element.pyx
**********************************************************************
File "src/sage/structure/element.pyx", line 2770, in sage.structure.element.Matrix.__div__
Failed example:
    a / b
Expected:
    [-0.159090909091 -0.295454545455]
    [  2.86363636364 -0.681818181818]
Got:
    [-0.15909090909090906 -0.29545454545454547]
    [   2.863636363636364  -0.6818181818181817]
**********************************************************************
1 item had failures:
   1 of  22 in sage.structure.element.Matrix.__div__
    [455 tests, 1 failure, 38.78 s]

comment:71 Changed 5 years ago by git

  • Commit changed from 6359d72733219492dccba27caf10a04b909eda97 to 6757e11b36f7bc647eae1adcf636bed76fb3e85d

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

6757e11Trac 16803: Fix a doctest for RDF matrices.

comment:72 Changed 5 years ago by mmasdeu

If I am not mistaken, it is the current output which is correct (I checked with Magma also). I had myself added this doctest last week, sorry!

comment:73 Changed 5 years ago by vbraun

  • Status changed from needs_work to positive_review

comment:74 Changed 5 years ago by vbraun

Numerical noise on 32-bit, should be marked with a suitable # rel tol

sage -t --long src/sage/structure/element.pyx
**********************************************************************
File "src/sage/structure/element.pyx", line 2770, in sage.structure.element.Matrix.__div__
Failed example:
    a / b
Expected:
    [-0.15909090909090906 -0.29545454545454547]
    [   2.863636363636364  -0.6818181818181817]
Got:
    [ -0.1590909090909091 -0.29545454545454547]
    [   2.863636363636364  -0.6818181818181815]
**********************************************************************
1 item had failures:

comment:75 Changed 5 years ago by vbraun

  • Status changed from positive_review to needs_work

comment:76 Changed 5 years ago by git

  • Commit changed from 6757e11b36f7bc647eae1adcf636bed76fb3e85d to 812a5099731dd5e819ebce0b19ee744de439a6ea

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

812a509Trac 16803: Added rel tol to doctest.

comment:77 Changed 5 years ago by mmasdeu

Done.

comment:78 Changed 5 years ago by mmasdeu

  • Status changed from needs_work to needs_review

comment:79 Changed 5 years ago by vbraun

  • Status changed from needs_review to positive_review

Linear algebra should really give the right answer to a few ulp, and not just 1e-10. But the test is about testing that the division operator works, so ok.

comment:80 Changed 5 years ago by vbraun

  • Branch changed from u/mmasdeu/16803-matflint-4 to 812a5099731dd5e819ebce0b19ee744de439a6ea
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:81 Changed 5 years ago by jdemeyer

  • Commit 812a5099731dd5e819ebce0b19ee744de439a6ea deleted

Follow-ups: #17080 and #17090.

comment:82 Changed 5 years ago by jdemeyer

Did any of the reviewers (Volker?) actually read the code on this ticket? I've started doing this and there are many things to be improved... perhaps this ticket was merged too fast.

comment:83 Changed 5 years ago by vbraun

Is it perfect? No. But better than before? IMHO yes. Faster? For sure.

comment:84 follow-up: Changed 5 years ago by was

Replying to jdemeyer:

Did any of the reviewers (Volker?) actually read the code on this ticket? I've started doing this and there are many things to be improved...

I certainly read the version of the code 5 weeks ago, and had many, many, many suggestions for improvement, found numerous major bugs, etc.

perhaps this ticket was merged too fast.

One has to balance that with the fact that for whatever reason, without this ticket a vast range of linear algebra in Sage is so utterly orders of magnitude too slow, as to make it painful and useless for research work. This new code is dramatically faster than before.

comment:85 in reply to: ↑ 84 ; follow-up: Changed 5 years ago by jdemeyer

Replying to was:

This new code is dramatically faster than before.

I only hope we're not trading correctness for speed.

comment:86 in reply to: ↑ 85 ; follow-ups: Changed 5 years ago by was

Replying to jdemeyer:

Replying to was:

This new code is dramatically faster than before.

I only hope we're not trading correctness for speed.

Well when I looked over the code 5 weeks ago there are numerous memory leaks and other causes for concern (e.g., unitialized memory assumed to be 0). But those were I think all addressed. Did you find anything that is actually incorrect? Or does the overall code just make you unhappy and nervous? Or are you not happy building on FLINT?

Last edited 5 years ago by was (previous) (diff)

comment:87 Changed 5 years ago by was

I just have to add jdemeyer that I'm really glad you're looking at this too, since you have an amazing eye for quality (similar to your recent remarks about the pexpect interfaces).

comment:88 in reply to: ↑ 86 Changed 5 years ago by jdemeyer

Replying to was:

Did you find anything that is actually incorrect?

The segmentation fault in _lmul_, which is the reason for creating #17090.

Or does the overall code just make you unhappy and nervous?

When looking at this code, there are some sloppy things (like why is the _det_pari() function duplicated?). Overall, I have less faith in sloppy code than clear code.

Or are you not happy building on FLINT?

That's not the problem. In general, I think it's good to rely on external packages if they do a good job (which seems to be the case here).

comment:89 in reply to: ↑ 86 Changed 5 years ago by jdemeyer

Last edited 5 years ago by jdemeyer (previous) (diff)

comment:90 Changed 5 years ago by jdemeyer

What's with the changes in padics in this ticket? Surely, yet another mistake...

comment:91 Changed 5 years ago by jdemeyer

Yet another follow-up: #17094

comment:92 follow-up: Changed 5 years ago by mmasdeu

I must say that IMO some of the comments (82 and 90) made above are not very constructive. jdemeyer questions the thoroughness of the reviewing process, while clearly he himself didn't read the comments at the beginning of the ticket (at least those that refer to the convenience of leaving in the legacy functions).

After the experience with this ticket I feel less inclined to take on the analogous job for matrix_rational_dense. Sorry!

comment:93 in reply to: ↑ 92 ; follow-up: Changed 5 years ago by jdemeyer

Replying to mmasdeu:

I must say that IMO some of the comments (82 and 90) made above are not very constructive.

82 was made out of frustration by seeing so many things here which are dubious. So okay, that's not constructive.

However in 90, I noticed a mistake, mentioned it in the comment and fixed it in #17090. How is that not constructive?

while clearly he himself didn't read the comments at the beginning of the ticket (at least those that refer to the convenience of leaving in the legacy functions).

That's true, but after William Stein's comment, I have fixed that back in #17090.

I feel less inclined to take on the analogous job for matrix_rational_dense.

Please don't. I think it's absolutely good to improve matrices like you did here. The fact that I am making so many comments and creating follow-up tickets means that I care. I make mistakes, you make mistakes, that's normal and that's why have this review process.

My only complaint (and I still stand behind this) is that this ticket was given positive_review too soon. Ideally, those 3 follow-up tickets should have been reviewer patches on this ticket.

comment:94 in reply to: ↑ 93 ; follow-ups: Changed 5 years ago by mmasdeu

Replying to jdemeyer:

82 was made out of frustration by seeing so many things here which are dubious. So okay, that's not constructive.

I can understand your frustration. But the bar is high enough (python + cython + sage oddities + legacy code + git + trac, and I am probably missing something here) that makes it very hard for people to contribute, and these comments do not help. But I take your last comment in good faith :-).

Also, if there are other things that are dubious please share them because I (and possibly others!) would like to learn from it.

However in 90, I noticed a mistake, mentioned it in the comment and fixed it in #17090. How is that not constructive?

Surely I did not go and touch p-adics for no reason. These changes weren't there in the first version of the my patch, and I had to add them because during the time it took to from the first version until this ticket got closed other changes were being merged, and these made the patch to stop compiling. Changing those lines in p-adics is how I managed to fix it, but the hackish idea was not mine, I saw it used elsewhere. I appreciate you having fixed it back.

comment:95 in reply to: ↑ 94 Changed 5 years ago by jdemeyer

Replying to mmasdeu:

Surely I did not go and touch p-adics for no reason. These changes weren't there in the first version of the my patch, and I had to add them because during the time it took to from the first version until this ticket got closed other changes were being merged, and these made the patch to stop compiling.

I think you simply did something compiling Sage, because those changes to p-adics are totally unrelated to matrices. I cannot imagine how the change to matrices could break building p-adics or vice versa. And indeed, I simply undid the changes in #17090 without any bad consequences.

comment:96 in reply to: ↑ 94 Changed 5 years ago by was

Replying to mmasdeu:

Replying to jdemeyer:

82 was made out of frustration by seeing so many things here which are dubious. So okay, that's not constructive.

I can understand your frustration. But the bar is high enough (python + cython + sage oddities + legacy code + git + trac, and I am probably missing something here) that makes it very hard for people to contribute, and these comments do not help. But I take your last comment in good faith :-).

Marc -- don't you want to be a bad ass computer programmer, in addition to mathematician? Just redo matrix_rational_dense, and you'll be way closer :-)

comment:97 in reply to: ↑ 94 ; follow-up: Changed 5 years ago by jdemeyer

Replying to mmasdeu:

Also, if there are other things that are dubious please share them because I (and possibly others!) would like to learn from it.

Well, you can look at the follow-up tickets #17090 and #17094 and see what I changed (and preferably review those tickets).

Some general things which haven't been mentioned yet:

  1. clean-up your code when you change things: when you remove code, try to remove everything related to the old code. When you change code, the old structure might not be applicable anymore, so refactor if needed (e.g. the sepatate branch for n <= 4 in determinant() no longer serves any purpose)
  2. also fix documentation: when you change something, also change the corresponding documentation. When you copy/paste something, be sure to update documentation and doctests (see the various set_unsafe methods).
  3. don't add unneeded import statements or includes, those add needless extra dependencies. Like, why did you add include "sage/libs/ntl/decl.pxi" in several places when you don't use any NTL code?
  4. try to use correct formatting: spacing, indentation, docstring formatting (this was mostly good, but a bit messed up in src/module_list.py). If anything, this makes your code look more professional and less sloppy.
  5. don't use assert for checking user input: use assert for when you know that something is true (but want to check just in case), not when you need to check that something is true.

comment:98 in reply to: ↑ 97 Changed 5 years ago by mmasdeu

Replying to jdemeyer:

Thanks! These comments are helpful indeed.

William -- I rather try to be a better mathematician...it has the advantage that although you get constantly evaluated, the evaluations are mostly anonymous :-). But that was in my to-do list. And it still is, despite my previous comments in this ticket!

Note: See TracTickets for help on using tickets.