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

@@ -1038,17 +1038,22 @@ const NCRingElement = Union{NCRingElem, RingElement}
 abstract type Mat{T} <: MatElem{T} end
 
 # not really a mathematical ring
-struct MatSpace{T <: NCRingElement} <: AbstractAlgebra.MatSpace{T}
+struct MatSpace{CoeffT, RingT} <: AbstractAlgebra.MatSpace{CoeffT}
+   base_ring::RingT
    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)
+   function MatSpace{CoeffT, RingT}(R::RingT, r::Int, c::Int) where {CoeffT, RingT <: NCRing}
+      @assert elem_type(RingT) === CoeffT
+      (r < 0 || c < 0) && throw(error_dim_negative)
+      return new{CoeffT, RingT}(R, r, c)
    end
 end
 
+function MatSpace(R::T, r::Int, c::Int) where {T <: NCRing}
+   return MatSpace{elem_type(T), T}(R, r, c)
+end
+
 struct MatSpaceElem{T <: NCRingElement} <: Mat{T}
    base_ring::NCRing
    entries::Matrix{T}
@@ -1073,7 +1078,10 @@ end
 
 # construct zero matrix
 function MatSpaceElem{T}(R::NCRing, r::Int, c::Int) where T <: NCRingElement
-   entries = fill(zero(R), r, c)
+   entries = Matrix{T}(undef, r, c)
+   for i in 1:r*c
+      entries[i] = zero(R)
+   end
    return MatSpaceElem{T}(R, entries)
 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, parent_type(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
@@ -43,6 +43,8 @@ dense_matrix_type(::Type{T}) where T <: NCRingElement = MatSpaceElem{T}
 #
 ###############################################################################
 
+base_ring(a::MatSpace) = a.base_ring
+
 @doc raw"""
     nrows(a::MatSpace)
 
@@ -157,63 +159,44 @@ end
 #
 ###############################################################################
 
+# create a zero matrix
 function (a::MatSpace{T})() where {T <: NCRingElement}
-   R = base_ring(a)
-   entries = Array{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 dense_matrix_type(T)(base_ring(a), nrows(a), ncols(a))
 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 = Array{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 = Array{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)))  # zero matrix
 end
 
 ###############################################################################
@@ -222,9 +205,6 @@ end
 #
 ###############################################################################
 
-function matrix_space(R::AbstractAlgebra.NCRing, r::Int, c::Int; cached::Bool = true)
-   # TODO: the 'cached' argument is ignored and mainly here for backwards compatibility
-   # (and perhaps future compatibility, in case we need it again)
-   T = elem_type(R)
-   return MatSpace{T}(R, r, c)
+function matrix_space(R::AbstractAlgebra.NCRing, r::Int, c::Int)
+   return MatSpace(R, r, c)
 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, F2}
 
 (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.nrows(a::F2Matrix) = nrows(a.m)
 AbstractAlgebra.ncols(a::F2Matrix) = ncols(a.m)
@@ -121,7 +117,7 @@ end
 
    @test elem_type(S) == Generic.MatSpaceElem{elem_type(R)}
    @test elem_type(Generic.MatSpace{elem_type(R)}) == Generic.MatSpaceElem{elem_type(R)}
-   @test parent_type(Generic.MatSpaceElem{elem_type(R)}) == Generic.MatSpace{elem_type(R)}
+   @test parent_type(Generic.MatSpaceElem{elem_type(R)}) == Generic.MatSpace{elem_type(R), typeof(R)}
    @test base_ring(S) == R
    @test nrows(S) == 3
    @test ncols(S) == 3
@@ -323,7 +319,7 @@ end
 
    @test elem_type(S) == Generic.MatSpaceElem{elem_type(R)}
    @test elem_type(Generic.MatSpace{elem_type(R)}) == Generic.MatSpaceElem{elem_type(R)}
-   @test parent_type(Generic.MatSpaceElem{elem_type(R)}) == Generic.MatSpace{elem_type(R)}
+   @test parent_type(Generic.MatSpaceElem{elem_type(R)}) == Generic.MatSpace{elem_type(R), typeof(R)}
 
    @test dense_matrix_type(R) == elem_type(S)
    @test dense_matrix_type(R(1)) == elem_type(S)