From a2b0a601aa28a67fc48ddeb11720ab2c15a4c33a Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Tue, 12 Aug 2025 09:50:31 +0200 Subject: [PATCH 1/8] extract Row rank profile from LU + Ensure FfpackSlabRecursive strategy in parallel version --- src/sage/libs/linbox/fflas.pxd | 20 ++++++++++++---- .../matrix/matrix_modn_dense_template.pxi | 24 ++++++++++++------- 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/src/sage/libs/linbox/fflas.pxd b/src/sage/libs/linbox/fflas.pxd index e8e6cd6de40..4f8b43b76e9 100644 --- a/src/sage/libs/linbox/fflas.pxd +++ b/src/sage/libs/linbox/fflas.pxd @@ -80,6 +80,12 @@ 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: + FfpackSlabRecursive + + void RankProfileFromLU (size_t* P, size_t N, size_t R, + size_t* rkprofile, FFPACK_LU_TAG LuTag) + # double bint IsSingular (Modular_double F, size_t nrows, size_t ncols, Modular_double.Element* A, @@ -104,11 +110,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 +167,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, diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 8e9e7bf0f89..dfd7d06332e 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -94,7 +94,7 @@ from cysignals.signals cimport sig_check, sig_on, sig_off from sage.libs.gmp.mpz cimport * from sage.libs.linbox.fflas cimport FFLAS_TRANSPOSE, FflasNoTrans, FflasTrans, \ - FflasRight, vector, list as std_list + FfpackSlabRecursive, FflasRight, vector, list as std_list, RankProfileFromLU from libcpp cimport bool from sage.parallel.parallelism import Parallelism @@ -176,12 +176,13 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr Return the reduced row echelon form of this matrix. """ if linbox_is_zero(modulus, entries, nrows, ncols): - return 0, [] + return 0, [], [] cdef Py_ssize_t i, j cdef ModField *F = new ModField(modulus) cdef size_t* P = check_allocarray(nrows, sizeof(size_t)) cdef size_t* Q = check_allocarray(ncols, sizeof(size_t)) + cdef size_t* rkprofile = check_allocarray(nrows, sizeof(size_t)) cdef Py_ssize_t r cdef size_t nbthreads @@ -190,9 +191,11 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr if nrows * ncols > 1000: sig_on() if nbthreads > 1 : - r = pReducedRowEchelonForm(F[0], nrows, ncols, entries, ncols, P, Q, transform, nbthreads) + r = pReducedRowEchelonForm(F[0], nrows, ncols, entries, + ncols, P, Q, transform, nbthreads, FfpackSlabRecursive) else : - r = ReducedRowEchelonForm(F[0], nrows, ncols, entries, ncols, P, Q) + r = ReducedRowEchelonForm(F[0], nrows, ncols, entries, ncols, P, Q, + transform, FfpackSlabRecursive) if nrows * ncols > 1000: sig_off() @@ -204,12 +207,16 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr applyP(F[0], FflasRight, FflasNoTrans, nrows, 0, r, entries, ncols, Q) + RankProfileFromLU(Q, nrows, r, rkprofile, FfpackSlabRecursive) + cdef list pivots = [int(Q[i]) for i in range(r)] + cdef list rrp = [int(rkprofile[i]) for i in range(r)] sig_free(P) sig_free(Q) + sig_free(rkprofile) del F - return r, pivots + return r, pivots, rrp cdef inline linbox_echelonize_efd(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols): # See trac #13878: This is to avoid sending invalid data to linbox, @@ -1595,7 +1602,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense): - ``gauss`` -- uses a custom slower `O(n^3)` Gauss elimination implemented in Sage - - ``all`` -- compute using both algorithms and verify that + - ``all`` -- compute using all algorithms and verify that the results are the same - ``**kwds`` -- these are all ignored @@ -1788,12 +1795,13 @@ cdef class Matrix_modn_dense_template(Matrix_dense): r, pivots = linbox_echelonize_efd(self.p, self._entries, self._nrows, self._ncols) else: - r, pivots = linbox_echelonize(self.p, self._entries, - self._nrows, self._ncols) + r, pivots, rrp = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols) verbose('done with echelonize', t) self.cache('in_echelon_form', True) self.cache('rank', r) self.cache('pivots', tuple(pivots)) + if not efd: + self.cache('pivot_rows', tuple(rrp)) def _echelon_in_place_classical(self): """ From 49dd4f5c40ecffe8e1d1b00e297be26ca5629239 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Thu, 14 Aug 2025 11:57:25 +0200 Subject: [PATCH 2/8] rely on TileRecursive for rank profile preservation; cache both pivots and pivot_rows when either of them is asked for --- src/sage/libs/linbox/fflas.pxd | 8 + src/sage/matrix/matrix2.pyx | 1 + .../matrix/matrix_modn_dense_template.pxi | 167 +++++++++++++++--- 3 files changed, 155 insertions(+), 21 deletions(-) diff --git a/src/sage/libs/linbox/fflas.pxd b/src/sage/libs/linbox/fflas.pxd index 4f8b43b76e9..ce8d35e465d 100644 --- a/src/sage/libs/linbox/fflas.pxd +++ b/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 @@ -81,11 +82,18 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFLAS": cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFPACK": ctypedef enum FFPACK_LU_TAG: + FfpackTileRecursive FfpackSlabRecursive 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) + + void LAPACKPerm2MathPerm (size_t * MathP, size_t * LapackP, size_t N) + # double bint IsSingular (Modular_double F, size_t nrows, size_t ncols, Modular_double.Element* A, diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index e5f2ebb453f..a6fbac8b102 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/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: diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index dfd7d06332e..e0f64360631 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -94,7 +94,9 @@ from cysignals.signals cimport sig_check, sig_on, sig_off from sage.libs.gmp.mpz cimport * from sage.libs.linbox.fflas cimport FFLAS_TRANSPOSE, FflasNoTrans, FflasTrans, \ - FfpackSlabRecursive, FflasRight, vector, list as std_list, RankProfileFromLU + FfpackTileRecursive, FfpackSlabRecursive, FflasLeft, FflasRight, vector, list as std_list, \ + RankProfileFromLU, PLUQtoEchelonPermutation, MathPerm2LAPACKPerm, LAPACKPerm2MathPerm + from libcpp cimport bool from sage.parallel.parallelism import Parallelism @@ -182,7 +184,13 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr cdef ModField *F = new ModField(modulus) cdef size_t* P = check_allocarray(nrows, sizeof(size_t)) cdef size_t* Q = check_allocarray(ncols, sizeof(size_t)) - cdef size_t* rkprofile = check_allocarray(nrows, sizeof(size_t)) + + cdef size_t* Qmath = check_allocarray(ncols, sizeof(size_t)) + cdef size_t* Qtmp = check_allocarray(ncols, sizeof(size_t)) + cdef size_t* Qperm = check_allocarray(ncols, sizeof(size_t)) + + cdef size_t* rrp = check_allocarray(nrows, sizeof(size_t)) + cdef size_t* crp = check_allocarray(ncols, sizeof(size_t)) cdef Py_ssize_t r cdef size_t nbthreads @@ -192,31 +200,61 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr sig_on() if nbthreads > 1 : r = pReducedRowEchelonForm(F[0], nrows, ncols, entries, - ncols, P, Q, transform, nbthreads, FfpackSlabRecursive) + ncols, P, Q, transform, nbthreads, FfpackTileRecursive) else : r = ReducedRowEchelonForm(F[0], nrows, ncols, entries, ncols, P, Q, - transform, FfpackSlabRecursive) + transform, FfpackTileRecursive) if nrows * ncols > 1000: sig_off() + # extract row and column rank profiles + RankProfileFromLU(P, nrows, r, rrp, FfpackTileRecursive) + RankProfileFromLU(Q, ncols, r, crp, FfpackTileRecursive) + + cdef list pivots = [int(crp[i]) for i in range(r)] + cdef list pivot_rows = [int(rrp[i]) for i in range(r)] + + # row permutation into echelon + PLUQtoEchelonPermutation(ncols, r, Q, Qperm) + applyP(F[0], FflasLeft, FflasNoTrans, ncols, 0, r, entries, ncols, Qperm) + + # fill top left with identity for i in range(nrows): for j in range(r): - (entries+i*ncols+j)[0] = 0 + (entries + i*ncols + j)[0] = 0 if ientries, ncols, Q) + (entries + i*ncols + i)[0] = 1 + + # column permutation to put pivots in correct columns + piv = 0 + ipiv = r + idx = 0 + while idx < r and ipiv < ncols: + if crp[idx] > piv: + crp[ipiv] = piv + ipiv += 1 + piv += 1 + else: + idx += 1 + piv += 1 + while ipiv < ncols: + crp[ipiv] = piv + ipiv += 1 + piv += 1 - RankProfileFromLU(Q, nrows, r, rkprofile, FfpackSlabRecursive) + MathPerm2LAPACKPerm(Q, crp, ncols) - cdef list pivots = [int(Q[i]) for i in range(r)] - cdef list rrp = [int(rkprofile[i]) for i in range(r)] + applyP(F[0], FflasRight, FflasNoTrans, nrows, 0, ncols, entries, ncols, Q) sig_free(P) sig_free(Q) - sig_free(rkprofile) + sig_free(rrp) + sig_free(crp) + sig_free(Qmath) + sig_free(Qtmp) + sig_free(Qperm) del F - return r, pivots, rrp + return r, pivots, pivot_rows cdef inline linbox_echelonize_efd(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols): # See trac #13878: This is to avoid sending invalid data to linbox, @@ -1607,14 +1645,20 @@ cdef class Matrix_modn_dense_template(Matrix_dense): - ``**kwds`` -- these are all ignored - OUTPUT: ``self`` is put in reduced row echelon form + OUTPUT: if ``self`` is known to be echelonized (the information is + cached), nothing is done; else, ``self`` is put in reduced row echelon + form and + + - the fact that ``self`` is now in echelon form is cached so future + calls to echelonize return immediately - the rank of ``self`` is computed and cached - - the pivot columns of ``self`` are computed and cached + - the pivot columns (a.k.a. column rank profile) of ``self`` are + computed and cached - - the fact that ``self`` is now in echelon form is recorded and - cached so future calls to echelonize return immediately + - only if using algorithm ``linbox_noefd``: return the pivot rows of + ``self`` before echelonization (a.k.a. its row rank profile) EXAMPLES:: @@ -1632,7 +1676,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense): ....: A = random_matrix(GF(13), 10, 10) sage: MS = parent(A) sage: B = A.augment(MS(1)) - sage: B.echelonize() + sage: rrp = B.echelonize() sage: A.rank() 10 sage: C = B.submatrix(0,10,10,10) @@ -1724,12 +1768,16 @@ cdef class Matrix_modn_dense_template(Matrix_dense): [0 0 0 0 0 0 0 0 0 0] sage: A = matrix(GF(97),3,4,range(12)) - sage: A.echelonize(); A + sage: A.echelonize() + (0, 1) + sage: A [ 1 0 96 95] [ 0 1 2 3] [ 0 0 0 0] sage: A.pivots() (0, 1) + sage: A.pivot_rows() + (0, 1) sage: for p in (3,17,97,127,1048573): ....: for i in range(10): @@ -1746,7 +1794,8 @@ cdef class Matrix_modn_dense_template(Matrix_dense): if algorithm == 'linbox': self._echelonize_linbox(efd=True) elif algorithm == 'linbox_noefd': - self._echelonize_linbox(efd=False) + rrp = self._echelonize_linbox(efd=False) + return rrp elif algorithm == 'gauss': self._echelon_in_place_classical() @@ -1777,6 +1826,9 @@ cdef class Matrix_modn_dense_template(Matrix_dense): ignore). However, ``efd=True`` uses more memory than FFLAS directly (default: ``True``) + OUTPUT: if ``efd`` is ``False``, return the row rank profile of + ``self`` + EXAMPLES:: sage: A = random_matrix(GF(7), 10, 20) @@ -1800,8 +1852,10 @@ cdef class Matrix_modn_dense_template(Matrix_dense): self.cache('in_echelon_form', True) self.cache('rank', r) self.cache('pivots', tuple(pivots)) + self.cache('pivot_rows', tuple(range(r))) + if not efd: - self.cache('pivot_rows', tuple(rrp)) + return tuple(rrp) def _echelon_in_place_classical(self): """ @@ -1857,8 +1911,79 @@ cdef class Matrix_modn_dense_template(Matrix_dense): start_row = start_row + 1 break self.cache('pivots', tuple(pivots)) + self.cache('pivot_rows', tuple(range(r))) self.cache('in_echelon_form', True) + def pivots(self): + """ + Return the pivot column positions of this matrix. + + OUTPUT: a tuple of Python integers: the position of the + first nonzero entry in each row of the echelon form. + + This returns a tuple so it is immutable; see :issue:`10752`. + + EXAMPLES:: + + sage: A = matrix(GF(7), 2, 2, range(4)) + sage: A.pivots() + (0, 1) + """ + v = self.fetch('pivots') + if v is not None: + return tuple(v) + + E = self.__copy__() + rrp = E.echelonize(algorithm="linbox_noefd") + E.set_immutable() + v = E.pivots() + self.cache('echelon_form', E) + self.cache('rank', E._cache['rank']) + self.cache('pivots', tuple(v)) + self.cache('pivot_rows', rrp) + return tuple(v) + + def pivot_rows(self): + """ + Return the pivot row positions for this matrix, which form the topmost + subset of the rows that span the row space and are linearly + independent. Also known as row rank profile. + + If this has already been computed and cached, this returns the cached + value. Otherwise, this computes an echelon form (using algorithm + ``linbox_noefd``), and deduces the pivot row positions to be returned. + In the latter case, this also caches other attributes at the same time + (see :meth:`Matrix_modn_dense_template.echelonize`). + + OUTPUT: a tuple of integers + + EXAMPLES:: + + sage: A = matrix(GF(3), [[1,0,1,0],[1,0,0,0],[1,0,0,0],[0,1,0,0]]) + sage: B = A.transpose() + sage: A + [1 0 1 0] + [1 0 0 0] + [1 0 0 0] + [0 1 0 0] + sage: A.pivot_rows() + (0, 1, 3) + sage: B.pivots() == A.pivot_rows() + True + """ + v = self.fetch('pivot_rows') + if v is not None: + return tuple(v) + + E = self.__copy__() + v = E.echelonize(algorithm="linbox_noefd") + E.set_immutable() + self.cache('echelon_form', E) + self.cache('rank', E._cache['rank']) + self.cache('pivots', E._cache['pivots']) + self.cache('pivot_rows', v) + return v + def right_kernel_matrix(self, algorithm='linbox', basis='echelon'): r""" Return a matrix whose rows form a basis for the right kernel From ae3f61417f10baa171037e083567fab5be2d9f53 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Thu, 14 Aug 2025 16:36:34 +0200 Subject: [PATCH 3/8] fix echelonize vs. special _echelonize_linbox: do not return value --- .../matrix/matrix_modn_dense_template.pxi | 22 +++++++++++-------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index e0f64360631..0425398f800 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -1657,8 +1657,9 @@ cdef class Matrix_modn_dense_template(Matrix_dense): - the pivot columns (a.k.a. column rank profile) of ``self`` are computed and cached - - only if using algorithm ``linbox_noefd``: return the pivot rows of - ``self`` before echelonization (a.k.a. its row rank profile) + - only if using algorithm ``linbox_noefd``: the pivot rows of ``self`` + before echelonization (a.k.a. its row rank profile) are computed and + cached EXAMPLES:: @@ -1768,9 +1769,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense): [0 0 0 0 0 0 0 0 0 0] sage: A = matrix(GF(97),3,4,range(12)) - sage: A.echelonize() - (0, 1) - sage: A + sage: A.echelonize(); A [ 1 0 96 95] [ 0 1 2 3] [ 0 0 0 0] @@ -1794,8 +1793,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense): if algorithm == 'linbox': self._echelonize_linbox(efd=True) elif algorithm == 'linbox_noefd': - rrp = self._echelonize_linbox(efd=False) - return rrp + self._echelonize_linbox(efd=False) elif algorithm == 'gauss': self._echelon_in_place_classical() @@ -1929,12 +1927,15 @@ cdef class Matrix_modn_dense_template(Matrix_dense): sage: A.pivots() (0, 1) """ + if not self.base_ring().is_field(): + raise NotImplementedError("Echelon form not implemented over '%s'." % self.base_ring()) + v = self.fetch('pivots') if v is not None: return tuple(v) E = self.__copy__() - rrp = E.echelonize(algorithm="linbox_noefd") + rrp = E._echelonize_linbox(efd=False) E.set_immutable() v = E.pivots() self.cache('echelon_form', E) @@ -1971,12 +1972,15 @@ cdef class Matrix_modn_dense_template(Matrix_dense): sage: B.pivots() == A.pivot_rows() True """ + if not self.base_ring().is_field(): + raise NotImplementedError("Echelon form not implemented over '%s'." % self.base_ring()) + v = self.fetch('pivot_rows') if v is not None: return tuple(v) E = self.__copy__() - v = E.echelonize(algorithm="linbox_noefd") + v = E._echelonize_linbox(efd=False) E.set_immutable() self.cache('echelon_form', E) self.cache('rank', E._cache['rank']) From c8a927f3881a1abd6778138ad1ef25dd06401c90 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 15 Aug 2025 17:11:14 +0200 Subject: [PATCH 4/8] improve documentation and add tests/examples --- .../matrix/matrix_modn_dense_template.pxi | 94 ++++++++++++++++--- 1 file changed, 80 insertions(+), 14 deletions(-) diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 0425398f800..ccf23fe3bb4 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -1914,18 +1914,59 @@ cdef class Matrix_modn_dense_template(Matrix_dense): def pivots(self): """ - Return the pivot column positions of this matrix. + Return the column pivot positions for this matrix, which form the + leftmost subset of the columns that span the column space and are + linearly independent. This coincides with the position of the first + nonzero entry in each row of the reduced row echelon form of ``self``, + and is also known as the column rank profile of ``self``. The returned + tuple is ordered increasingly. - OUTPUT: a tuple of Python integers: the position of the - first nonzero entry in each row of the echelon form. + If this has already been computed and cached, this returns the cached + value. Otherwise, this computes an echelon form (using algorithm + ``linbox_noefd``), and deduces the pivot positions to be cached and + returned. In the latter case, this also caches other attributes at the + same time: the reduced row echelon form of ``self`` and the rank of + ``self``, as well as its row pivot positions (also known as row rank + profile). + + OUTPUT: a tuple of `r` integers where `r` is the rank of ``self`` - This returns a tuple so it is immutable; see :issue:`10752`. + .. SEEALSO:: + + The method :meth:`Matrix_modn_dense_template.pivot_rows` computes + the row pivot positions, also known as row rank profile. EXAMPLES:: sage: A = matrix(GF(7), 2, 2, range(4)) sage: A.pivots() (0, 1) + + sage: A = matrix(GF(3), [[1,1,1,0],[0,0,0,1],[1,0,0,0]]) + sage: A + [1 1 1 0] + [0 0 0 1] + [1 0 0 0] + sage: A.pivots() == (0, 1, 3) + True + + sage: A = matrix(GF(65537), 5, 3, + ....: [ 223, 669, 21130, + ....: 13996, 41988, 21387, + ....: 39034, 51565, 40500, + ....: 14660, 43980, 3899, + ....: 12016, 36048, 9308]) + sage: A[:,1] == 3 * A[:,0] + True + sage: A.pivots() == (0, 2) + True + sage: A = matrix(GF(7), 3, 5, [2, 2, 4, 1, 4, + ....: 4, 4, 1, 4, 6, + ....: 5, 2, 5, 6, 6]) + sage: A[:,0] == 2 * A[:,1] + 3 * A[:,2] + True + sage: A.pivots() == (0, 1, 3) + True """ if not self.base_ring().is_field(): raise NotImplementedError("Echelon form not implemented over '%s'." % self.base_ring()) @@ -1946,30 +1987,55 @@ cdef class Matrix_modn_dense_template(Matrix_dense): def pivot_rows(self): """ - Return the pivot row positions for this matrix, which form the topmost + Return the row pivot positions for this matrix, which form the topmost subset of the rows that span the row space and are linearly - independent. Also known as row rank profile. + independent. This coincides with the position of the first nonzero + entry in each column of the reduced column echelon form of ``self``, + and is also known as the row rank profile of ``self``. The returned + tuple is ordered increasingly. If this has already been computed and cached, this returns the cached value. Otherwise, this computes an echelon form (using algorithm - ``linbox_noefd``), and deduces the pivot row positions to be returned. - In the latter case, this also caches other attributes at the same time - (see :meth:`Matrix_modn_dense_template.echelonize`). + ``linbox_noefd``), and deduces the pivot row positions to be cached and + returned. In the latter case, this also caches other attributes at the + same time: the reduced row echelon form of ``self`` and the rank of + ``self``, as well as its pivot indices (also known as column rank + profile). + + OUTPUT: a tuple of `r` integers where `r` is the rank of ``self`` - OUTPUT: a tuple of integers + .. SEEALSO:: + + The method :meth:`Matrix_modn_dense_template.pivots` computes + the column pivot positions, also known as column rank profile. EXAMPLES:: sage: A = matrix(GF(3), [[1,0,1,0],[1,0,0,0],[1,0,0,0],[0,1,0,0]]) - sage: B = A.transpose() sage: A [1 0 1 0] [1 0 0 0] [1 0 0 0] [0 1 0 0] - sage: A.pivot_rows() - (0, 1, 3) - sage: B.pivots() == A.pivot_rows() + sage: A.pivot_rows() == (0, 1, 3) + True + + sage: A = matrix(GF(65537), 3, 5, + ....: [ 223, 13996, 39034, 14660, 12016, + ....: 669, 41988, 51565, 43980, 36048, + ....: 21130, 21387, 40500, 3899, 9308]) + sage: A[1,:] == 3 * A[0,:] + True + sage: A.pivot_rows() == (0, 2) + True + sage: A = matrix(GF(7), 5, 3, [2, 4, 5, + ....: 2, 4, 2, + ....: 4, 1, 5, + ....: 1, 4, 6, + ....: 4, 6, 6]) + sage: A[0,:] == 2 * A[1,:] + 3 * A[2,:] + True + sage: A.pivot_rows() == (0, 1, 3) True """ if not self.base_ring().is_field(): From a1109fc8001d004acfcd68ae00beeea75a3af75b Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Tue, 19 Aug 2025 16:38:05 +0200 Subject: [PATCH 5/8] suggested fixes --- src/sage/matrix/matrix_modn_dense_template.pxi | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index ccf23fe3bb4..8935ed1e834 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -184,9 +184,6 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr cdef ModField *F = new ModField(modulus) cdef size_t* P = check_allocarray(nrows, sizeof(size_t)) cdef size_t* Q = check_allocarray(ncols, sizeof(size_t)) - - cdef size_t* Qmath = check_allocarray(ncols, sizeof(size_t)) - cdef size_t* Qtmp = check_allocarray(ncols, sizeof(size_t)) cdef size_t* Qperm = check_allocarray(ncols, sizeof(size_t)) cdef size_t* rrp = check_allocarray(nrows, sizeof(size_t)) @@ -248,11 +245,9 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr sig_free(P) sig_free(Q) + sig_free(Qperm) sig_free(rrp) sig_free(crp) - sig_free(Qmath) - sig_free(Qtmp) - sig_free(Qperm) del F return r, pivots, pivot_rows @@ -1677,7 +1672,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense): ....: A = random_matrix(GF(13), 10, 10) sage: MS = parent(A) sage: B = A.augment(MS(1)) - sage: rrp = B.echelonize() + sage: B.echelonize() sage: A.rank() 10 sage: C = B.submatrix(0,10,10,10) From cc9299436fabc6ecfc49b45d8eeab97e71176dfd Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Tue, 19 Aug 2025 16:44:43 +0200 Subject: [PATCH 6/8] documentation for linbox_echelonize --- src/sage/matrix/matrix_modn_dense_template.pxi | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 8935ed1e834..08222ac27e4 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -175,7 +175,12 @@ cdef inline bint linbox_is_zero(celement modulus, celement* entries, Py_ssize_t cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols): """ - Return the reduced row echelon form of this matrix. + In-place transform this matrix into its reduced row echelon form, and return + the rank `r` of this matrix as well as two lists of length `r`, sorted + increasingly. The first list gives the column rank profile of this matrix + (which is also that of its reduced row echelon form) while the second list + gives the row rank profile of the input matrix (which may differ from that + of its reduced row echelon form). """ if linbox_is_zero(modulus, entries, nrows, ncols): return 0, [], [] @@ -252,6 +257,12 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr return r, pivots, pivot_rows cdef inline linbox_echelonize_efd(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols): + """ + In-place transform this matrix into its reduced row echelon form, and return + the rank `r` of this matrix as well as a list of length `r`, sorted + increasingly. This list gives the column rank profile of this matrix + (which is also that of its reduced row echelon form). + """ # See trac #13878: This is to avoid sending invalid data to linbox, # which would yield a segfault in Sage's debug version. TODO: Fix # that bug upstream. From 2ff6016c5858cd66e1224d402243df30e9cb85fa Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Tue, 19 Aug 2025 16:49:31 +0200 Subject: [PATCH 7/8] fix/augment documentation for echelonize_linbox --- src/sage/matrix/matrix_modn_dense_template.pxi | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 08222ac27e4..79a24b9bd89 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -1690,6 +1690,14 @@ cdef class Matrix_modn_dense_template(Matrix_dense): sage: ~A == C True + sage: B = A.augment(MS(1)) + sage: B.echelonize(algorithm='linbox') + sage: A.rank() + 10 + sage: C = B.submatrix(0,10,10,10) + sage: ~A == C + True + :: sage: A = random_matrix(Integers(10), 10, 20) @@ -1786,7 +1794,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense): sage: for p in (3,17,97,127,1048573): ....: for i in range(10): - ....: A = random_matrix(GF(3), 100, 100) + ....: A = random_matrix(GF(p), 100, 100) ....: A.echelonize(algorithm='all') """ x = self.fetch('in_echelon_form') From 0dfb9490f77a9fa2a51329d9d11c6664e6da1d10 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Wed, 20 Aug 2025 09:57:06 +0200 Subject: [PATCH 8/8] remove unused imports --- src/sage/libs/linbox/fflas.pxd | 3 --- src/sage/matrix/matrix_modn_dense_template.pxi | 11 ++++++----- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/sage/libs/linbox/fflas.pxd b/src/sage/libs/linbox/fflas.pxd index ce8d35e465d..aaedc00c0b0 100644 --- a/src/sage/libs/linbox/fflas.pxd +++ b/src/sage/libs/linbox/fflas.pxd @@ -83,7 +83,6 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFLAS": cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFPACK": ctypedef enum FFPACK_LU_TAG: FfpackTileRecursive - FfpackSlabRecursive void RankProfileFromLU (size_t* P, size_t N, size_t R, size_t* rkprofile, FFPACK_LU_TAG LuTag) @@ -92,8 +91,6 @@ cdef extern from "fflas-ffpack/fflas-ffpack.h" namespace "FFPACK": void MathPerm2LAPACKPerm (size_t * LapackP, size_t * MathP, size_t N) - void LAPACKPerm2MathPerm (size_t * MathP, size_t * LapackP, size_t N) - # double bint IsSingular (Modular_double F, size_t nrows, size_t ncols, Modular_double.Element* A, diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 79a24b9bd89..16e66197ecf 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -94,8 +94,8 @@ from cysignals.signals cimport sig_check, sig_on, sig_off from sage.libs.gmp.mpz cimport * from sage.libs.linbox.fflas cimport FFLAS_TRANSPOSE, FflasNoTrans, FflasTrans, \ - FfpackTileRecursive, FfpackSlabRecursive, FflasLeft, FflasRight, vector, list as std_list, \ - RankProfileFromLU, PLUQtoEchelonPermutation, MathPerm2LAPACKPerm, LAPACKPerm2MathPerm + FfpackTileRecursive, FflasLeft, FflasRight, vector, list as std_list, \ + RankProfileFromLU, PLUQtoEchelonPermutation, MathPerm2LAPACKPerm from libcpp cimport bool from sage.parallel.parallelism import Parallelism @@ -1838,8 +1838,8 @@ cdef class Matrix_modn_dense_template(Matrix_dense): ignore). However, ``efd=True`` uses more memory than FFLAS directly (default: ``True``) - OUTPUT: if ``efd`` is ``False``, return the row rank profile of - ``self`` + OUTPUT: if ``efd`` is ``False``, return the pivot rows (a.k.a. row rank + profile) of the matrix ``self`` before echelonization EXAMPLES:: @@ -1859,7 +1859,8 @@ cdef class Matrix_modn_dense_template(Matrix_dense): r, pivots = linbox_echelonize_efd(self.p, self._entries, self._nrows, self._ncols) else: - r, pivots, rrp = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols) + r, pivots, rrp = linbox_echelonize(self.p, self._entries, + self._nrows, self._ncols) verbose('done with echelonize', t) self.cache('in_echelon_form', True) self.cache('rank', r)