Row rank profile / row pivots: direct extraction from echelon form over prime fields

[vneiger]

Fixes #40572

• ensure the use of the LUdivine / LQUP strategy with the flag FfpackSlabRecursive
• extract row rank profile using FFPACK's RankProfileFromLU

Concerning the first point, this strategy was already the one used in the single threaded case, since it is the default method in the FFPACK implementation. However it was not the default in the multi threaded case, which means there was actually a bug in Sage when calling _echelonize_linbox(efd=False) with 2 or more threads:

sage: mat = matrix(GF(3), [[1,0,1,0],[1,0,0,0],[1,0,0,0],[0,1,0,0]])
sage: mat.echelon_form()
[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 0]
sage: Parallelism().set("linbox", nproc=2)
sage: mat = matrix(GF(3), [[1,0,1,0],[1,0,0,0],[1,0,0,0],[0,1,0,0]])
sage: mat.echelon_form()
[1 0 0 0]
[0 0 1 0]
[0 1 0 0]
[0 0 0 0]

(the latter does not meet Sage's requirements of the reduced row echelon form)

Some more complete tests concerning the row rank profile are to be done before merging.

[vneiger]

@Biffo89 : to add to the tests, could you please check if that solves your issue in PR #40508 ? thanks.

[github-actions]

Documentation preview for this PR (built with commit 0dfb949; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[vneiger]

The previous version was not working correctly due to the row rank profile not being preserved by the used Gaussian elimination. The latest commits switched to the PLUQ method for which FFPACK guarantees to preserve both rank profiles. The class also overrides pivots and pivot_rows to ensure that when either of them is called, both rank profiles are cached (to avoid an unnecessary call to echelon_form when the other rank profile is requested).

It seems that the echelonization methods could be simplified and could also easily provide access to the transformation, but this is another issue.

[Biffo89]

@Biffo89 : to add to the tests, could you please check if that solves your issue in PR #40508 ? thanks.

I've manually tested all m x n matrices over GF(3), with m, n <= 4. All calls to pivots and pivot_rows are consistent. The following code terminates successfully (although of course takes a long time to do so):

for i in range(1,5):
     for j in range(1,5):
         for m in MatrixSpace(GF(3),i,j):
             m._clear_cache()
             p1 = m.pivots()
             q1 = m.pivot_rows()
             m._clear_cache()
             q2 = m.pivot_rows()
             p2 = m.pivots()
             m._clear_cache()
             mt = m.transpose()
             mt._clear_cache()
             p3 = mt.pivot_rows()
             q3 = mt.pivots()
             mt._clear_cache()
             q4 = mt.pivots()
             p4 = mt.pivot_rows()
             if p1 != p2 or p2 != p3 or p3 != p4 or q1 != q2 or q2 != q3 or q3 != q4:
                 print(m)

[vneiger]

@Biffo89 : to add to the tests, could you please check if that solves your issue in PR #40508 ? thanks.

I've manually tested all m x n matrices over GF(3), with m, n <= 4. All calls to pivots and pivot_rows are consistent. The following code terminates successfully (although of course takes a long time to do so):

Thanks for the feedback! Apparently I cannot ask you to be a reviewer on github; I think you need to ask to join the sagemath github organization for that (or maybe the triage team). I also added @xcaruso as a reviewer in case he wishes to review the changes.

[xcaruso]

I do not know very much about linbox but the PR looks good to me, expect the two comments below.

[xcaruso]

Also, I underline that the documentation of the function linbox_echelonize is not really correct.
Maybe, you could fix it (and add doctests for linbox_echelonize_efd as well) in the PR.

[xcaruso]

Do you want to update the doctests or not?

[vneiger]

Do you want to update the doctests or not?

Doing it, should be ready in a few minutes.

[vneiger]

Do you want to update the doctests or not?

I added/fixed documentation for linbox_echelonize and linbox_echelonize_efd. For the doctests, note that the echelonize variants are tested against each other and against a naive version via the last doctest in echelonize (which tests on random matrices of size 100 x 100 over a few prime fields).

[xcaruso]

Looks very nice! Btw, what does efd mean?

[vneiger]

Looks very nice! Btw, what does efd mean?

Echelon Form Domain. This name is not used much in LinBox anymore, it probably was at the time where this was integrated. Also, the description of the two variants

          - ``linbox`` -- uses the LinBox library (wrapping fflas-ffpack)
          - ``linbox_noefd`` -- uses the FFPACK directly, less memory and faster (default)

and particularly the more detailed comment

        - ``efd`` -- if ``True`` LinBox's ``EchelonFormDomain``
          implementation is used, which is faster than the direct
          ``LinBox::FFPACK`` implementation, since the latter also
          computes the transformation matrix (which we
          ignore). However, ``efd=True`` uses more memory than FFLAS
          directly (default: ``True``)

seem contradictory (and probably outdated with respect to the current ReducedRowEchelonForm in FFLAS-FFPACK). I'm not sure having the two variants makes sense anymore, however we should definitely incorporate FFPACK's mechanism to efficiently recover the transformation (related to #37480 ). I keep this in mind for a future enhancement... I did not want to include too much in this PR to avoid breaking things that work.

[xcaruso]

Everything looks good to me and failures do not seem related to this PR. I give a positive review.