gh-36059: Speed up the creation of submatrices of `Matrix_modn_dense_template` matrices
    
<!-- ^^^^^
Please provide a concise, informative and self-explanatory title.
Don't put issue numbers in there, do this in the PR body below.
For example, instead of "Fixes #1234" use "Introduce new method to
calculate 1+1"
-->
<!-- Describe your changes here in detail -->
Modified the method `submatrix()` to do a systematic call to `memcpy()`.
Overrode  methods `matrix_from_rows()`, `matrix_from_columns()` and
`matrix_from_rows_and_columns()`.
<!-- Why is this change required? What problem does it solve? -->
In the current version of `submatrix()`, only the case when the
submatrix has as many columns as the matrix,  relies on the contiguous
memory storage of `Matrix_modn_dense_template` and call `memcpy()`.
In fact, this can be applied to any set of arguments given to the method
and enables better performances.

Overriding the other methods mentioned above allows to avoid calling
`get_unsafe()` and `set_unsafe()` that are doing time-consuming
casts/conversions that are unnecessary in this case.
<!-- If this PR resolves an open issue, please link to it here. For
example "Fixes #12345". -->
<!-- If your change requires a documentation PR, please link it
appropriately. -->

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->
<!-- If your change requires a documentation PR, please link it
appropriately -->
<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
<!-- Feel free to remove irrelevant items. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [ ] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [ ] I have updated the documentation accordingly.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on
- #12345: short description why this is a dependency
- #34567: ...
-->

<!-- If you're unsure about any of these, don't hesitate to ask. We're
here to help! -->
    
URL: #36059
Reported by: Marie Bonboire
Reviewer(s): Vincent Neiger

Author: Release Manager

src/sage/matrix/matrix1.pyx

@@ -2176,9 +2176,7 @@ cdef class Matrix(Matrix0):
             [5 4]
 
         For example here we take from row 1 columns 2 then 0 twice, and do
-        this 3 times.
-
-        ::
+        this 3 times::
 
             sage: A.matrix_from_rows_and_columns([1,1,1],[2,0,0])
             [5 3 3]

src/sage/matrix/matrix_integer_dense.pyx

@@ -1614,7 +1614,7 @@ cdef class Matrix_integer_dense(Matrix_dense):
             return self._mod_two()
         elif p < MAX_MODULUS_FLOAT:
             res_f = Matrix_modn_dense_float.__new__(Matrix_modn_dense_float,
-                                                    matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None)
+                                                    matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None, zeroed_alloc=False)
             for i from 0 <= i < self._nrows:
                 res_row_f = res_f._matrix[i]
                 for j from 0 <= j < self._ncols:
@@ -1623,7 +1623,7 @@ cdef class Matrix_integer_dense(Matrix_dense):
 
         elif p < MAX_MODULUS_DOUBLE:
             res_d = Matrix_modn_dense_double.__new__(Matrix_modn_dense_double,
-                                                     matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None)
+                                                     matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False), None, None, None, zeroed_alloc=False)
             for i from 0 <= i < self._nrows:
                 res_row_d = res_d._matrix[i]
                 for j from 0 <= j < self._ncols:
@@ -1649,11 +1649,11 @@ cdef class Matrix_integer_dense(Matrix_dense):
             if p < MAX_MODULUS_FLOAT:
                 res.append( Matrix_modn_dense_float.__new__(Matrix_modn_dense_float,
                                                             matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False),
-                                                            None, None, None) )
+                                                            None, None, None, zeroed_alloc=False) )
             elif p < MAX_MODULUS_DOUBLE:
                 res.append( Matrix_modn_dense_double.__new__(Matrix_modn_dense_double,
                                                              matrix_space.MatrixSpace(IntegerModRing(p), self._nrows, self._ncols, sparse=False),
-                                                             None, None, None) )
+                                                             None, None, None, zeroed_alloc=False) )
             else:
                 raise ValueError("p=%d too big."%p)
 

src/sage/matrix/matrix_modn_dense_template.pxi

@@ -90,7 +90,7 @@ We test corner cases for multiplication::
 from libc.stdint cimport uint64_t
 from cpython.bytes cimport *
 
-from cysignals.memory cimport check_malloc, check_allocarray, sig_malloc, sig_free
+from cysignals.memory cimport check_malloc, check_allocarray, check_calloc, sig_malloc, sig_free
 from cysignals.signals cimport sig_check, sig_on, sig_off
 
 from sage.libs.gmp.mpz cimport *
@@ -123,7 +123,7 @@ from sage.structure.proof.proof import get_flag as get_proof_flag
 from sage.structure.richcmp cimport rich_to_bool
 from sage.misc.randstate cimport randstate, current_randstate
 import sage.matrix.matrix_space as matrix_space
-from .args cimport MatrixArgs_init
+from .args cimport SparseEntry, MatrixArgs_init
 
 
 from sage.cpython.string cimport char_to_str
@@ -441,15 +441,18 @@ cpdef __matrix_from_rows_of_matrices(X):
 
 
 cdef class Matrix_modn_dense_template(Matrix_dense):
-    def __cinit__(self):
+    def __cinit__(self, *args, bint zeroed_alloc=True, **kwds):
         cdef long p = self._base_ring.characteristic()
         self.p = p
         if p >= MAX_MODULUS:
             raise OverflowError("p (=%s) must be < %s."%(p, MAX_MODULUS))
 
-        self._entries = <celement *>check_allocarray(self._nrows * self._ncols, sizeof(celement))
+        if zeroed_alloc:
+            self._entries = <celement *>check_calloc(self._nrows * self._ncols, sizeof(celement))
+        else:
+            self._entries = <celement *>check_allocarray(self._nrows * self._ncols, sizeof(celement))
+            
         self._matrix = <celement **>check_allocarray(self._nrows, sizeof(celement*))
-
         cdef unsigned int k
         cdef Py_ssize_t i
         k = 0
@@ -518,27 +521,28 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
         """
         ma = MatrixArgs_init(parent, entries)
         cdef long i, j
-        it = ma.iter(False)
+        it = ma.iter(convert=False, sparse=True)
         R = ma.base
         p = R.characteristic()
-        for i in range(ma.nrows):
-            v = self._matrix[i]
-            for j in range(ma.ncols):
-                x = next(it)
-                if type(x) is int:
-                    tmp = (<long>x) % p
-                    v[j] = tmp + (tmp<0)*p
-                elif type(x) is IntegerMod_int and (<IntegerMod_int>x)._parent is R:
-                    v[j] = <celement>(<IntegerMod_int>x).ivalue
-                elif type(x) is Integer:
-                    if coerce:
-                        v[j] = mpz_fdiv_ui((<Integer>x).value, p)
-                    else:
-                        v[j] = mpz_get_ui((<Integer>x).value)
-                elif coerce:
-                    v[j] = R(x)
+        
+        for t in it:
+            se = <SparseEntry>t
+            x = se.entry
+            v = self._matrix[se.i]
+            if type(x) is int:
+                tmp = (<long>x) % p
+                v[se.j] = tmp + (tmp<0)*p
+            elif type(x) is IntegerMod_int and (<IntegerMod_int>x)._parent is R:
+                v[se.j] = <celement>(<IntegerMod_int>x).ivalue
+            elif type(x) is Integer:
+                if coerce:
+                    v[se.j] = mpz_fdiv_ui((<Integer>x).value, p)
                 else:
-                    v[j] = <celement>x
+                    v[se.j] = mpz_get_ui((<Integer>x).value)
+            elif coerce:
+                v[se.j] = R(x)
+            else:
+                v[se.j] = <celement>x
 
     cdef long _hash_(self) except -1:
         """
@@ -786,7 +790,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
         cdef Matrix_modn_dense_template M
         cdef celement p = self.p
 
-        M = self.__class__.__new__(self.__class__, self._parent,None,None,None)
+        M = self.__class__.__new__(self.__class__, self._parent,None,None,None, zeroed_alloc=False)
 
         sig_on()
         for i in range(self._nrows*self._ncols):
@@ -825,7 +829,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
         cdef celement p = self.p
         cdef celement a = left
 
-        M = self.__class__.__new__(self.__class__, self._parent,None,None,None)
+        M = self.__class__.__new__(self.__class__, self._parent,None,None,None,zeroed_alloc=False)
 
         sig_on()
         for i in range(self._nrows*self._ncols):
@@ -844,7 +848,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
             False
         """
         cdef Matrix_modn_dense_template A
-        A = self.__class__.__new__(self.__class__, self._parent, 0, 0, 0)
+        A = self.__class__.__new__(self.__class__,self._parent,None,None,None,zeroed_alloc=False)
         memcpy(A._entries, self._entries, sizeof(celement)*self._nrows*self._ncols)
         if self._subdivisions is not None:
             A.subdivide(*self.subdivisions())
@@ -883,7 +887,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
         cdef celement k, p
         cdef Matrix_modn_dense_template M
 
-        M = self.__class__.__new__(self.__class__, self._parent,None,None,None)
+        M = self.__class__.__new__(self.__class__, self._parent,None,None,None,zeroed_alloc=False)
         p = self.p
         cdef celement* other_ent = (<Matrix_modn_dense_template>right)._entries
 
@@ -920,7 +924,7 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
         cdef celement k, p
         cdef Matrix_modn_dense_template M
 
-        M = self.__class__.__new__(self.__class__, self._parent, None, None, None)
+        M = self.__class__.__new__(self.__class__, self._parent, None, None, None, zeroed_alloc=False)
         p = self.p
         cdef celement* other_ent = (<Matrix_modn_dense_template>right)._entries
 
@@ -3012,14 +3016,21 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
         if nrows == -1:
             nrows = self._nrows - row
 
-        if col != 0 or ncols != self._ncols:
-            return self.matrix_from_rows_and_columns(range(row, row+nrows), range(col, col+ncols))
-
         if nrows < 0 or row < 0 or row + nrows > self._nrows:
             raise IndexError("rows out of range")
+        if ncols < 0 or col < 0 or col + ncols > self._ncols:
+            raise IndexError("columns out of range")
+
+        cdef Matrix_modn_dense_template M = self.new_matrix(nrows=nrows, ncols=ncols)
+
+        if col == 0 and ncols == self._ncols:
+            memcpy(M._entries, self._matrix[row], sizeof(celement)*ncols*nrows)
+            return M
+
+        cdef Py_ssize_t i,r
+        for i,r in enumerate(range(row, row+nrows)) :
+            memcpy(M._matrix[i], self._matrix[r]+col, sizeof(celement)*ncols)
 
-        cdef Matrix_modn_dense_template M = self.new_matrix(nrows=nrows, ncols=self._ncols)
-        memcpy(M._entries, self._entries+row*ncols, sizeof(celement)*ncols*nrows)
         return M
 
     def _matrices_from_rows(self, Py_ssize_t nrows, Py_ssize_t ncols):
@@ -3065,6 +3076,122 @@ cdef class Matrix_modn_dense_template(Matrix_dense):
             ans.append(M)
         return ans
 
+    def matrix_from_columns(self, columns):
+        """
+        Return the matrix constructed from self using columns with indices
+        in the columns list.
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(Integers(8),3,3)
+            sage: A = M(range(9)); A
+            [0 1 2]
+            [3 4 5]
+            [6 7 0]
+            sage: A.matrix_from_columns([2,1])
+            [2 1]
+            [5 4]
+            [0 7]
+        """
+        cdef Py_ssize_t ncols = len(columns)
+
+        # Construct new matrix
+        cdef Matrix_modn_dense_template A = self.new_matrix(ncols=ncols)
+        cdef Py_ssize_t i, j, col
+        for j, col in enumerate(columns):
+            if col < 0 or col >= self._ncols:
+                raise IndexError("column index out of range")
+            for i in range(self._nrows):
+                A._matrix[i][j] = self._matrix[i][col]
+
+        return A
+
+    def matrix_from_rows(self, rows):
+        """
+        Return the matrix constructed from self using rows with indices in
+        the rows list.
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(Integers(8),3,3)
+            sage: A = M(range(9)); A
+            [0 1 2]
+            [3 4 5]
+            [6 7 0]
+            sage: A.matrix_from_rows([2,1])
+            [6 7 0]
+            [3 4 5]
+        """
+        cdef Py_ssize_t nrows = len(rows)
+
+        # Construct new matrix
+        cdef Matrix_modn_dense_template A = self.new_matrix(nrows=nrows)
+
+        cdef Py_ssize_t i, row
+        for i, row in enumerate(rows):
+            if row < 0 or row >= self._nrows:
+                raise IndexError("row index out of range")
+            memcpy(A._matrix[i], self._matrix[row], sizeof(celement)*self._ncols)
+
+        return A
+
+    def matrix_from_rows_and_columns(self, rows, columns):
+        """
+        Return the matrix constructed from self from the given rows and
+        columns.
+
+        EXAMPLES::
+
+            sage: M = MatrixSpace(Integers(8),3,3)
+            sage: A = M(range(9)); A
+            [0 1 2]
+            [3 4 5]
+            [6 7 0]
+            sage: A.matrix_from_rows_and_columns([1], [0,2])
+            [3 5]
+            sage: A.matrix_from_rows_and_columns([1,2], [1,2])
+            [4 5]
+            [7 0]
+
+        Note that row and column indices can be reordered or repeated::
+
+            sage: A.matrix_from_rows_and_columns([2,1], [2,1])
+            [0 7]
+            [5 4]
+
+        For example here we take from row 1 columns 2 then 0 twice, and do
+        this 3 times::
+
+            sage: A.matrix_from_rows_and_columns([1,1,1],[2,0,0])
+            [5 3 3]
+            [5 3 3]
+            [5 3 3]
+
+        AUTHORS:
+
+        - Jaap Spies (2006-02-18)
+
+        - Didier Deshommes: some Pyrex speedups implemented
+        """
+        cdef Py_ssize_t ncols = len(columns)
+        cdef Py_ssize_t nrows = len(rows)
+
+        # Check whether column indices are valid
+        cdef Py_ssize_t i, j, row, col
+        for col in columns:
+            if col < 0 or col >= self._ncols:
+                raise IndexError("column index out of range")
+
+        # Construct new matrix
+        cdef Matrix_modn_dense_template A = self.new_matrix(nrows=nrows, ncols=ncols)
+        for i, row in enumerate(rows):
+            if row < 0 or row >= self._nrows:
+                raise IndexError("row index out of range")
+            for j, col in enumerate(columns):
+                A._matrix[i][j] = self._matrix[row][col]
+
+        return A
+
     def __bool__(self):
         """
         Test whether this matrix is zero.