Make Generic.MatSpace truly generic

docs/src/matrix_interface.md

@@ -44,7 +44,8 @@ It also provides two abstract types for matrix algebras and their elements:
 Note that these abstract types are parameterised. The type `T` should usually be the
 type of elements of the matrices.
 
-Matrix spaces and matrix algebras should be made unique on the system by caching parent
+Matrix spaces and matrix algebras should be made unique on the system by either making
+them `struct` types, or by caching parent
 objects (unless an optional `cache` parameter is set to `false`). Matrix spaces and
 algebras should at least be distinguished based on their base (coefficient) ring and the
 dimensions of the matrices in the space.

src/generic/GenericTypes.jl

@@ -1054,13 +1054,15 @@ abstract type Mat{T} <: MatElem{T} end
 
 # not really a mathematical ring
 struct MatSpace{T <: NCRingElement} <: AbstractAlgebra.MatSpace{T}
+   base_ring::NCRing
    nrows::Int
    ncols::Int
-   base_ring::NCRing
 
    function MatSpace{T}(R::NCRing, r::Int, c::Int, cached::Bool = true) where T <: NCRingElement
-       # TODO/FIXME: `cached` is ignored and only exists for backwards compatibility
-       new{T}(r, c, R)
+      # TODO/FIXME: `cached` is ignored and only exists for backwards compatibility
+      @assert elem_type(R) === T
+      (r < 0 || c < 0) && error("Dimensions must be non-negative")
+      return new{T}(R, r, c)
    end
 end
 

src/generic/Matrix.jl

@@ -10,16 +10,16 @@
 #
 ###############################################################################
 
-parent_type(::Type{S}) where {T <: NCRingElement, S <: Mat{T}} = MatSpace{T}
+parent_type(::Type{<:MatElem{T}}) where {T <: NCRingElement} = MatSpace{T}
 
-elem_type(::Type{MatSpace{T}}) where {T <: NCRingElement} = MatSpaceElem{T}
+elem_type(::Type{MatSpace{T}}) where {T <: NCRingElement} = dense_matrix_type(T)
 
 @doc raw"""
-    parent(a::AbstractAlgebra.MatElem{T}) where T <: NCRingElement
+    parent(a::AbstractAlgebra.MatElem)
 
 Return the parent object of the given matrix.
 """
-parent(a::Mat{T}) where T <: NCRingElement = MatSpace{T}(a.base_ring, size(a.entries)...)
+parent(a::MatElem) = matrix_space(base_ring(a), nrows(a), ncols(a))
 
 @doc raw"""
     dense_matrix_type(::Type{T}) where T<:NCRingElement
@@ -148,63 +148,44 @@ end
 #
 ###############################################################################
 
+# create a zero matrix
 function (a::MatSpace{T})() where {T <: NCRingElement}
-   R = base_ring(a)
-   entries = Matrix{T}(undef, a.nrows, a.ncols)
-   for i = 1:a.nrows
-      for j = 1:a.ncols
-         entries[i, j] = zero(R)
-      end
-   end
-   z = MatSpaceElem{T}(R, entries)
-   return z
+   return zero_matrix(base_ring(a), nrows(a), ncols(a))::dense_matrix_type(T)
 end
 
-function (a::MatSpace{T})(b::S) where {S <: NCRingElement, T <: NCRingElement}
+# create a matrix with b on the diagonal
+function (a::AbstractAlgebra.Generic.MatSpace)(b::NCRingElement)
+   M = a()  # zero matrix
    R = base_ring(a)
-   entries = Matrix{T}(undef, a.nrows, a.ncols)
    rb = R(b)
-   for i = 1:a.nrows
-      for j = 1:a.ncols
-         if i != j
-            entries[i, j] = zero(R)
-         else
-            entries[i, j] = rb
-         end
-      end
+   for i in 1:min(nrows(a), ncols(a))
+      M[i, i] = rb
    end
-   z = MatSpaceElem{T}(R, entries)
-   return z
+   return M
 end
 
-function (a::MatSpace{T})(b::Matrix{T}) where T <: NCRingElement
+# convert a Julia matrix
+function (a::MatSpace{T})(b::AbstractMatrix{S}) where {T <: NCRingElement, S}
+   _check_dim(nrows(a), ncols(a), b)
    R = base_ring(a)
-   _check_dim(a.nrows, a.ncols, b)
-   if !isempty(b)
-      R != parent(b[1, 1]) && error("Unable to coerce matrix")
+
+   # minor optimization for MatSpaceElem
+   if S === T && dense_matrix_type(T) === MatSpaceElem{T} && all(x -> R === parent(x), b)
+      return MatSpaceElem{T}(R, b)
    end
-   z = MatSpaceElem{T}(R, b)
-   return z
-end
 
-function (a::MatSpace{T})(b::AbstractMatrix{S}) where {S <: NCRingElement, T <: NCRingElement}
-   R = base_ring(a)
-   _check_dim(a.nrows, a.ncols, b)
-   entries = Matrix{T}(undef, a.nrows, a.ncols)
-   for i = 1:a.nrows
-      for j = 1:a.ncols
-         entries[i, j] = R(b[i, j])
-      end
+   # generic code
+   M = a()  # zero matrix
+   for i = 1:nrows(a), j = 1:ncols(a)
+      M[i, j] = R(b[i, j])
    end
-   z = MatSpaceElem{T}(R, entries)
-   return z
+   return M
 end
 
-function (a::MatSpace{T})(b::AbstractVector{S}) where {S <: NCRingElement, T <: NCRingElement}
-   _check_dim(a.nrows, a.ncols, b)
-   b = Matrix{S}(transpose(reshape(b, a.ncols, a.nrows)))
-   z = a(b)
-   return z
+# convert a Julia vector
+function (a::MatSpace{T})(b::AbstractVector) where T <: NCRingElement
+   _check_dim(nrows(a), ncols(a), b)
+   return a(transpose(reshape(b, a.ncols, a.nrows)))
 end
 
 ###############################################################################

test/generic/Matrix-test.jl

@@ -65,10 +65,7 @@ AbstractAlgebra.divexact(x::F2Elem, y::F2Elem) = y.x ? x : throw(DivideError())
 Random.rand(rng::AbstractRNG, sp::Random.SamplerTrivial{F2}) = F2Elem(rand(rng, Bool))
 Random.gentype(::Type{F2}) = F2Elem
 
-struct F2MatSpace <: AbstractAlgebra.MatSpace{F2Elem}
-   nrows::Int
-   ncols::Int
-end
+const F2MatSpace = AbstractAlgebra.Generic.MatSpace{F2Elem}
 
 (S::F2MatSpace)() = zero_matrix(F2(), S.nrows, S.ncols)
 
@@ -81,8 +78,7 @@ AbstractAlgebra.parent_type(::Type{F2Matrix}) = F2MatSpace
 
 AbstractAlgebra.base_ring(::F2MatSpace) = F2()
 AbstractAlgebra.dense_matrix_type(::Type{F2}) = F2Matrix
-AbstractAlgebra.parent(a::F2Matrix) = F2MatSpace(nrows(a), ncols(a))
-AbstractAlgebra.matrix_space(::F2, r::Int, c::Int) = F2MatSpace(r, c)
+AbstractAlgebra.matrix_space(::F2, r::Int, c::Int) = F2MatSpace(F2(), r, c)
 
 AbstractAlgebra.number_of_rows(a::F2Matrix) = nrows(a.m)
 AbstractAlgebra.number_of_columns(a::F2Matrix) = ncols(a.m)