diff --git a/Project.toml b/Project.toml
index 0d75506b..5ad72414 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-Aqua = "0.5"
+Aqua = "0.5, 0.6"
 julia = "1.6"
 
 [extras]
diff --git a/README.md b/README.md
index 997b3fec..86d7137e 100644
--- a/README.md
+++ b/README.md
@@ -14,12 +14,8 @@ as well as identity matrices.  This package exports the following types:
 
 
 The primary purpose of this package is to present a unified way of constructing
-matrices. For example, to construct a 5-by-5 `CLArray` of all zeros, one would use
-```julia
-julia> CLArray(Zeros(5,5))
-```
-Because `Zeros` is lazy, this can be accomplished on the GPU with no memory transfer.
-Similarly, to construct a 5-by-5 `BandedMatrix` of all zeros with bandwidths `(1,2)`, one would use  
+matrices. 
+For example, to construct a 5-by-5 `BandedMatrix` of all zeros with bandwidths `(1,2)`, one would use  
 ```julia
 julia> BandedMatrix(Zeros(5,5), (1, 2))
 ```
diff --git a/src/FillArrays.jl b/src/FillArrays.jl
index bf090490..da758c54 100644
--- a/src/FillArrays.jl
+++ b/src/FillArrays.jl
@@ -6,7 +6,7 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
     +, -, *, /, \, diff, sum, cumsum, maximum, minimum, sort, sort!,
     any, all, axes, isone, iterate, unique, allunique, permutedims, inv,
     copy, vec, setindex!, count, ==, reshape, _throw_dmrs, map, zero,
-    show, view, in, mapreduce, one, reverse, promote_op
+    show, view, in, mapreduce, one, reverse, promote_op, promote_rule
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
@@ -18,7 +18,7 @@ import Base.Broadcast: broadcasted, DefaultArrayStyle, broadcast_shape
 import Statistics: mean, std, var, cov, cor
 
 
-export Zeros, Ones, Fill, Eye, Trues, Falses
+export Zeros, Ones, Fill, Eye, Trues, Falses, OneElement
 
 import Base: oneto
 
@@ -34,6 +34,7 @@ const AbstractFillVecOrMat{T} = Union{AbstractFillVector{T},AbstractFillMatrix{T
 
 ==(a::AbstractFill, b::AbstractFill) = axes(a) == axes(b) && getindex_value(a) == getindex_value(b)
 
+
 @inline function _fill_getindex(F::AbstractFill, kj::Integer...)
     @boundscheck checkbounds(F, kj...)
     getindex_value(F)
@@ -147,6 +148,27 @@ Fill{T,0}(x::T, ::Tuple{}) where T = Fill{T,0,Tuple{}}(x, ()) # ambiguity fix
 @inline axes(F::Fill) = F.axes
 @inline size(F::Fill) = map(length, F.axes)
 
+"""
+    getindex_value(F::AbstractFill)
+
+Return the value that `F` is filled with.
+
+# Examples
+
+```jldoctest
+julia> f = Ones(3);
+
+julia> FillArrays.getindex_value(f)
+1.0
+
+julia> g = Fill(2, 10);
+
+julia> FillArrays.getindex_value(g)
+2
+```
+"""
+getindex_value
+
 @inline getindex_value(F::Fill) = F.value
 
 AbstractArray{T}(F::Fill{V,N}) where {T,V,N} = Fill{T}(convert(T, F.value)::T, F.axes)
@@ -155,7 +177,12 @@ AbstractFill{T}(F::AbstractFill) where T = AbstractArray{T}(F)
 
 copy(F::Fill) = Fill(F.value, F.axes)
 
-""" Throws an error if `arr` does not contain one and only one unique value. """
+"""
+    unique_value(arr::AbstractArray)
+
+Return `only(unique(arr))` without intermediate allocations.
+Throws an error if `arr` does not contain one and only one unique value.
+"""
 function unique_value(arr::AbstractArray)
     if isempty(arr) error("Cannot convert empty array to Fill") end
     val = first(arr)
@@ -262,6 +289,7 @@ for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
         @inline $Typ{T,N}(A::AbstractArray{V,N}) where{T,V,N} = $Typ{T,N}(size(A))
         @inline $Typ{T}(A::AbstractArray) where{T} = $Typ{T}(size(A))
         @inline $Typ(A::AbstractArray) = $Typ{eltype(A)}(A)
+        @inline $Typ(::Type{T}, m...) where T = $Typ{T}(m...)
 
         @inline axes(Z::$Typ) = Z.axes
         @inline size(Z::$Typ) = length.(Z.axes)
@@ -273,6 +301,14 @@ for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
         copy(F::$Typ) = F
 
         getindex(F::$Typ{T,0}) where T = getindex_value(F)
+
+        promote_rule(::Type{$Typ{T, N, Axes}}, ::Type{$Typ{V, N, Axes}}) where {T,V,N,Axes} = $Typ{promote_type(T,V),N,Axes}
+        function convert(::Type{$Typ{T,N,Axes}}, A::$Typ{V,N,Axes}) where {T,V,N,Axes}
+            convert(T, getindex_value(A)) # checks that the types are convertible
+            $Typ{T,N,Axes}(axes(A))
+        end
+        convert(::Type{$Typ{T,N}}, A::$Typ{V,N,Axes}) where {T,V,N,Axes} = convert($Typ{T,N,Axes}, A)
+        convert(::Type{$Typ{T}}, A::$Typ{V,N,Axes}) where {T,V,N,Axes} = convert($Typ{T,N,Axes}, A)
     end
 end
 
@@ -284,6 +320,8 @@ for TYPE in (:Fill, :AbstractFill, :Ones, :Zeros), STYPE in (:AbstractArray, :Ab
     end
 end
 
+promote_rule(::Type{<:AbstractFill{T, N, Axes}}, ::Type{<:AbstractFill{V, N, Axes}}) where {T,V,N,Axes} = Fill{promote_type(T,V),N,Axes}
+
 """
     fillsimilar(a::AbstractFill, axes)
 
@@ -426,7 +464,8 @@ end
 
 
 ## Array
-Base.Array{T,N}(F::AbstractFill{V,N}) where {T,V,N} = fill(convert(T, getindex_value(F)), size(F))
+Base.Array{T,N}(F::AbstractFill{V,N}) where {T,V,N} =
+    convert(Array{T,N}, fill(convert(T, getindex_value(F)), size(F)))
 
 # These are in case `zeros` or `ones` are ever faster than `fill`
 for (Typ, funcs, func) in ((:Zeros, :zeros, :zero), (:Ones, :ones, :one))
@@ -437,7 +476,7 @@ end
 
 # temporary patch. should be a PR(#48895) to LinearAlgebra
 Diagonal{T}(A::AbstractFillMatrix) where T = Diagonal{T}(diag(A))
-function convert(::Type{T}, A::AbstractFillMatrix) where T<:Diagonal 
+function convert(::Type{T}, A::AbstractFillMatrix) where T<:Diagonal
     checksquare(A)
     isdiag(A) ? T(A) : throw(InexactError(:convert, T, A))
 end
@@ -496,14 +535,14 @@ sum(x::AbstractFill) = getindex_value(x)*length(x)
 sum(f, x::AbstractFill) = length(x) * f(getindex_value(x))
 sum(x::Zeros) = getindex_value(x)
 
-cumsum(x::AbstractFill{<:Any,1}) = range(getindex_value(x); step=getindex_value(x),
-                                                    length=length(x))
+# needed to support infinite case
+steprangelen(st...) = StepRangeLen(st...)
+cumsum(x::AbstractFill{<:Any,1}) = steprangelen(getindex_value(x), getindex_value(x), length(x))
 
 cumsum(x::ZerosVector) = x
 cumsum(x::ZerosVector{Bool}) = x
 cumsum(x::OnesVector{II}) where II<:Integer = convert(AbstractVector{II}, oneto(length(x)))
 cumsum(x::OnesVector{Bool}) = oneto(length(x))
-cumsum(x::AbstractFillVector{Bool}) = cumsum(AbstractFill{Int}(x))
 
 
 #########
@@ -717,4 +756,6 @@ Base.@propagate_inbounds function view(A::AbstractFill{<:Any,N}, I::Vararg{Real,
     fillsimilar(A)
 end
 
+include("oneelement.jl")
+
 end # module
diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl
index 2dec1b61..800e803a 100644
--- a/src/fillalgebra.jl
+++ b/src/fillalgebra.jl
@@ -86,7 +86,6 @@ end
 *(a::ZerosMatrix, b::AbstractMatrix) = mult_zeros(a, b)
 *(a::AbstractMatrix, b::ZerosVector) = mult_zeros(a, b)
 *(a::AbstractMatrix, b::ZerosMatrix) = mult_zeros(a, b)
-*(a::ZerosVector, b::AbstractVector) = mult_zeros(a, b)
 *(a::ZerosMatrix, b::AbstractVector) = mult_zeros(a, b)
 *(a::AbstractVector, b::ZerosMatrix) = mult_zeros(a, b)
 
diff --git a/src/fillbroadcast.jl b/src/fillbroadcast.jl
index 574fb145..82ba41e2 100644
--- a/src/fillbroadcast.jl
+++ b/src/fillbroadcast.jl
@@ -102,10 +102,6 @@ _broadcasted_zeros(f, a, b) = Zeros{Base.Broadcast.combine_eltypes(f, (a, b))}(b
 _broadcasted_ones(f, a, b) = Ones{Base.Broadcast.combine_eltypes(f, (a, b))}(broadcast_shape(axes(a), axes(b)))
 _broadcasted_nan(f, a, b) = Fill(convert(Base.Broadcast.combine_eltypes(f, (a, b)), NaN), broadcast_shape(axes(a), axes(b)))
 
-# TODO: remove at next breaking version
-_broadcasted_zeros(a, b) = _broadcasted_zeros(+, a, b)
-_broadcasted_ones(a, b) = _broadcasted_ones(+, a, b)
-
 broadcasted(::DefaultArrayStyle, ::typeof(+), a::Zeros, b::Zeros) = _broadcasted_zeros(+, a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(+), a::Ones, b::Zeros) = _broadcasted_ones(+, a, b)
 broadcasted(::DefaultArrayStyle, ::typeof(+), a::Zeros, b::Ones) = _broadcasted_ones(+, a, b)
@@ -247,3 +243,8 @@ broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{
 broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, r::Ones{T,N}, ::Base.RefValue{Val{k}}) where {T,N,k} = Ones{T}(axes(r))
 broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, r::Zeros{T,N}, ::Base.RefValue{Val{0}}) where {T,N} = Ones{T}(axes(r))
 broadcasted(::DefaultArrayStyle{N}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, r::Zeros{T,N}, ::Base.RefValue{Val{k}}) where {T,N,k} = Zeros{T}(axes(r))
+
+# supports structured broadcast
+if isdefined(LinearAlgebra, :fzero)
+    LinearAlgebra.fzero(x::Zeros) = zero(eltype(x))
+end
\ No newline at end of file
diff --git a/src/oneelement.jl b/src/oneelement.jl
new file mode 100644
index 00000000..a51a3077
--- /dev/null
+++ b/src/oneelement.jl
@@ -0,0 +1,58 @@
+"""
+    OneElement(val, ind, axesorsize) <: AbstractArray
+
+Represents an array with the specified axes (if its a tuple of `AbstractUnitRange`s)
+or size (if its a tuple of `Integer`s), with a single entry set to `val` and all others equal to zero,
+specified by `ind``.
+"""
+struct OneElement{T,N,I,A} <: AbstractArray{T,N}
+  val::T
+  ind::I
+  axes::A
+  OneElement(val::T, ind::I, axes::A) where {T<:Number, I<:NTuple{N,Int}, A<:NTuple{N,AbstractUnitRange}} where {N} = new{T,N,I,A}(val, ind, axes)
+end
+
+OneElement(val, inds::NTuple{N,Int}, sz::NTuple{N,Integer}) where N = OneElement(val, inds, oneto.(sz))
+"""
+    OneElement(val, ind::Int, n::Int)
+
+Creates a length `n` vector where the `ind` entry is equal to `val`, and all other entries are zero.
+"""
+OneElement(val, ind::Int, len::Int) = OneElement(val, (ind,), (len,))
+"""
+    OneElement(ind::Int, n::Int)
+
+Creates a length `n` vector where the `ind` entry is equal to `1`, and all other entries are zero.
+"""
+OneElement(inds::Int, sz::Int) = OneElement(1, inds, sz)
+OneElement{T}(val, inds::NTuple{N,Int}, sz::NTuple{N,Integer}) where {T,N} = OneElement(convert(T,val), inds, oneto.(sz))
+OneElement{T}(val, inds::Int, sz::Int) where T = OneElement{T}(val, (inds,), (sz,))
+
+"""
+    OneElement{T}(val, ind::Int, n::Int)
+
+Creates a length `n` vector where the `ind` entry is equal to `one(T)`, and all other entries are zero.
+"""
+OneElement{T}(inds::Int, sz::Int) where T = OneElement(one(T), inds, sz)
+
+Base.size(A::OneElement) = map(length, A.axes)
+Base.axes(A::OneElement) = A.axes
+function Base.getindex(A::OneElement{T,N}, kj::Vararg{Int,N}) where {T,N}
+    @boundscheck checkbounds(A, kj...)
+    ifelse(kj == A.ind, A.val, zero(T))
+end
+
+Base.replace_in_print_matrix(o::OneElement{<:Any,2}, k::Integer, j::Integer, s::AbstractString) =
+    o.ind == (k,j) ? s : Base.replace_with_centered_mark(s)
+
+function Base.setindex(A::Zeros{T,N}, v, kj::Vararg{Int,N}) where {T,N}
+    @boundscheck checkbounds(A, kj...)
+    OneElement(convert(T, v), kj, axes(A))
+end
+
+
+Base.@propagate_inbounds function view(A::RectOrDiagonal{<:Any,<:AbstractFill}, kr::AbstractRange, j::Integer)
+    @boundscheck checkbounds(A, kr, j)
+    k = findfirst(isequal(j), kr)
+    OneElement(getindex_value(A.diag), isnothing(k) ? 0 : something(k), length(kr))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 81ba61c5..6df2edb5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,6 +20,7 @@ include("infinitearrays.jl")
 
             for T in (Int, Float64)
                 Z = $Typ{T}(5)
+                @test $Typ(T, 5) ≡ Z
                 @test eltype(Z) == T
                 @test Array(Z) == $funcs(T,5)
                 @test Array{T}(Z) == $funcs(T,5)
@@ -28,12 +29,14 @@ include("infinitearrays.jl")
                 @test convert(AbstractArray,Z) ≡ Z
                 @test convert(AbstractArray{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                 @test convert(AbstractVector{T},Z) ≡ AbstractVector{T}(Z) ≡ Z
+                @test convert(AbstractFill{T},Z) ≡ AbstractFill{T}(Z) ≡ Z
 
                 @test $Typ{T,1}(2ones(T,5)) == Z
                 @test $Typ{T}(2ones(T,5)) == Z
                 @test $Typ(2ones(T,5)) == Z
 
                 Z = $Typ{T}(5, 5)
+                @test $Typ(T, 5, 5) ≡ Z
                 @test eltype(Z) == T
                 @test Array(Z) == $funcs(T,5,5)
                 @test Array{T}(Z) == $funcs(T,5,5)
@@ -194,6 +197,23 @@ include("infinitearrays.jl")
         A = FillArrays.RectDiagonal(Int[], (1:0, 1:4))
         @test !(0 in A)
     end
+
+    @testset "promotion" begin
+        Z = Zeros{Int}(5)
+        Zf = Zeros(5)
+        O = Ones{Int}(5)
+        Of = Ones{Float64}(5)
+        @test [Z,O] isa Vector{Fill{Int,1,Tuple{Base.OneTo{Int}}}}
+        @test [Z,Of] isa Vector{Fill{Float64,1,Tuple{Base.OneTo{Int}}}}
+        @test [O,O] isa Vector{Ones{Int,1,Tuple{Base.OneTo{Int}}}}
+        @test [O,Of] isa Vector{Ones{Float64,1,Tuple{Base.OneTo{Int}}}}
+        @test [Z,Zf] isa Vector{Zeros{Float64,1,Tuple{Base.OneTo{Int}}}}
+
+        @test convert(Ones{Int}, Of) ≡ convert(Ones{Int,1}, Of) ≡ convert(typeof(O), Of) ≡ O
+        @test convert(Zeros{Int}, Zf) ≡ convert(Zeros{Int,1}, Zf) ≡ convert(typeof(Z), Zf) ≡ Z
+
+        @test_throws MethodError convert(Zeros{SVector{2,Int}}, Zf)
+    end
 end
 
 @testset "indexing" begin
@@ -508,9 +528,9 @@ end
     @test_throws MethodError [1,2,3]*Zeros(3) # Not defined for [1,2,3]*[0,0,0] either
 
     @testset "Check multiplication by Adjoint vectors works as expected." begin
-        @test randn(4, 3)' * Zeros(4) === Zeros(3)
-        @test randn(4)' * Zeros(4) === zero(Float64)
-        @test [1, 2, 3]' * Zeros{Int}(3) === zero(Int)
+        @test randn(4, 3)' * Zeros(4) ≡ Zeros(3)
+        @test randn(4)' * Zeros(4) ≡ transpose(randn(4)) * Zeros(4) ≡ zero(Float64)
+        @test [1, 2, 3]' * Zeros{Int}(3) ≡ zero(Int)
         @test [SVector(1,2)', SVector(2,3)', SVector(3,4)']' * Zeros{Int}(3) === SVector(0,0)
         @test_throws DimensionMismatch randn(4)' * Zeros(3)
         @test Zeros(5)' * randn(5,3) ≡ Zeros(5)'*Zeros(5,3) ≡ Zeros(5)'*Ones(5,3) ≡ Zeros(3)'
@@ -533,6 +553,7 @@ end
         @test randn(5) * transpose(Zeros(6)) ≡ randn(5,1) * transpose(Zeros(6)) ≡ Zeros(5,6)
         @test Zeros(5) * transpose(randn(6)) ≡ Zeros(5,6)
         @test transpose(randn(5)) * Zeros(5) ≡ 0.0
+        @test transpose(randn(5) .+ im) * Zeros(5) ≡ 0.0 + 0im
 
         @test transpose([[1,2]]) * Zeros{SVector{2,Int}}(1) ≡ 0
         @test_broken transpose([[1,2,3]]) * Zeros{SVector{2,Int}}(1)
@@ -600,12 +621,13 @@ end
 
     @testset "Tests for ranges." begin
         X = randn(5)
-        @test !(Zeros(5) + X === X)
-        @test Zeros{Int}(5) + (1:5) === (1:5) && (1:5) + Zeros{Int}(5) === (1:5)
-        @test Zeros(5) + (1:5) === (1.0:1.0:5.0) && (1:5) + Zeros(5) === (1.0:1.0:5.0)
-        @test (1:5) - Zeros{Int}(5) === (1:5)
-        @test Zeros{Int}(5) - (1:5) === -1:-1:-5
-        @test Zeros(5) - (1:5) === -1.0:-1.0:-5.0
+        @test !(Zeros(5) + X ≡ X)
+        @test Zeros{Int}(5) + (1:5) ≡ (1:5) + Zeros{Int}(5) ≡ (1:5)
+        @test Zeros(5) + (1:5) ≡ (1:5) + Zeros(5) ≡ (1.0:1.0:5.0)
+        @test (1:5) - Zeros{Int}(5) ≡ (1:5)
+        @test Zeros{Int}(5) - (1:5) ≡ -1:-1:-5
+        @test Zeros(5) - (1:5) ≡ -1.0:-1.0:-5.0
+        @test Zeros{Int}(5) + (1.0:5) ≡ (1.0:5) + Zeros{Int}(5) ≡ 1.0:5
     end
 
     @testset "test Base.zero" begin
@@ -634,11 +656,11 @@ end
     @test sum(Fill(3,10)) ≡ 30
     @test reduce(+, Fill(3,10)) ≡ 30
     @test sum(x -> x + 1, Fill(3,10)) ≡ 40
-    @test cumsum(Fill(3,10)) ≡ 3:3:30
+    @test cumsum(Fill(3,10)) ≡ StepRangeLen(3,3,10)
 
     @test sum(Ones(10)) ≡ 10.0
     @test sum(x -> x + 1, Ones(10)) ≡ 20.0
-    @test cumsum(Ones(10)) ≡ 1.0:10.0
+    @test cumsum(Ones(10)) ≡ StepRangeLen(1.0, 1.0, 10)
 
     @test sum(Ones{Int}(10)) ≡ 10
     @test sum(x -> x + 1, Ones{Int}(10)) ≡ 20
@@ -654,7 +676,7 @@ end
 
     @test cumsum(Zeros{Bool}(10)) ≡ Zeros{Bool}(10)
     @test cumsum(Ones{Bool}(10)) ≡ Base.OneTo{Int}(10)
-    @test cumsum(Fill(true,10)) ≡ 1:1:10
+    @test cumsum(Fill(true,10)) ≡ StepRangeLen(true, true, 10)
 
     @test diff(Fill(1,10)) ≡ Zeros{Int}(9)
     @test diff(Ones{Float64}(10)) ≡ Zeros{Float64}(9)
@@ -759,6 +781,8 @@ end
         @test broadcast(*, rnge, Zeros(10, 10)) ≡ Zeros{Float64}(10, 10)
         @test broadcast(*, Ones{Int}(10), rnge) ≡ rnge
         @test broadcast(*, rnge, Ones{Int}(10)) ≡ rnge
+        @test broadcast(*, Ones(10), -5:4) ≡ broadcast(*, -5:4, Ones(10)) ≡ rnge
+        @test broadcast(*, Ones(10), -5:1:4) ≡ broadcast(*, -5:1:4, Ones(10)) ≡ rnge
         @test_throws DimensionMismatch broadcast(*, Fill(5.0, 11), rnge)
         @test broadcast(*, rnge, Fill(5.0, 10)) == broadcast(*, rnge, 5.0)
         @test_throws DimensionMismatch broadcast(*, rnge, Fill(5.0, 11))
@@ -1338,6 +1362,8 @@ end
     @test stringmime("text/plain", Fill(7,2,3)) == "2×3 Fill{$Int}, with entries equal to 7"
     @test stringmime("text/plain", Fill(8.0,1)) == "1-element Fill{Float64}, with entry equal to 8.0"
     @test stringmime("text/plain", Eye(5)) == "5×5 Eye{Float64}"
+    # used downstream in LazyArrays.jl to deduce sparsity
+    @test Base.replace_in_print_matrix(Zeros(5,3), 1, 2, "0.0") == " ⋅ "
 
     # 2-arg show, compact printing
     @test repr(Zeros(3)) == "Zeros(3)"
@@ -1469,3 +1495,40 @@ end
     @test cor(Fill(3, 4, 5)) ≈ cor(fill(3, 4, 5)) nans=true
     @test cor(Fill(3, 4, 5), dims=2) ≈ cor(fill(3, 4, 5), dims=2) nans=true
 end
+
+@testset "Structured broadcast" begin
+    D = Diagonal(1:5)
+    @test D + Zeros(5,5) isa Diagonal
+    @test D - Zeros(5,5) isa Diagonal
+    @test D .+ Zeros(5,5) isa Diagonal
+    @test D .- Zeros(5,5) isa Diagonal
+    @test D .* Zeros(5,5) isa Diagonal
+    @test Zeros(5,5) .* D isa Diagonal
+    @test Zeros(5,5) - D isa Diagonal
+    @test Zeros(5,5) + D isa Diagonal
+    @test Zeros(5,5) .- D isa Diagonal
+    @test Zeros(5,5) .+ D isa Diagonal
+    f = (x,y) -> x+1
+    @test f.(D, Zeros(5,5)) isa Matrix
+end
+
+@testset "OneElement" begin
+    e₁ = OneElement(2, 5)
+    @test e₁ == [0,1,0,0,0]
+    @test_throws BoundsError e₁[6]
+
+    e₁ = OneElement{Float64}(2, 5)
+    @test e₁ == [0,1,0,0,0]
+
+    v = OneElement{Float64}(2, 3, 4)
+    @test v == [0,0,2,0]
+
+    V = OneElement(2, (2,3), (3,4))
+    @test V == [0 0 0 0; 0 0 2 0; 0 0 0 0]
+
+    @test stringmime("text/plain", V) == "3×4 OneElement{$Int, 2, Tuple{$Int, $Int}, Tuple{Base.OneTo{$Int}, Base.OneTo{$Int}}}:\n ⋅  ⋅  ⋅  ⋅\n ⋅  ⋅  2  ⋅\n ⋅  ⋅  ⋅  ⋅"
+
+    @test Base.setindex(Zeros(5), 2, 2) ≡ OneElement(2.0, 2, 5)
+    @test Base.setindex(Zeros(5,3), 2, 2, 3) ≡ OneElement(2.0, (2,3), (5,3))
+    @test_throws BoundsError Base.setindex(Zeros(5), 2, 6)
+end
\ No newline at end of file