Trac #23706: allow several implementations of matrices in MatrixSpace

The aim of this ticket is to revamp the `MatrixSpace` constructor in
order to be able to choose the implementation.

For example, with the attached branch one can have access to 2
implementations of integer matrices
{{{
sage: MatrixSpace(ZZ, 3, implementation='flint')
Full MatrixSpace of 3 by 3 dense matrices over Integer Ring
sage: MatrixSpace(ZZ, 3, implementation='generic')
Full MatrixSpace of 3 by 3 dense matrices over Integer Ring (using
Matrix_generic_dense)
}}}

One important change is the move of the method `_get_matrix_class` of
`MatrixSpace` as an independent function.

A problem arose with `CartanMatrix` (in
`sage/combinat/root_system/cartan_matrix.py`). The class `CartanMatrix`
inherits from `Matrix_integer_sparse` but does not properly redefines
the parent. As a consequence of the changes here, many matrix operations
are not valid anymore for inherited types (mostly submatrices operations
such as `__getitem__` with slices or `stack`). A fix is provided by
defining `CartanMatrix.matrix_space`.

follow up: #23714 (implementation of gap matrices)

URL: https://trac.sagemath.org/23706
Reported by: vdelecroix
Ticket author(s): Vincent Delecroix
Reviewer(s): Jean-Philippe Labbé, Travis Scrimshaw

src/sage/algebras/clifford_algebra.py

@@ -1907,7 +1907,7 @@ def lifted_form(x, y):
                     if m != n:
                         continue
                     MS = MatrixSpace(R, n, n)
-                    MC = MS._get_matrix_class()
+                    MC = MS._matrix_class
                     matrix_list = [M[mx[i], my[j]]
                                    for i in range(n)
                                    for j in range(n)]

src/sage/combinat/root_system/cartan_matrix.py

@@ -295,6 +295,42 @@ def __classcall_private__(cls, data=None, index_set=None,
         mat._subdivisions = subdivisions
         return mat
 
+    def matrix_space(self, nrows=None, ncols=None, sparse=None):
+        r"""
+        Return a matrix space over the integers.
+
+        INPUT:
+
+        - ``nrows`` - number of rows
+
+        - ``ncols`` - number of columns
+
+        - ``sparse`` - (boolean) sparseness
+
+        EXAMPLES::
+
+            sage: cm = CartanMatrix(['A', 3])
+            sage: cm.matrix_space()
+            Full MatrixSpace of 3 by 3 sparse matrices over Integer Ring
+            sage: cm.matrix_space(2, 2)
+            Full MatrixSpace of 2 by 2 sparse matrices over Integer Ring
+            sage: cm[:2,1:]   # indirect doctest
+            [-1  0]
+            [ 2 -1]
+        """
+        if nrows is None:
+            nrows = self.nrows()
+        if ncols is None:
+            ncols = self.ncols()
+        if sparse is None:
+            sparse = True
+
+        if nrows == self.nrows() and ncols == self.ncols() and sparse:
+            return self.parent()
+        else:
+            from sage.matrix.matrix_space import MatrixSpace
+            return MatrixSpace(ZZ, nrows, ncols, sparse is None or bool(sparse))
+
     def _CM_init(self, cartan_type, index_set, cartan_type_check):
         """
         Initialize ``self`` as a Cartan matrix.
@@ -674,11 +710,11 @@ def row_with_indices(self, i):
     def is_finite(self):
         """
         Return ``True`` if ``self`` is a finite type or ``False`` otherwise.
-        
+
         A generalized Cartan matrix is finite if the determinant of all its
-        principal submatrices (see :meth:`principal_submatrices`) is positive. 
-        Such matrices have a positive definite symmetrized matrix. Note that a 
-        finite matrix may consist of multiple blocks of Cartan matrices each 
+        principal submatrices (see :meth:`principal_submatrices`) is positive.
+        Such matrices have a positive definite symmetrized matrix. Note that a
+        finite matrix may consist of multiple blocks of Cartan matrices each
         having finite Cartan type.
 
         EXAMPLES::
@@ -696,16 +732,16 @@ def is_finite(self):
         if self._cartan_type is None:
             if not self.is_symmetrizable():
                 return False
-            return self.symmetrized_matrix().is_positive_definite() 
+            return self.symmetrized_matrix().is_positive_definite()
         return self._cartan_type.is_finite()
 
     @cached_method
     def is_affine(self):
         """
         Return ``True`` if ``self`` is an affine type or ``False`` otherwise.
-        
-        A generalized Cartan matrix is affine if all of its indecomposable 
-        blocks are either finite (see :meth:`is_finite`) or have zero 
+
+        A generalized Cartan matrix is affine if all of its indecomposable
+        blocks are either finite (see :meth:`is_finite`) or have zero
         determinant with all proper principal minors positive.
 
         EXAMPLES::
@@ -725,26 +761,26 @@ def is_affine(self):
                 return False
             for b in self.indecomposable_blocks():
                 if b.det() < 0 or not all(
-                    a.det() > 0 for a in b.principal_submatrices(proper=True)): 
+                    a.det() > 0 for a in b.principal_submatrices(proper=True)):
                     return False
             return True
         return self._cartan_type.is_affine()
-    
+
     @cached_method
     def is_hyperbolic(self, compact=False):
         """
-        Return if ``True`` if ``self`` is a (compact) hyperbolic type 
+        Return if ``True`` if ``self`` is a (compact) hyperbolic type
         or ``False`` otherwise.
-        
+
         An indecomposable generalized Cartan matrix is hyperbolic if it has
         negative determinant and if any proper connected subdiagram of its
         Dynkin diagram is of finite or affine type. It is compact hyperbolic
         if any proper connected subdiagram has finite type.
-        
+
         INPUT:
 
-        - ``compact`` -- if ``True``, check if matrix is compact hyperbolic  
-        
+        - ``compact`` -- if ``True``, check if matrix is compact hyperbolic
+
         EXAMPLES::
 
             sage: M = CartanMatrix([[2,-2,0],[-2,2,-1],[0,-1,2]])
@@ -761,7 +797,7 @@ def is_hyperbolic(self, compact=False):
         """
         if not self.is_indefinite() or not self.is_indecomposable():
             return False
-        
+
         D = self.dynkin_diagram()
         verts = tuple(D.vertex_iterator())
         for v in verts:
@@ -772,15 +808,15 @@ def is_hyperbolic(self, compact=False):
             elif not subg.is_finite() and not subg.is_affine():
                 return False
         return True
-    
+
     @cached_method
     def is_lorentzian(self):
         """
         Return ``True`` if ``self`` is a Lorentzian type or ``False`` otherwise.
-        
+
         A generalized Cartan matrix is Lorentzian if it has negative determinant
         and exactly one negative eigenvalue.
-        
+
         EXAMPLES::
 
             sage: M = CartanMatrix([[2,-3],[-3,2]])
@@ -793,12 +829,12 @@ def is_lorentzian(self):
         if self.det() >= 0:
             return False
         return sum(1 for x in self.eigenvalues() if x < 0) == 1
-        
-    @cached_method        
+
+    @cached_method
     def is_indefinite(self):
         """
         Return if ``self`` is an indefinite type or ``False`` otherwise.
-        
+
         EXAMPLES::
 
            sage: M = CartanMatrix([[2,-3],[-3,2]])
@@ -809,12 +845,12 @@ def is_indefinite(self):
            False
         """
         return not self.is_finite() and not self.is_affine()
-                
+
     @cached_method
     def is_indecomposable(self):
         """
         Return if ``self`` is an indecomposable matrix or ``False`` otherwise.
-        
+
         EXAMPLES::
 
             sage: M = CartanMatrix(['A',5])
@@ -831,11 +867,11 @@ def is_indecomposable(self):
     def principal_submatrices(self, proper=False):
         """
         Return a list of all principal submatrices of ``self``.
-        
+
         INPUT:
 
-        - ``proper`` -- if ``True``, return only proper submatrices 
-        
+        - ``proper`` -- if ``True``, return only proper submatrices
+
         EXAMPLES::
 
             sage: M = CartanMatrix(['A',2])
@@ -846,20 +882,20 @@ def principal_submatrices(self, proper=False):
             ]
             sage: M.principal_submatrices(proper=True)
             [[], [2], [2]]
-            
+
         """
         iset = list(range(self.ncols()))
         ret = []
         for l in powerset(iset):
             if not proper or (proper and l != iset):
                 ret.append(self.matrix_from_rows_and_columns(l,l))
         return ret
-    
+
     @cached_method
     def indecomposable_blocks(self):
         """
         Return a tuple of all indecomposable blocks of ``self``.
-        
+
         EXAMPLES::
 
             sage: M = CartanMatrix(['A',2])
@@ -871,7 +907,7 @@ def indecomposable_blocks(self):
             sage: M = CartanMatrix([['A',2,1],['A',3,1]])
             sage: M.indecomposable_blocks()
             (
-            [ 2 -1  0 -1]            
+            [ 2 -1  0 -1]
             [-1  2 -1  0]  [ 2 -1 -1]
             [ 0 -1  2 -1]  [-1  2 -1]
             [-1  0 -1  2], [-1 -1  2]

src/sage/combinat/root_system/coxeter_type.py

@@ -411,7 +411,7 @@ def val(x):
                     return R(x)
 
         MS = MatrixSpace(R, n, sparse=True)
-        MC = MS._get_matrix_class()
+        MC = MS._matrix_class
 
         bilinear = MC(MS, entries={(i, j): val(mat[i, j])
                                    for i in range(n) for j in range(n)

src/sage/geometry/hyperplane_arrangement/arrangement.py

@@ -1981,7 +1981,7 @@ def face_semigroup_algebra(self, field=None, names='e'):
         # Some hackery to generate a matrix quickly and without
         # unnecessary sanitization/ducktyping:
         MS = MatrixSpace(field, N, N)
-        MC = MS._get_matrix_class()
+        MC = MS._matrix_class
         table = []
         for j, sj in enumerate(Fs):
             matrix_j = []

src/sage/groups/matrix_gps/coxeter_group.py

@@ -278,7 +278,7 @@ def __init__(self, coxeter_matrix, base_ring, index_set):
         n = coxeter_matrix.rank()
         # Compute the matrix with entries `2 \cos( \pi / m_{ij} )`.
         MS = MatrixSpace(base_ring, n, sparse=True)
-        MC = MS._get_matrix_class()
+        MC = MS._matrix_class
         # FIXME: Hack because there is no ZZ \cup \{ \infty \}: -1 represents \infty
         E = UniversalCyclotomicField().gen
         if base_ring is UniversalCyclotomicField():

src/sage/matrix/action.pyx

@@ -117,42 +117,84 @@ cdef class MatrixMulAction(Action):
 
 
 cdef class MatrixMatrixAction(MatrixMulAction):
+    """
+    Action of a matrix on another matrix.
+
+    EXAMPLES:
+
+    By :trac:`715`, there only is a weak reference on the underlying set,
+    so that it can be garbage collected if only the action itself is
+    explicitly referred to. Hence, we first assign the involved matrix
+    spaces to a variable::
+
+        sage: R.<x> = ZZ[]
+        sage: MSR = MatrixSpace(R, 3, 3)
+        sage: MSQ = MatrixSpace(QQ, 3, 2)
+        sage: from sage.matrix.action import MatrixMatrixAction
+        sage: A = MatrixMatrixAction(MSR, MSQ); A
+        Left action by Full MatrixSpace of 3 by 3 dense matrices over Univariate Polynomial Ring in x over Integer Ring on Full MatrixSpace of 3 by 2 dense matrices over Rational Field
+        sage: A.codomain()
+        Full MatrixSpace of 3 by 2 dense matrices over Univariate Polynomial Ring in x over Rational Field
+        sage: A(matrix(R, 3, 3, x), matrix(QQ, 3, 2, range(6)))
+        [  0   x]
+        [2*x 3*x]
+        [4*x 5*x]
+
+    .. NOTE::
+
+        The :func:`MatrixSpace` function caches the object it creates.
+        Therefore, the underlying set ``MSZ`` in the above example will not
+        be garbage collected, even if it is not strongly ref'ed.
+        Nonetheless, there is no guarantee that the set that is acted upon
+        will always be cached in such a way, so that following the above
+        example is good practice.
+    """
     def __init__(self, G, S):
         """
-        EXAMPLES:
-
-        By :trac:`715`, there only is a weak reference on the underlying set,
-        so that it can be garbage collected if only the action itself is
-        explicitly referred to. Hence, we first assign the involved matrix
-        spaces to a variable::
-
-            sage: R.<x> = ZZ[]
-            sage: MSR = MatrixSpace(R, 3, 3)
-            sage: MSQ = MatrixSpace(QQ, 3, 2)
-            sage: from sage.matrix.action import MatrixMatrixAction
-            sage: A = MatrixMatrixAction(MSR, MSQ); A
-            Left action by Full MatrixSpace of 3 by 3 dense matrices over Univariate Polynomial Ring in x over Integer Ring on Full MatrixSpace of 3 by 2 dense matrices over Rational Field
-            sage: A.codomain()
-            Full MatrixSpace of 3 by 2 dense matrices over Univariate Polynomial Ring in x over Rational Field
-            sage: A(matrix(R, 3, 3, x), matrix(QQ, 3, 2, range(6)))
-            [  0   x]
-            [2*x 3*x]
-            [4*x 5*x]
-
-        .. NOTE::
-
-            The :func:`MatrixSpace` function caches the object it creates.
-            Therefore, the underlying set ``MSZ`` in the above example will not
-            be garbage collected, even if it is not strongly ref'ed.
-            Nonetheless, there is no guarantee that the set that is acted upon
-            will always be cached in such a way, so that following the above
-            example is good practice.
-
+        TESTS:
+
+        Check that multiplication for matrices with different backends are not allowed::
+
+            sage: M1 = MatrixSpace(ZZ, 2, implementation='flint')
+            sage: M2 = MatrixSpace(ZZ, 2, implementation='generic')
+            sage: M3 = MatrixSpace(ZZ, 2, sparse=True)
+            sage: M = [M1, M2, M3]
+
+            sage: coercions = ''
+            sage: for i in range(3):
+            ....:     for j in range(3):
+            ....:         try:
+            ....:             s = M[i].an_element() * M[j].an_element()
+            ....:             coercions += 'X'
+            ....:         except TypeError:
+            ....:             coercions += ' '
+            ....:     coercions += '\n'
+            sage: print(coercions)
+            X X
+             X
+            X X
         """
         if not is_MatrixSpace(S):
             raise TypeError("Not a matrix space: %s" % S)
+
         MatrixMulAction.__init__(self, G, S, True)
 
+        # disallow multiplication on different backends (same size and rings)
+        if G.base_ring() is S.base_ring() and \
+           G.is_sparse() == S.is_sparse() and \
+           G._matrix_class != S._matrix_class:
+            raise TypeError("no matrix multiplication between different implementations")
+
+        # disallow multiplication (sparse) x (dense) when the densification is not the default
+        # implementation
+        if self.fix_sparseness:
+            if G.is_sparse():
+                if not S._has_default_implementation():
+                    raise TypeError("matrix multiplication not allowed")
+            else:
+                if not G._has_default_implementation():
+                    raise TypeError("matrix multiplication not allowed")
+
     def _create_codomain(self, base):
         """
         EXAMPLES:
@@ -252,6 +294,7 @@ cdef class MatrixMatrixAction(MatrixMulAction):
                 B = B.dense_matrix()
             else:
                 A = A.dense_matrix()
+        assert type(A) == type(B), (type(A), type(B))
         prod = A._matrix_times_matrix_(B)
         if A._subdivisions is not None or B._subdivisions is not None:
             Asubs = A.subdivisions()

src/sage/matrix/matrix0.pyx

@@ -869,6 +869,19 @@ cdef class Matrix(sage.structure.element.Matrix):
             ...
             IndexError: row indices must be integers
 
+        Check that submatrices with a specified implementation have the
+        same implementation::
+
+            sage: M = MatrixSpace(GF(2), 3, 3, implementation='generic')
+            sage: m = M(range(9))
+            sage: type(m)
+            <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
+            sage: parent(m)
+            Full MatrixSpace of 3 by 3 dense matrices over Finite Field of size 2 (using Matrix_generic_dense)
+            sage: type(m[:2,:2])
+            <type 'sage.matrix.matrix_generic_dense.Matrix_generic_dense'>
+            sage: parent(m[:2,:2])
+            Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 2 (using Matrix_generic_dense)
         """
         cdef list row_list
         cdef list col_list