gh-40579: Row rank profile / row pivots: direct extraction from echelon form over prime fields
    
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.
    
URL: #40579
Reported by: Vincent Neiger
Reviewer(s): Vincent Neiger, Xavier Caruso

Author: Release Manager

src/sage/libs/linbox/fflas.pxd

@@ -29,6 +29,7 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFLAS":
         FflasTrans
 
     ctypedef enum FFLAS_SIDE:
+        FflasLeft
         FflasRight
 
     # double
@@ -80,6 +81,16 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFLAS":
              size_t C_stride, size_t numthreads)
 
 cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFPACK":
+    ctypedef enum FFPACK_LU_TAG:
+        FfpackTileRecursive
+
+    void RankProfileFromLU (size_t* P, size_t N, size_t R,
+                            size_t* rkprofile, FFPACK_LU_TAG LuTag)
+
+    void PLUQtoEchelonPermutation (size_t N, size_t R, size_t * P, size_t * outPerm)
+
+    void MathPerm2LAPACKPerm (size_t * LapackP, size_t * MathP, size_t N)
+
     # double
     bint IsSingular (Modular_double F,
                      size_t nrows, size_t ncols, Modular_double.Element* A,
@@ -104,11 +115,14 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFPACK":
 
     size_t ReducedRowEchelonForm (Modular_double F, size_t a, size_t b,
                                   Modular_double.Element* matrix,
-                                  size_t s, size_t* P, size_t* Q)
+                                  size_t s, size_t* P, size_t* Q,
+                                  bool transform, FFPACK_LU_TAG LuTag)
 
     size_t pReducedRowEchelonForm (Modular_double F, size_t a, size_t b,
                                   Modular_double.Element* matrix,
-                                  size_t s, size_t* P, size_t* Q, bool transform, size_t numthreads)
+                                  size_t s, size_t* P, size_t* Q,
+                                  bool transform, size_t numthreads,
+                                  FFPACK_LU_TAG LuTag)
 
     Modular_double.Element* Solve (Modular_double F, size_t M,
            Modular_double.Element* A, size_t lda,
@@ -158,11 +172,14 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFPACK":
 
     size_t ReducedRowEchelonForm (Modular_float F, size_t a, size_t b,
                                   Modular_float.Element* matrix,
-                                  size_t s, size_t* P, size_t* Q)
+                                  size_t s, size_t* P, size_t* Q,
+                                  bool transform, FFPACK_LU_TAG LuTag)
 
     size_t pReducedRowEchelonForm (Modular_float F, size_t a, size_t b,
                                   Modular_float.Element* matrix,
-                                  size_t s, size_t* P, size_t* Q, bool transform, size_t numthreads)
+                                  size_t s, size_t* P, size_t* Q,
+                                  bool transform, size_t numthreads,
+                                  FFPACK_LU_TAG LuTag)
 
     Modular_float.Element* Solve (Modular_float F, size_t M,
            Modular_float.Element* A, size_t lda,

src/sage/matrix/matrix2.pyx

@@ -8649,6 +8649,7 @@ cdef class Matrix(Matrix1):
         if v is not None:
             self.cache('echelon_transformation', v)
         self.cache('pivots', E.pivots())
+
         if transformation and v is not None:
             return (E, v)
         else: