Matrix optimizations

src/meson.build

@@ -288,6 +288,7 @@ add_project_arguments('--embed-positions', language: 'cython')
 add_project_arguments('-X fast_getattr=True', language: 'cython')
 #add_project_arguments('-X language_level="3"', language : 'cython')
 add_project_arguments('-X legacy_implicit_noexcept=True', language: 'cython')
+add_project_arguments('--annotate', language: 'cython')
 add_project_arguments(
   '-X preliminary_late_includes_cy28=True',
   language: 'cython',

src/sage/matrix/docs.py

@@ -391,11 +391,11 @@ class derived from Matrix). They can be either sparse or dense, and
        * cdef _list -- list of underlying elements (need not be a copy)
        * cdef _dict -- sparse dictionary of underlying elements
        * cdef _add_ -- add two matrices with identical parents
-       * _matrix_times_matrix_c_impl -- multiply two matrices with compatible dimensions and
+       * _matrix_times_matrix_ -- multiply two matrices with compatible dimensions and
                                         identical base rings (both sparse or both dense)
        * cpdef _richcmp_ -- compare two matrices with identical parents
-       * cdef _lmul_c_impl -- multiply this matrix on the right by a scalar, i.e., self * scalar
-       * cdef _rmul_c_impl -- multiply this matrix on the left by a scalar, i.e., scalar * self
+       * cdef _lmul_ -- multiply this matrix on the right by a scalar, i.e., self * scalar
+       * cdef _rmul_ -- multiply this matrix on the left by a scalar, i.e., scalar * self
        * __copy__
        * __neg__
 

src/sage/matrix/matrix0.pxd

@@ -24,6 +24,11 @@ cdef class Matrix(sage.structure.element.Matrix):
     cpdef _add_(self, other)
     cpdef _sub_(self, other)
 
+    # These are bint instead of void because exception handling is faster for methods that return int than void methods
+    # TODO: Check if that is actually true. We only care about performance in the happy path.
+    cdef bint _check_matrix_multiplication_sizes(self, sage.structure.element.Matrix right) except -1
+    cdef bint _check_set_to_matrix_product(self, sage.structure.element.Matrix left, sage.structure.element.Matrix right) except -1
+
     cdef bint _will_use_strassen(self, Matrix right) except -2
     cdef bint _will_use_strassen_echelon(self) except -2
     cdef int _strassen_default_cutoff(self, Matrix right) except -2

src/sage/matrix/matrix0.pyx

@@ -5566,6 +5566,31 @@ cdef class Matrix(sage.structure.element.Matrix):
                 ans.set_unsafe(r, c, self.get_unsafe(r, c) * x)
         return ans
 
+    cdef bint _check_matrix_multiplication_sizes(self, sage.structure.element.Matrix right) except -1:
+        if self._ncols != right._nrows:
+            raise ArithmeticError("number of columns of self must equal number of rows of right")
+        return 0
+
+    cdef bint _check_set_to_matrix_product(self, sage.structure.element.Matrix left, sage.structure.element.Matrix right) except -1:
+        if self._nrows != left._nrows or self._ncols != right._ncols:
+            raise ArithmeticError("size of self is not equal to size of left * right")
+        if left._ncols != right._nrows:
+            # We don't call _check_matrix_multiplication_sizes because we want to avoid function-call overhead
+            raise ArithmeticError("number of columns of self must equal number of rows of right")
+        if self._is_immutable:
+            raise ValueError("cannot set the value of an immutable matrix")
+        if self is left or self is right:
+            raise ValueError("cannot set matrix to product involving itself")
+
+        return 0
+
+    cdef _set_matrix_times_matrix_(self, sage.structure.element.Matrix left, sage.structure.element.Matrix right):
+        # Both self and right are matrices with compatible dimensions and base ring.
+        if (<Matrix> left)._will_use_strassen(right):
+            self._set_multiply_strassen(left, right)
+        else:
+            self._set_multiply_classical(left, right)
+
     cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix right):
         r"""
         Return the product of two matrices.
@@ -5734,11 +5759,9 @@ cdef class Matrix(sage.structure.element.Matrix):
             [             0     -x*y + y*x]
             [             0 -x*y^2 + y^2*x]
         """
-        # Both self and right are matrices with compatible dimensions and base ring.
-        if self._will_use_strassen(right):
-            return self._multiply_strassen(right)
-        else:
-            return self._multiply_classical(right)
+        cdef Matrix output = self.new_matrix(self._nrows, right._ncols)
+        output._set_matrix_times_matrix_(self, right)
+        return output
 
     cdef bint _will_use_strassen(self, Matrix right) except -2:
         """

src/sage/matrix/matrix2.pxd

@@ -22,3 +22,4 @@ cdef class Matrix(Matrix1):
     cpdef _echelon(self, str algorithm)
     cpdef _echelon_in_place(self, str algorithm)
     cpdef matrix_window(self, Py_ssize_t row=*, Py_ssize_t col=*, Py_ssize_t nrows=*, Py_ssize_t ncols=*, bint check=*)
+    cpdef _multiply_strassen(self, Matrix right, int cutoff=*)

src/sage/matrix/matrix2.pyx

@@ -108,7 +108,6 @@ from sage.matrix.matrix_misc import permanental_minor_polynomial
 
 from sage.misc.misc_c import prod
 
-
 SymbolicRing = LazyImport('sage.symbolic.ring', 'SymbolicRing')
 
 
@@ -9399,7 +9398,7 @@ cdef class Matrix(Matrix1):
     # Precise algorithms invented and implemented by David Harvey and Robert Bradshaw
     # at William Stein's MSRI 2006 Summer Workshop on Modular Forms.
     #####################################################################################
-    def _multiply_strassen(self, Matrix right, int cutoff=0):
+    cpdef _multiply_strassen(self, Matrix right, int cutoff=0):
         """
         Multiply ``self`` by the matrix right using a Strassen-based
         asymptotically fast arithmetic algorithm.
@@ -9421,10 +9420,9 @@ cdef class Matrix(Matrix1):
             [248 286 324 362]
             [344 398 452 506]
         """
-        if self._ncols != right._nrows:
-            raise ArithmeticError("Number of columns of self must equal number of rows of right.")
+        self._check_matrix_multiplication_sizes(right)
         if self._base_ring is not right.base_ring():
-            raise TypeError("Base rings must be the same.")
+            raise TypeError("base rings must be the same")
 
         if cutoff == 0:
             cutoff = self._strassen_default_cutoff(right)

src/sage/matrix/matrix_complex_ball_dense.pyx

@@ -496,17 +496,26 @@ cdef class Matrix_complex_ball_dense(Matrix_dense):
         """
         return self._lmul_(a)
 
-    cdef _matrix_times_matrix_(self, Matrix other):
+    cdef _set_matrix_times_matrix_(self, Matrix left, Matrix other):
         r"""
         TESTS::
 
             sage: matrix(CBF, [[1,2]])*matrix([[3], [4]]) # indirect doctest
             [11.00000000000000]
         """
-        cdef Matrix_complex_ball_dense res = self._new(self._nrows, other._ncols)
         sig_on()
-        acb_mat_mul(res.value, self.value, (<Matrix_complex_ball_dense> other).value, prec(self))
+        acb_mat_mul(self.value, (<Matrix_complex_ball_dense> left).value, (<Matrix_complex_ball_dense> other).value, prec(left))
         sig_off()
+
+    cdef Matrix _matrix_times_matrix_(self, Matrix other):
+        r"""
+        TESTS::
+
+            sage: matrix(CBF, [[1,2]])*matrix([[3], [4]]) # indirect doctest
+            [11.00000000000000]
+        """
+        cdef Matrix_complex_ball_dense res = self._new(self._nrows, other._ncols)
+        res._set_matrix_times_matrix_(self, other)
         return res
 
     cpdef _pow_int(self, n):

src/sage/matrix/matrix_cyclo_dense.pyx

@@ -635,7 +635,7 @@ cdef class Matrix_cyclo_dense(Matrix_dense):
             A._matrix = T * self._matrix
         return A
 
-    cdef _matrix_times_matrix_(self, baseMatrix right):
+    cdef _set_matrix_times_matrix_(self, baseMatrix left, baseMatrix right):
         """
         Return the product of two cyclotomic dense matrices.
 
@@ -694,40 +694,43 @@ cdef class Matrix_cyclo_dense(Matrix_dense):
             sage: (-m)*n
             [-23250]
         """
-        A, denom_self = self._matrix._clear_denom()
-        B, denom_right = (<Matrix_cyclo_dense>right)._matrix._clear_denom()
+        A, denom_left = (<Matrix_cyclo_dense> left)._matrix._clear_denom()
+        B, denom_right = (<Matrix_cyclo_dense> right)._matrix._clear_denom()
 
         # conservative but correct estimate: 2 is there to account for the
         # sign of the entries
-        bound = 1 + 2 * A.height() * B.height() * self._ncols
+        bound = 1 + 2 * A.height() * B.height() * left._ncols
 
-        n = self._base_ring._n()
+        n = left._base_ring._n()
         p = previous_prime(MAX_MODULUS)
         prod = 1
         v = []
         while prod <= bound:
-            while (n >= 2 and p % n != 1) or denom_self % p == 0 or denom_right % p == 0:
+            while (n >= 2 and p % n != 1) or denom_left % p == 0 or denom_right % p == 0:
                 if p == 2:
-                    raise RuntimeError("we ran out of primes in matrix multiplication.")
+                    raise RuntimeError("we ran out of primes in matrix multiplication")
                 p = previous_prime(p)
             prod *= p
-            Amodp, _ = self._reductions(p)
+            Amodp, _ = left._reductions(p)
             Bmodp, _ = right._reductions(p)
-            _, S = self._reduction_matrix(p)
+            _, S = left._reduction_matrix(p)
             X = Amodp[0]._matrix_from_rows_of_matrices([Amodp[i] * Bmodp[i] for i in range(len(Amodp))])
             v.append(S*X)
             p = previous_prime(p)
-        M = matrix(ZZ, self._base_ring.degree(), self._nrows*right.ncols())
+        M = matrix(ZZ, left._base_ring.degree(), left._nrows*right.ncols())
         _lift_crt(M, v)
-        d = denom_self * denom_right
+        d = denom_left * denom_right
         if d == 1:
             M = M.change_ring(QQ)
         else:
             M = (1/d)*M
+        self._matrix = M
+
+    cdef Matrix_cyclo_dense _matrix_times_matrix_(self, baseMatrix right):
         cdef Matrix_cyclo_dense C = Matrix_cyclo_dense.__new__(Matrix_cyclo_dense,
-                    MatrixSpace(self._base_ring, self._nrows, right.ncols()),
+                                                               MatrixSpace(self._base_ring, self._nrows, right.ncols()),
                                                                None, None, None)
-        C._matrix = M
+        C._set_matrix_times_matrix_(self, right)
         return C
 
     cdef Py_hash_t _hash_(self) except -1:

src/sage/matrix/matrix_dense.pxd

@@ -2,3 +2,6 @@ from sage.matrix.matrix cimport Matrix
 
 cdef class Matrix_dense(Matrix):
     cdef void set_unsafe_int(self, Py_ssize_t i, Py_ssize_t j, int value) noexcept
+    cdef _set_multiply_strassen(self, Matrix left, Matrix right, int cutoff=*)
+    cdef _set_multiply_classical(self, Matrix left, Matrix right)
+    cpdef _multiply_classical(self, Matrix right)

src/sage/matrix/matrix_dense.pyx

@@ -8,14 +8,12 @@ TESTS::
     sage: TestSuite(m).run(skip='_test_minpoly')
 """
 
-cimport sage.matrix.matrix as matrix
-
 from sage.structure.richcmp cimport richcmp_item, rich_to_bool
 import sage.matrix.matrix_space
 import sage.structure.sequence
 
 
-cdef class Matrix_dense(matrix.Matrix):
+cdef class Matrix_dense(Matrix):
     cdef bint is_sparse_c(self) noexcept:
         return 0
 
@@ -297,9 +295,15 @@ cdef class Matrix_dense(matrix.Matrix):
             image.subdivide(*self.subdivisions())
         return image
 
-    def _multiply_classical(left, matrix.Matrix right):
+    cpdef _multiply_classical(self, Matrix right):
+        self._check_matrix_multiplication_sizes(right)
+        cdef Matrix_dense res = self.new_matrix(nrows=self._nrows, ncols=right._ncols)
+        res._set_multiply_classical(self, right)
+        return res
+
+    cdef _set_multiply_classical(self, Matrix left, Matrix right):
         """
-        Multiply the matrices left and right using the classical `O(n^3)`
+        Multiply the matrices self and right using the classical `O(n^3)`
         algorithm.
 
         This method will almost always be overridden either by the
@@ -325,17 +329,36 @@ cdef class Matrix_dense(matrix.Matrix):
             sage: Matrix_dense._multiply_classical(matrix(2, 1), matrix(2, 0))
             Traceback (most recent call last):
             ...
-            ArithmeticError: number of columns of left must equal number of rows of right
+            ArithmeticError: number of columns of self must equal number of rows of right
         """
         cdef Py_ssize_t i, j, k
-        if left._ncols != right._nrows:
-            raise ArithmeticError("number of columns of left must equal number of rows of right")
-        zero = left.base_ring().zero()
-        cdef matrix.Matrix res = left.new_matrix(nrows=left._nrows, ncols=right._ncols)
-        for i in range(left._nrows):
+        zero = self.base_ring().zero()
+        for i in range(self._nrows):
             for j in range(right._ncols):
                 dotp = zero
-                for k in range(left._ncols):
-                    dotp += left.get_unsafe(i, k) * right.get_unsafe(k, j)
-                res.set_unsafe(i, j, dotp)
-        return res
+                for k in range(self._ncols):
+                    dotp += self.get_unsafe(i, k) * right.get_unsafe(k, j)
+                self.set_unsafe(i, j, dotp)
+
+    cpdef _multiply_strassen(self, Matrix right, int cutoff=0):
+        self._check_matrix_multiplication_sizes(right)
+        if self._base_ring is not right.base_ring():
+            raise TypeError("base rings must be the same")
+
+        cdef Matrix_dense output = self.new_matrix(self._nrows, right._ncols)
+        output._set_multiply_strassen(self, right, cutoff)
+        return output
+
+    cdef _set_multiply_strassen(self, Matrix left, Matrix right, int cutoff=0):
+        if cutoff == 0:
+            cutoff = self._strassen_default_cutoff(right)
+
+        if cutoff <= 0:
+            raise ValueError("cutoff must be at least 1")
+
+        self_window = self.matrix_window()
+        left_window = left.matrix_window()
+        right_window = right.matrix_window()
+
+        from sage.matrix import strassen
+        strassen.strassen_window_multiply(self_window, left_window, right_window, cutoff)

src/sage/matrix/matrix_double_dense.pyx

@@ -220,7 +220,19 @@ cdef class Matrix_double_dense(Matrix_numpy_dense):
     # def _pickle(self):                        #unsure how to implement
     # def _unpickle(self, data, int version):   # use version >= 0 #unsure how to implement
     ######################################################################
-    cdef sage.structure.element.Matrix _matrix_times_matrix_(self, sage.structure.element.Matrix right):
+
+    cdef _set_matrix_times_matrix_(self, sage.structure.element.Matrix left, sage.structure.element.Matrix right):
+        if left._nrows == 0 or left._ncols == 0 or right._nrows == 0 or right._ncols == 0:  # TODO: Do we even need this check?
+            self._matrix_numpy.fill(0)
+            return
+
+        global numpy
+        if numpy is None:
+            import numpy
+
+        self._matrix_numpy = numpy.dot((<Matrix_double_dense> left)._matrix_numpy, (<Matrix_double_dense> right)._matrix_numpy)
+
+    cdef Matrix_double_dense _matrix_times_matrix_(self, sage.structure.element.Matrix right):
         r"""
         Multiply ``self * right`` as matrices.
 
@@ -248,24 +260,9 @@ cdef class Matrix_double_dense(Matrix_numpy_dense):
             [0.0 0.0 0.0]
             [0.0 0.0 0.0]
         """
-        if self._ncols != right._nrows:
-            raise IndexError("Number of columns of self must equal number of rows of right")
-
-        cdef Matrix_double_dense M, _right, _left
-
-        if self._nrows == 0 or self._ncols == 0 or right._nrows == 0 or right._ncols == 0:
-            M = self._new(self._nrows, right._ncols)
-            M._matrix_numpy.fill(0)
-            return M
-
+        self._check_matrix_multiplication_sizes(right)
         M = self._new(self._nrows, right._ncols)
-        _right = right
-        _left = self
-        global numpy
-        if numpy is None:
-            import numpy
-
-        M._matrix_numpy = numpy.dot(_left._matrix_numpy, _right._matrix_numpy)
+        M._set_matrix_times_matrix_(self, right)
         return M
 
     def __invert__(self):

src/sage/matrix/matrix_gap.pyx

@@ -347,7 +347,10 @@ cdef class Matrix_gap(Matrix_dense):
         ans._libgap = left._libgap - (<Matrix_gap> right)._libgap
         return ans
 
-    cdef Matrix _matrix_times_matrix_(left, Matrix right):
+    cdef _set_matrix_times_matrix_(self, Matrix left, Matrix right):
+        self._libgap = <Matrix_gap> ((<Matrix_gap> left)._libgap * (<Matrix_gap> right)._libgap)
+
+    cdef Matrix_gap _matrix_times_matrix_(self, Matrix right):
         r"""
         TESTS::
 
@@ -358,10 +361,10 @@ cdef class Matrix_gap(Matrix_dense):
             [ 1 -1]
             [ 7 -7]
         """
-        if left._ncols != right._nrows:
-            raise IndexError("Number of columns of self must equal number of rows of right.")
-        cdef Matrix_gap M = left._new(left._nrows, right._ncols)
-        M._libgap = <Matrix_gap> ((<Matrix_gap> left)._libgap * (<Matrix_gap> right)._libgap)
+        self._check_matrix_multiplication_sizes(right)
+        cdef Matrix_gap M = self._new(self._nrows, right._ncols)
+        M._set_matrix_times_matrix_(self, right)
+        #M._libgap = <Matrix_gap> ((<Matrix_gap> self)._libgap * (<Matrix_gap> right)._libgap)
         return M
 
     def transpose(self):