From 58215a890171d031a03dc4a2c5ce2aade958b264 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Wed, 3 Jan 2024 22:57:20 +0530
Subject: [PATCH 01/35] Broadcast arithmetic operators for OneElement (#331)

---
 src/oneelement.jl | 23 +++++++++++++++++++++++
 test/runtests.jl  | 14 ++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index 869f56c8..46b2b471 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -151,3 +151,26 @@ end
 
 adjoint(A::OneElementMatrix) = OneElement(adjoint(A.val), reverse(A.ind), reverse(A.axes))
 transpose(A::OneElementMatrix) = OneElement(transpose(A.val), reverse(A.ind), reverse(A.axes))
+
+# broadcast
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(conj), r::OneElement{<:Any,N}) where {N}
+    OneElement(conj(r.val), r.ind, axes(r))
+end
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(real), r::OneElement{<:Any,N}) where {N}
+    OneElement(real(r.val), r.ind, axes(r))
+end
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(imag), r::OneElement{<:Any,N}) where {N}
+    OneElement(imag(r.val), r.ind, axes(r))
+end
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(^), r::OneElement{<:Any,N}, x::Number) where {N}
+    OneElement(r.val^x, r.ind, axes(r))
+end
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(*), r::OneElement{<:Any,N}, x::Number) where {N}
+    OneElement(r.val*x, r.ind, axes(r))
+end
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(/), r::OneElement{<:Any,N}, x::Number) where {N}
+    OneElement(r.val/x, r.ind, axes(r))
+end
+function broadcasted(::DefaultArrayStyle{N}, ::typeof(\), x::Number, r::OneElement{<:Any,N}) where {N}
+    OneElement(x \ r.val, r.ind, axes(r))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0de16a98..293115d8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2218,6 +2218,20 @@ end
         @test A * (2 + 3.0im) === OneElement(val * (2 + 3.0im), (2,2), (3,4)) == B * (2 + 3.0im)
         @test A / (2 + 3.0im) === OneElement(val / (2 + 3.0im), (2,2), (3,4)) == B / (2 + 3.0im)
     end
+
+    @testset "broadcasting" begin
+        for v in (OneElement(2, 3, 4), OneElement(2im, (1,2), (3,4)))
+            w = Array(v)
+            n = 2
+            @test real.(v) == real.(w)
+            @test imag.(v) == imag.(w)
+            @test conj.(v) == conj.(w)
+            @test v .^ n == w .^ n
+            @test v .* n == w .* n
+            @test v ./ n == w ./ n
+            @test n .\ v == n .\ w
+        end
+    end
 end
 
 @testset "repeat" begin

From 617fffdf7dc561d1c1845ee84df0a510487a363f Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 9 Jan 2024 14:46:16 +0530
Subject: [PATCH 02/35] Disable downstream tests against broken packages (#338)

* Disable downstream tests against broken packages

LazyArrays and LazyBandedMatrices have known breakages on v1.10, so these tests need to be disabled until these are fixed.

* Disable InfiniteLinearAlgebra test

This also has a known breakage
---
 .github/workflows/downstream.yml | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index 598c4739..af31ed12 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -36,12 +36,12 @@ jobs:
         package:
           - {repo: Distributions.jl, group: JuliaStats}
           - {repo: BlockArrays.jl, group: JuliaArrays}
-          - {repo: LazyArrays.jl, group: JuliaArrays}
+          # - {repo: LazyArrays.jl, group: JuliaArrays}
           - {repo: ArrayLayouts.jl, group: JuliaLinearAlgebra}
-          - {repo: LazyBandedMatrices.jl, group: JuliaLinearAlgebra}
+          # - {repo: LazyBandedMatrices.jl, group: JuliaLinearAlgebra}
           - {repo: BandedMatrices.jl, group: JuliaLinearAlgebra}
           - {repo: BlockBandedMatrices.jl, group: JuliaLinearAlgebra}
-          - {repo: InfiniteLinearAlgebra.jl, group: JuliaLinearAlgebra}
+          # - {repo: InfiniteLinearAlgebra.jl, group: JuliaLinearAlgebra}
           - {repo: Optim.jl, group: JuliaNLSolvers}
 
     steps:

From 6b69b81961cc4e5499c6bb35aa3842935b73ed19 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 9 Jan 2024 15:18:25 +0530
Subject: [PATCH 03/35] Change downstream fail fast to false (#337)

---
 .github/workflows/downstream.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index af31ed12..b804e802 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -29,7 +29,7 @@ jobs:
     name: ${{ matrix.package.group }}/${{ matrix.package.repo }}/${{ matrix.julia-version }}
     runs-on: ${{ matrix.os }}
     strategy:
-      fail-fast: true
+      fail-fast: false
       matrix:
         julia-version: ['1']
         os: [ubuntu-latest]

From e9077c3197cf542fa76f51954a7fce1345ea85e2 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 9 Jan 2024 15:18:42 +0530
Subject: [PATCH 04/35] Remove explicit test on v1.10 (#336)

---
 .github/workflows/ci.yml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 72bf50f0..40668356 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -30,7 +30,6 @@ jobs:
         version:
           - '1.6'
           - '1'
-          - '~1.10.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

From 9f5c6aaf2e9b7db2eb12c93dd85b502f5a65fc1e Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 11 Jan 2024 20:47:03 +0530
Subject: [PATCH 05/35] Remove Random as a dependendency

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 06ba0f13..d980a8a5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "1.10.0"
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
@@ -42,6 +41,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Quaternions = "94ee1d12-ae83-5a48-8b1c-48b8ff168ae0"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -49,4 +49,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test", "Base64", "Infinities", "PDMats", "ReverseDiff", "SparseArrays", "StaticArrays", "Statistics", "Quaternions", "Documenter"]
+test = ["Aqua", "Test", "Base64", "Infinities", "PDMats", "ReverseDiff", "SparseArrays", "StaticArrays", "Statistics", "Quaternions", "Documenter", "Random"]

From 7f60221afb45c709ed19e1866a88ecc2caf4adab Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Fri, 12 Jan 2024 10:29:59 +0530
Subject: [PATCH 06/35] Move Aqua tests to a separate file (#339)

---
 test/aqua.jl     | 9 +++++++++
 test/runtests.jl | 8 +-------
 2 files changed, 10 insertions(+), 7 deletions(-)
 create mode 100644 test/aqua.jl

diff --git a/test/aqua.jl b/test/aqua.jl
new file mode 100644
index 00000000..3ad9d957
--- /dev/null
+++ b/test/aqua.jl
@@ -0,0 +1,9 @@
+using Aqua
+using FillArrays
+using Test
+@testset "Project quality" begin
+    Aqua.test_all(FillArrays;
+        # https://github.com/JuliaArrays/FillArrays.jl/issues/105#issuecomment-1582516319
+        ambiguities=(; broken=true),
+    )
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 293115d8..be4b8fec 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,13 +5,7 @@ using Documenter
 DocMeta.setdocmeta!(FillArrays, :DocTestSetup, :(using FillArrays))
 doctest(FillArrays; manual = false)
 
-using Aqua
-@testset "Project quality" begin
-    Aqua.test_all(FillArrays;
-        # https://github.com/JuliaArrays/FillArrays.jl/issues/105#issuecomment-1582516319
-        ambiguities=(; broken=true),
-    )
-end
+include("aqua.jl")
 
 include("infinitearrays.jl")
 import .InfiniteArrays

From 1788108049c6c105fd4a2ea76490930be812355a Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Mon, 22 Jan 2024 08:56:27 +0530
Subject: [PATCH 07/35] Bump actions/cache from 3 to 4 (#342)

Bumps [actions/cache](https://github.com/actions/cache) from 3 to 4.
- [Release notes](https://github.com/actions/cache/releases)
- [Changelog](https://github.com/actions/cache/blob/main/RELEASES.md)
- [Commits](https://github.com/actions/cache/compare/v3...v4)

---
updated-dependencies:
- dependency-name: actions/cache
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
---
 .github/workflows/docs.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 26905e08..c2b3b768 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -23,7 +23,7 @@ jobs:
         with:
           version: ${{ matrix.julia-version }}
       - name: Cache artifacts
-        uses: actions/cache@v3
+        uses: actions/cache@v4
         env:
           cache-name: cache-artifacts
         with:

From 4d4173f9281a9fbf8cd5faabd8c79ebcfa13b3c2 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Fri, 2 Feb 2024 02:30:13 +0530
Subject: [PATCH 08/35] Pretty-printing for OneElement

---
 src/oneelement.jl | 7 +++++++
 test/runtests.jl  | 9 +++++++++
 2 files changed, 16 insertions(+)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index 46b2b471..4937ee34 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -174,3 +174,10 @@ end
 function broadcasted(::DefaultArrayStyle{N}, ::typeof(\), x::Number, r::OneElement{<:Any,N}) where {N}
     OneElement(x \ r.val, r.ind, axes(r))
 end
+
+# show
+_maybesize(t::Tuple{Base.OneTo{Int}, Vararg{Base.OneTo{Int}}}) = size.(t,1)
+_maybesize(t) = t
+Base.show(io::IO, A::OneElement) = print(io, OneElement, "(", A.val, ", ", A.ind, ", ", _maybesize(axes(A)), ")")
+Base.show(io::IO, A::OneElement{<:Any,1,Tuple{Int},Tuple{Base.OneTo{Int}}}) =
+    print(io, OneElement, "(", A.val, ", ", A.ind[1], ", ", size(A,1), ")")
diff --git a/test/runtests.jl b/test/runtests.jl
index be4b8fec..b6fb3341 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2226,6 +2226,15 @@ end
             @test n .\ v == n .\ w
         end
     end
+
+    @testset "show" begin
+        B = OneElement(2, (1, 2), (3, 4))
+        @test repr(B) == "OneElement(2, (1, 2), (3, 4))"
+        B = OneElement(2, 1, 3)
+        @test repr(B) == "OneElement(2, 1, 3)"
+        B = OneElement(2, (1, 2), (Base.IdentityUnitRange(1:1), Base.IdentityUnitRange(2:2)))
+        @test repr(B) == "OneElement(2, (1, 2), (Base.IdentityUnitRange(1:1), Base.IdentityUnitRange(2:2)))"
+    end
 end
 
 @testset "repeat" begin

From 99278fa79f0e108d290e45718fea3c539706e2e3 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Fri, 2 Feb 2024 23:02:55 +0530
Subject: [PATCH 09/35] Add missing matrix multiplication methods involving
 OneElement (#347)

* Add missing matrix multiplication methods involving OneElement

* multiplications with Diagonal

* Add suggested comment to `__muloneel!`

Co-authored-by: Frames White <oxinabox@ucc.asn.au>

---------

Co-authored-by: Frames White <oxinabox@ucc.asn.au>
---
 src/oneelement.jl | 177 +++++++++++++++++++++++++++++++++++++++++++---
 test/runtests.jl  |  87 +++++++++++++++++++++--
 2 files changed, 248 insertions(+), 16 deletions(-)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index 4937ee34..f1e2fd8e 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -93,7 +93,58 @@ function *(A::OneElementMatrix, B::AbstractFillVector)
     OneElement(val, A.ind[1], size(A,1))
 end
 
-@inline function __mulonel!(y, A, x, alpha, beta)
+# Special matrix types
+
+function *(A::OneElementMatrix, D::Diagonal)
+    check_matmul_sizes(A, D)
+    nzcol = A.ind[2]
+    val = if nzcol in axes(D,1)
+        A.val * D[nzcol, nzcol]
+    else
+        A.val * zero(eltype(D))
+    end
+    OneElement(val, A.ind, size(A))
+end
+function *(D::Diagonal, A::OneElementMatrix)
+    check_matmul_sizes(D, A)
+    nzrow = A.ind[1]
+    val = if nzrow in axes(D,2)
+        D[nzrow, nzrow] * A.val
+    else
+        zero(eltype(D)) * A.val
+    end
+    OneElement(val, A.ind, size(A))
+end
+
+# Inplace multiplication
+
+# We use this for out overloads for _mul! for OneElement because its more efficient
+# due to how efficient 2 arg mul is when one or more of the args are OneElement
+function __mulonel!(C, A, B, alpha, beta)
+    ABα = A * B * alpha
+    if iszero(beta)
+        C .= ABα
+    else
+        C .= ABα .+ C .* beta
+    end
+    return C
+end
+# These methods remove the ambituity in _mul!. This isn't strictly necessary, but this makes Aqua happy.
+function _mul!(C::AbstractVector, A::OneElementMatrix, B::OneElementVector, alpha, beta)
+    __mulonel!(C, A, B, alpha, beta)
+end
+function _mul!(C::AbstractMatrix, A::OneElementMatrix, B::OneElementMatrix, alpha, beta)
+    __mulonel!(C, A, B, alpha, beta)
+end
+
+function mul!(C::AbstractMatrix, A::OneElementMatrix, B::OneElementMatrix, alpha::Number, beta::Number)
+    _mul!(C, A, B, alpha, beta)
+end
+function mul!(C::AbstractVector, A::OneElementMatrix, B::OneElementVector, alpha::Number, beta::Number)
+    _mul!(C, A, B, alpha, beta)
+end
+
+@inline function __mul!(y, A::AbstractMatrix, x::OneElement, alpha, beta)
     αx = alpha * x.val
     ind1 = x.ind[1]
     if iszero(beta)
@@ -104,19 +155,19 @@ end
     return y
 end
 
-function _mulonel!(y, A, x::OneElementVector, alpha::Number, beta::Number)
+function _mul!(y::AbstractVector, A::AbstractMatrix, x::OneElementVector, alpha, beta)
     check_matmul_sizes(y, A, x)
-    if x.ind[1] ∉ axes(x,1) # in this case x is all zeros
+    if iszero(getindex_value(x))
         mul!(y, A, Zeros{eltype(x)}(axes(x)), alpha, beta)
         return y
     end
-    __mulonel!(y, A, x, alpha, beta)
+    __mul!(y, A, x, alpha, beta)
     y
 end
 
-function _mulonel!(C, A, B::OneElementMatrix, alpha::Number, beta::Number)
+function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::OneElementMatrix, alpha, beta)
     check_matmul_sizes(C, A, B)
-    if B.ind[1] ∉ axes(B,1) || B.ind[2] ∉ axes(B,2) # in this case x is all zeros
+    if iszero(getindex_value(B))
         mul!(C, A, Zeros{eltype(B)}(axes(B)), alpha, beta)
         return C
     end
@@ -127,24 +178,128 @@ function _mulonel!(C, A, B::OneElementMatrix, alpha::Number, beta::Number)
         view(C, :, B.ind[2]+1:size(C,2)) .*= beta
     end
     y = view(C, :, B.ind[2])
-    __mulonel!(y, A, B, alpha, beta)
+    __mul!(y, A, B, alpha, beta)
+    C
+end
+function _mul!(C::AbstractMatrix, A::Diagonal, B::OneElementMatrix, alpha, beta)
+    check_matmul_sizes(C, A, B)
+    if iszero(getindex_value(B))
+        mul!(C, A, Zeros{eltype(B)}(axes(B)), alpha, beta)
+        return C
+    end
+    if iszero(beta)
+        C .= zero(eltype(C))
+    else
+        view(C, :, 1:B.ind[2]-1) .*= beta
+        view(C, :, B.ind[2]+1:size(C,2)) .*= beta
+    end
+    ABα = A * B * alpha
+    nzrow, nzcol = B.ind
+    if iszero(beta)
+        C[B.ind...] = ABα[B.ind...]
+    else
+        y = view(C, :, nzcol)
+        y .= view(ABα, :, nzcol) .+ y .* beta
+    end
+    C
+end
+
+function _mul!(C::AbstractMatrix, A::OneElementMatrix, B::AbstractMatrix, alpha, beta)
+    check_matmul_sizes(C, A, B)
+    if iszero(getindex_value(A))
+        mul!(C, Zeros{eltype(A)}(axes(A)), B, alpha, beta)
+        return C
+    end
+    if iszero(beta)
+        C .= zero(eltype(C))
+    else
+        view(C, 1:A.ind[1]-1, :) .*= beta
+        view(C, A.ind[1]+1:size(C,1), :) .*= beta
+    end
+    y = view(C, A.ind[1], :)
+    ind2 = A.ind[2]
+    Aval = A.val
+    if iszero(beta)
+        y .= Aval .* view(B, ind2, :) .* alpha
+    else
+        y .= Aval .* view(B, ind2, :) .* alpha .+ y .* beta
+    end
+    C
+end
+function _mul!(C::AbstractMatrix, A::OneElementMatrix, B::Diagonal, alpha, beta)
+    check_matmul_sizes(C, A, B)
+    if iszero(getindex_value(A))
+        mul!(C, Zeros{eltype(A)}(axes(A)), B, alpha, beta)
+        return C
+    end
+    if iszero(beta)
+        C .= zero(eltype(C))
+    else
+        view(C, 1:A.ind[1]-1, :) .*= beta
+        view(C, A.ind[1]+1:size(C,1), :) .*= beta
+    end
+    ABα = A * B * alpha
+    nzrow, nzcol = A.ind
+    if iszero(beta)
+        C[A.ind...] = ABα[A.ind...]
+    else
+        y = view(C, nzrow, :)
+        y .= view(ABα, nzrow, :) .+ y .* beta
+    end
+    C
+end
+
+function _mul!(C::AbstractVector, A::OneElementMatrix, B::AbstractVector, alpha, beta)
+    check_matmul_sizes(C, A, B)
+    if iszero(getindex_value(A))
+        mul!(C, Zeros{eltype(A)}(axes(A)), B, alpha, beta)
+        return C
+    end
+    nzrow, nzcol = A.ind
+    if iszero(beta)
+        C .= zero(eltype(C))
+    else
+        view(C, 1:nzrow-1) .*= beta
+        view(C, nzrow+1:size(C,1)) .*= beta
+    end
+    Aval = A.val
+    if iszero(beta)
+        C[nzrow] = Aval * B[nzcol] * alpha
+    else
+        C[nzrow] = Aval * B[nzcol] * alpha + C[nzrow] * beta
+    end
     C
 end
 
 for MT in (:StridedMatrix, :(Transpose{<:Any, <:StridedMatrix}), :(Adjoint{<:Any, <:StridedMatrix}))
     @eval function mul!(y::StridedVector, A::$MT, x::OneElementVector, alpha::Number, beta::Number)
-        _mulonel!(y, A, x, alpha, beta)
+        _mul!(y, A, x, alpha, beta)
     end
+end
+for MT in (:StridedMatrix, :(Transpose{<:Any, <:StridedMatrix}), :(Adjoint{<:Any, <:StridedMatrix}),
+            :Diagonal)
     @eval function mul!(C::StridedMatrix, A::$MT, B::OneElementMatrix, alpha::Number, beta::Number)
-        _mulonel!(C, A, B, alpha, beta)
+        _mul!(C, A, B, alpha, beta)
     end
+    @eval function mul!(C::StridedMatrix, A::OneElementMatrix, B::$MT, alpha::Number, beta::Number)
+        _mul!(C, A, B, alpha, beta)
+    end
+end
+function mul!(C::StridedVector, A::OneElementMatrix, B::StridedVector, alpha::Number, beta::Number)
+    _mul!(C, A, B, alpha, beta)
 end
 
 function mul!(y::AbstractVector, A::AbstractFillMatrix, x::OneElementVector, alpha::Number, beta::Number)
-    _mulonel!(y, A, x, alpha, beta)
+    _mul!(y, A, x, alpha, beta)
 end
 function mul!(C::AbstractMatrix, A::AbstractFillMatrix, B::OneElementMatrix, alpha::Number, beta::Number)
-    _mulonel!(C, A, B, alpha, beta)
+    _mul!(C, A, B, alpha, beta)
+end
+function mul!(C::AbstractVector, A::OneElementMatrix, B::AbstractFillVector, alpha::Number, beta::Number)
+    _mul!(C, A, B, alpha, beta)
+end
+function mul!(C::AbstractMatrix, A::OneElementMatrix, B::AbstractFillMatrix, alpha::Number, beta::Number)
+    _mul!(C, A, B, alpha, beta)
 end
 
 # adjoint/transpose
diff --git a/test/runtests.jl b/test/runtests.jl
index b6fb3341..9d6d5378 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2113,14 +2113,16 @@ end
 
     @testset "matmul" begin
         A = reshape(Float64[1:9;], 3, 3)
+        v = reshape(Float64[1:3;], 3)
         testinds(w::AbstractArray) = testinds(size(w))
         testinds(szw::Tuple{Int}) = (szw .- 1, szw .+ 1)
         function testinds(szA::Tuple{Int,Int})
             (szA .- 1, szA .+ (-1,0), szA .+ (0,-1), szA .+ 1, szA .+ (1,-1), szA .+ (-1,1))
         end
-        function test_A_mul_OneElement(A, (w, w2))
-            @testset for ind in testinds(w)
-                x = OneElement(3, ind, size(w))
+        # test matvec if w is a vector, or matmat if w is a matrix
+        function test_mat_mul_OneElement(A, (w, w2), sz)
+            @testset for ind in testinds(sz)
+                x = OneElement(3, ind, sz)
                 xarr = Array(x)
                 Axarr = A * xarr
                 Aadjxarr = A' * xarr
@@ -2143,15 +2145,69 @@ end
                 @test mul!(w2, F, x, 1.0, 1.0) ≈ Array(F) * xarr .+ 1
             end
         end
+        function test_OneElementMatrix_mul_mat(A, (w, w2), sz)
+            @testset for ind in testinds(sz)
+                O = OneElement(3, ind, sz)
+                Oarr = Array(O)
+                OarrA = Oarr * A
+                OarrAadj = Oarr * A'
+
+                @test O * A ≈ OarrA
+                @test O * A' ≈ OarrAadj
+                @test O * transpose(A) ≈ Oarr * transpose(A)
+
+                @test mul!(w, O, A) ≈ OarrA
+                # check columnwise to ensure zero columns
+                @test all(((c1, c2),) -> c1 ≈ c2, zip(eachcol(w), eachcol(OarrA)))
+                @test mul!(w, O, A') ≈ OarrAadj
+                w .= 1
+                @test mul!(w, O, A, 1.0, 2.0) ≈ OarrA .+ 2
+                w .= 1
+                @test mul!(w, O, A', 1.0, 2.0) ≈ OarrAadj .+ 2
+
+                F = Fill(3, size(A))
+                w2 .= 1
+                @test mul!(w2, O, F, 1.0, 1.0) ≈ Oarr * Array(F) .+ 1
+            end
+        end
+        function test_OneElementMatrix_mul_vec(v, (w, w2), sz)
+            @testset for ind in testinds(sz)
+                O = OneElement(3, ind, sz)
+                Oarr = Array(O)
+                Oarrv = Oarr * v
+
+                @test O * v == Oarrv
+
+                @test mul!(w, O, v) == Oarrv
+                # check rowwise to ensure zero rows
+                @test all(((r1, r2),) -> r1 == r2, zip(eachrow(w), eachrow(Oarrv)))
+                w .= 1
+                @test mul!(w, O, v, 1.0, 2.0) == Oarrv .+ 2
+
+                F = Fill(3, size(v))
+                w2 .= 1
+                @test mul!(w2, O, F, 1.0, 1.0) == Oarr * Array(F) .+ 1
+            end
+        end
         @testset "Matrix * OneElementVector" begin
             w = zeros(size(A,1))
             w2 = MVector{length(w)}(w)
-            test_A_mul_OneElement(A, (w, w2))
+            test_mat_mul_OneElement(A, (w, w2), size(w))
         end
         @testset "Matrix * OneElementMatrix" begin
             C = zeros(size(A))
             C2 = MMatrix{size(C)...}(C)
-            test_A_mul_OneElement(A, (C, C2))
+            test_mat_mul_OneElement(A, (C, C2), size(C))
+        end
+        @testset "OneElementMatrix * Vector" begin
+            w = zeros(size(v))
+            w2 = MVector{size(v)...}(v)
+            test_OneElementMatrix_mul_vec(v, (w, w2), size(A))
+        end
+        @testset "OneElementMatrix * Matrix" begin
+            C = zeros(size(A))
+            C2 = MMatrix{size(C)...}(C)
+            test_OneElementMatrix_mul_mat(A, (C, C2), size(A))
         end
         @testset "OneElementMatrix * OneElement" begin
             @testset for ind in testinds(A)
@@ -2159,10 +2215,14 @@ end
                 v = OneElement(4, ind[2], size(A,1))
                 @test O * v isa OneElement
                 @test O * v == Array(O) * Array(v)
+                @test mul!(ones(size(O,1)), O, v) == O * v
+                @test mul!(ones(size(O,1)), O, v, 2, 1) == 2 * O * v .+ 1
 
                 B = OneElement(4, ind, size(A))
                 @test O * B isa OneElement
                 @test O * B == Array(O) * Array(B)
+                @test mul!(ones(size(O,1), size(B,2)), O, B) == O * B
+                @test mul!(ones(size(O,1), size(B,2)), O, B, 2, 1) == 2 * O * B .+ 1
             end
 
             @test OneElement(3, (2,3), (5,4)) * OneElement(2, 2, 4) == Zeros(5)
@@ -2191,6 +2251,23 @@ end
             B = Zeros(4)
             @test A * B === Zeros(5)
         end
+        @testset "Diagonal and OneElementMatrix" begin
+            for ind in ((2,3), (2,2), (10,10))
+                O = OneElement(3, ind, (4,3))
+                Oarr = Array(O)
+                C = zeros(size(O))
+                D = Diagonal(axes(O,1))
+                @test D * O == D * Oarr
+                @test mul!(C, D, O) == D * O
+                C .= 1
+                @test mul!(C, D, O, 2, 2) == 2 * D * O .+ 2
+                D = Diagonal(axes(O,2))
+                @test O * D == Oarr * D
+                @test mul!(C, O, D) == O * D
+                C .= 1
+                @test mul!(C, O, D, 2, 2) == 2 * O * D .+ 2
+            end
+        end
     end
 
     @testset "multiplication/division by a number" begin

From b04c00d1644d27c8db57837ed3de04d5431bb78d Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Mon, 5 Feb 2024 08:26:10 +0530
Subject: [PATCH 10/35] Bump codecov/codecov-action from 3 to 4 (#349)

Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 3 to 4.
- [Release notes](https://github.com/codecov/codecov-action/releases)
- [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md)
- [Commits](https://github.com/codecov/codecov-action/compare/v3...v4)

---
updated-dependencies:
- dependency-name: codecov/codecov-action
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
---
 .github/workflows/ci.yml         | 2 +-
 .github/workflows/downstream.yml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 40668356..b7917334 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -46,6 +46,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index b804e802..7338cf17 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -74,6 +74,6 @@ jobs:
             exit(0)  # Exit immediately, as a success
           end
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
           files: lcov.info

From 05804e400cbdb6019069ef4c3a4381bf34f72ae8 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 5 Feb 2024 17:53:17 +0530
Subject: [PATCH 11/35] Short-circuit `all` for `AbstractFill` (#329)

* Short-circuit all for AbstractFill

* expand the cases in all

* Handle the fallback case through all

* Simplify implementation
---
 src/FillArrays.jl | 14 ++++++++++++--
 test/runtests.jl  |  8 ++++++++
 2 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/src/FillArrays.jl b/src/FillArrays.jl
index 3418d640..ad479a32 100644
--- a/src/FillArrays.jl
+++ b/src/FillArrays.jl
@@ -643,8 +643,18 @@ end
 
 # In particular, these make iszero(Eye(n))  efficient.
 # use any/all on scalar to get Boolean error message
-any(f::Function, x::AbstractFill) = !isempty(x) && any(f(getindex_value(x)))
-all(f::Function, x::AbstractFill) = isempty(x) || all(f(getindex_value(x)))
+function any(f::Function, x::AbstractFill)
+    isempty(x) && return false
+    # If the condition is true for one value, then it's true for all
+    fval = f(getindex_value(x))
+    any((fval,))
+end
+function all(f::Function, x::AbstractFill)
+    isempty(x) && return true
+    # If the condition is true for one value, then it's true for all
+    fval = f(getindex_value(x))
+    return all((fval,))
+end
 any(x::AbstractFill) = any(identity, x)
 all(x::AbstractFill) = all(identity, x)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 9d6d5378..3b718309 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1330,6 +1330,9 @@ end
     @test iszero(Fill(SMatrix{2,2}(0,0,0,0), 2))
     @test iszero(Fill(SMatrix{2,2}(0,0,0,1), 0))
 
+    # compile-time evaluation
+    @test @inferred((Z -> Val(iszero(Z)))(Zeros(3,3))) == Val(true)
+
     @testset "all/any" begin
         @test any(Ones{Bool}(10)) === all(Ones{Bool}(10)) === any(Fill(true,10)) === all(Fill(true,10)) === true
         @test any(Zeros{Bool}(10)) === all(Zeros{Bool}(10)) === any(Fill(false,10)) === all(Fill(false,10)) === false
@@ -1339,6 +1342,11 @@ end
         @test all(Fill(2,0))
         @test !any(Fill(2,0))
         @test any(Trues(2,0)) == any(trues(2,0))
+        @test_throws TypeError all(Fill(2,2))
+        @test all(iszero, Fill(missing,0)) === all(iszero, fill(missing,0)) === true
+        @test all(iszero, Fill(missing,2)) === all(iszero, fill(missing,2)) === missing
+        @test any(iszero, Fill(missing,0)) === any(iszero, fill(missing,0)) === false
+        @test any(iszero, Fill(missing,2)) === any(iszero, fill(missing,2)) === missing
     end
 
     @testset "Error" begin

From 9482273404fd63b6e04f31d6dbad5a7ecec6b0d4 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Wed, 7 Feb 2024 23:37:38 +0530
Subject: [PATCH 12/35] Use codecov token for the codecov action in CI (#350)

---
 .github/workflows/ci.yml         | 1 +
 .github/workflows/downstream.yml | 1 +
 2 files changed, 2 insertions(+)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index b7917334..0e6d351b 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -48,4 +48,5 @@ jobs:
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info
diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index 7338cf17..3a25d8f0 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -76,4 +76,5 @@ jobs:
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           files: lcov.info

From 2784902d54e0ad32aca871cf3ad9b69bab42b480 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 27 Feb 2024 16:30:03 +0530
Subject: [PATCH 13/35] Use axes function instead of field access (#320)

* Use axes function instead of field access

* Don't change for FillMatrix

* Add tests for custom structs

* Typo
---
 src/fillalgebra.jl   |  2 +-
 src/fillbroadcast.jl |  8 ++++----
 test/runtests.jl     | 27 +++++++++++++++++++++++++++
 3 files changed, 32 insertions(+), 5 deletions(-)

diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl
index 0fef9a97..ea454fdf 100644
--- a/src/fillalgebra.jl
+++ b/src/fillalgebra.jl
@@ -19,7 +19,7 @@ for OP in (:transpose, :adjoint)
 end
 
 permutedims(a::AbstractFillVector) = fillsimilar(a, (1, length(a)))
-permutedims(a::AbstractFillMatrix) = fillsimilar(a, reverse(a.axes))
+permutedims(a::AbstractFillMatrix) = fillsimilar(a, reverse(axes(a)))
 
 function permutedims(B::AbstractFill, perm)
     dimsB = size(B)
diff --git a/src/fillbroadcast.jl b/src/fillbroadcast.jl
index b2cb0755..d286418a 100644
--- a/src/fillbroadcast.jl
+++ b/src/fillbroadcast.jl
@@ -85,10 +85,10 @@ broadcasted(::DefaultArrayStyle, ::typeof(+), r::AbstractOnes) = r
 
 broadcasted(::DefaultArrayStyle{N}, ::typeof(conj), r::AbstractZeros{T,N}) where {T,N} = r
 broadcasted(::DefaultArrayStyle{N}, ::typeof(conj), r::AbstractOnes{T,N}) where {T,N} = r
-broadcasted(::DefaultArrayStyle{N}, ::typeof(real), r::AbstractZeros{T,N}) where {T,N} = Zeros{real(T)}(r.axes)
-broadcasted(::DefaultArrayStyle{N}, ::typeof(real), r::AbstractOnes{T,N}) where {T,N} = Ones{real(T)}(r.axes)
-broadcasted(::DefaultArrayStyle{N}, ::typeof(imag), r::AbstractZeros{T,N}) where {T,N} = Zeros{real(T)}(r.axes)
-broadcasted(::DefaultArrayStyle{N}, ::typeof(imag), r::AbstractOnes{T,N}) where {T,N} = Zeros{real(T)}(r.axes)
+broadcasted(::DefaultArrayStyle{N}, ::typeof(real), r::AbstractZeros{T,N}) where {T,N} = Zeros{real(T)}(axes(r))
+broadcasted(::DefaultArrayStyle{N}, ::typeof(real), r::AbstractOnes{T,N}) where {T,N} = Ones{real(T)}(axes(r))
+broadcasted(::DefaultArrayStyle{N}, ::typeof(imag), r::AbstractZeros{T,N}) where {T,N} = Zeros{real(T)}(axes(r))
+broadcasted(::DefaultArrayStyle{N}, ::typeof(imag), r::AbstractOnes{T,N}) where {T,N} = Zeros{real(T)}(axes(r))
 
 ### Binary broadcasting
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 3b718309..8531b461 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2012,6 +2012,33 @@ end
         @test FillArrays.getindex_value(transpose(a)) == FillArrays.unique_value(transpose(a)) == 2.0 + im
         @test convert(Fill, transpose(a)) ≡ Fill(2.0+im,1,5)
     end
+
+    @testset "custom AbstractFill types" begin
+        # implicit axes
+        struct StaticZerosVec{L,T} <: FillArrays.AbstractZeros{T,1,Tuple{SOneTo{L}}} end
+        Base.size(::StaticZerosVec{L}) where {L} = (L,)
+        Base.axes(::StaticZerosVec{L}) where {L} = (SOneTo(L),)
+        S = StaticZerosVec{3,Int}()
+        @test real.(S) == S
+        @test imag.(S) == S
+
+        struct StaticOnesVec{L,T} <: FillArrays.AbstractOnes{T,1,Tuple{SOneTo{L}}} end
+        Base.size(::StaticOnesVec{L}) where {L} = (L,)
+        Base.axes(::StaticOnesVec{L}) where {L} = (SOneTo(L),)
+        S = StaticOnesVec{3,Int}()
+        @test real.(S) == S
+        @test imag.(S) == zero(S)
+
+        struct StaticFill{S1,S2,T} <: FillArrays.AbstractFill{T,2,Tuple{SOneTo{S1},SOneTo{S2}}}
+            x :: T
+        end
+        StaticFill{S1,S2}(x::T) where {S1,S2,T} = StaticFill{S1,S2,T}(x)
+        Base.size(::StaticFill{S1,S2}) where {S1,S2} = (S1,S2)
+        Base.axes(::StaticFill{S1,S2}) where {S1,S2} = (SOneTo(S1), SOneTo(S2))
+        FillArrays.getindex_value(S::StaticFill) = S.x
+        S = StaticFill{2,3}(2)
+        @test permutedims(S) == Fill(2, reverse(size(S)))
+    end
 end
 
 @testset "Statistics" begin

From abb2481806c5f54951bd925980398fecdf5f45e6 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Wed, 13 Mar 2024 17:22:45 +0530
Subject: [PATCH 14/35] Point to `getindex_value` in `AbstractFill` docstring
 (#323)

* Point to getindex_value in AbstractFill docstring

* Update phrasing
---
 src/FillArrays.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/FillArrays.jl b/src/FillArrays.jl
index ad479a32..581e8252 100644
--- a/src/FillArrays.jl
+++ b/src/FillArrays.jl
@@ -23,7 +23,8 @@ import Base: oneto
 """
     AbstractFill{T, N, Axes} <: AbstractArray{T, N}
 
-Supertype for lazy array types whose entries are all equal to constant.
+Supertype for lazy array types whose entries are all equal.
+Subtypes of `AbstractFill` should implement [`FillArrays.getindex_value`](@ref) to return the value of the entries.
 """
 abstract type AbstractFill{T, N, Axes} <: AbstractArray{T, N} end
 const AbstractFillVector{T} = AbstractFill{T,1}
@@ -141,7 +142,7 @@ Fill{T,0}(x, ::Tuple{}) where T = Fill{T,0,Tuple{}}(convert(T, x)::T, ()) # ambi
 @inline size(F::Fill) = map(length, F.axes)
 
 """
-    getindex_value(F::AbstractFill)
+    FillArrays.getindex_value(F::AbstractFill)
 
 Return the value that `F` is filled with.
 

From 3369ca06d4e1bd084edb02f05ce847210716ee66 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 2 Apr 2024 11:30:01 +0530
Subject: [PATCH 15/35] Make pkgeval badge point to Nanosoldier report (#354)

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index dd4c6cfc..0ec1788a 100644
--- a/README.md
+++ b/README.md
@@ -6,7 +6,7 @@
 [![codecov](https://codecov.io/gh/JuliaArrays/FillArrays.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaArrays/FillArrays.jl)
 [![deps](https://juliahub.com/docs/FillArrays/deps.svg)](https://juliahub.com/ui/Packages/FillArrays/2Dg1l?t=2)
 [![version](https://juliahub.com/docs/FillArrays/version.svg)](https://juliahub.com/ui/Packages/FillArrays/2Dg1l)
-[![pkgeval](https://juliahub.com/docs/FillArrays/pkgeval.svg)](https://juliahub.com/ui/Packages/FillArrays/2Dg1l)
+[![pkgeval](https://juliahub.com/docs/FillArrays/pkgeval.svg)](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/report.html)
 [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
 Julia package to lazily represent matrices filled with a single entry,

From f174c7b0bcac5023751a5ac07770d3d52f906d9b Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 2 Apr 2024 12:09:42 +0530
Subject: [PATCH 16/35] Ignore non-code files in CI (#355)

---
 .github/workflows/Invalidations.yml |  4 ++++
 .github/workflows/ci.yml            | 10 ++++++++++
 .github/workflows/docs.yml          |  8 ++++++++
 .github/workflows/downstream.yml    |  4 ++++
 4 files changed, 26 insertions(+)

diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml
index 28b9ce2f..e4ce0dc8 100644
--- a/.github/workflows/Invalidations.yml
+++ b/.github/workflows/Invalidations.yml
@@ -2,6 +2,10 @@ name: Invalidations
 
 on:
   pull_request:
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
 
 concurrency:
   # Skip intermediate builds: always.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0e6d351b..a75d6423 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -3,7 +3,17 @@ on:
   push:
     branches:
       - master
+    tags: [v*]
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
   pull_request:
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
+
 
 concurrency:
   group: build-${{ github.event.pull_request.number || github.ref }}-${{ github.workflow }}
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index c2b3b768..1a916082 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -2,11 +2,19 @@ name: Documentation
 
 on:
   pull_request:
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
   push:
     branches:
       - 'master'
       - 'release-'
     tags: '*'
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
   release:
     types: [published]
 
diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index 3a25d8f0..454f4777 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -3,6 +3,10 @@ on:
   push:
     branches: [master]
     tags: [v*]
+    paths-ignore:
+      - 'LICENSE'
+      - 'README.md'
+      - '.github/workflows/TagBot.yml'
   pull_request:
     paths-ignore:
       - 'LICENSE'

From 002f56c7562badc7c6bce8d58637f02799b537fa Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Wed, 3 Apr 2024 16:46:27 +0530
Subject: [PATCH 17/35] Test floating-point arrays using isapprox (#318)

---
 test/runtests.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 8531b461..380c15f8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -378,22 +378,22 @@ end
 # type, and produce numerically correct results.
 as_array(x::AbstractArray) = Array(x)
 as_array(x::UniformScaling) = x
-equal_or_undef(a::Number, b::Number) = (a == b) || isequal(a, b)
-equal_or_undef(a, b) = all(equal_or_undef.(a, b))
+isapprox_or_undef(a::Number, b::Number) = (a ≈ b) || isequal(a, b)
+isapprox_or_undef(a, b) = all(((x,y),) -> isapprox_or_undef(x,y), zip(a, b))
 function test_addition_subtraction_dot(As, Bs, Tout::Type)
     for A in As, B in Bs
         @testset "$(typeof(A)) and $(typeof(B))" begin
             @test @inferred(A + B) isa Tout{promote_type(eltype(A), eltype(B))}
-            @test equal_or_undef(as_array(A + B), as_array(A) + as_array(B))
+            @test isapprox_or_undef(as_array(A + B), as_array(A) + as_array(B))
 
             @test @inferred(A - B) isa Tout{promote_type(eltype(A), eltype(B))}
-            @test equal_or_undef(as_array(A - B), as_array(A) - as_array(B))
+            @test isapprox_or_undef(as_array(A - B), as_array(A) - as_array(B))
 
             @test @inferred(B + A) isa Tout{promote_type(eltype(B), eltype(A))}
-            @test equal_or_undef(as_array(B + A), as_array(B) + as_array(A))
+            @test isapprox_or_undef(as_array(B + A), as_array(B) + as_array(A))
 
             @test @inferred(B - A) isa Tout{promote_type(eltype(B), eltype(A))}
-            @test equal_or_undef(as_array(B - A), as_array(B) - as_array(A))
+            @test isapprox_or_undef(as_array(B - A), as_array(B) - as_array(A))
 
             # Julia 1.6 doesn't support dot(UniformScaling)
             if VERSION < v"1.6.0" || VERSION >= v"1.8.0"

From cc810ac5db4c8e65b29fe3876ff0a5e6a1d91474 Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Mon, 8 Apr 2024 09:42:31 +0530
Subject: [PATCH 18/35] Bump julia-actions/setup-julia from 1 to 2 (#357)

Bumps [julia-actions/setup-julia](https://github.com/julia-actions/setup-julia) from 1 to 2.
- [Release notes](https://github.com/julia-actions/setup-julia/releases)
- [Commits](https://github.com/julia-actions/setup-julia/compare/v1...v2)

---
updated-dependencies:
- dependency-name: julia-actions/setup-julia
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
---
 .github/workflows/CompatHelper.yml  | 2 +-
 .github/workflows/Invalidations.yml | 2 +-
 .github/workflows/ci.yml            | 2 +-
 .github/workflows/downstream.yml    | 2 +-
 4 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 8dd16866..eb40a561 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -15,7 +15,7 @@ jobs:
         run: which julia
         continue-on-error: true
       - name: Install Julia, but only if it is not already available in the PATH
-        uses: julia-actions/setup-julia@v1
+        uses: julia-actions/setup-julia@v2
         with:
           version: '1'
           arch: ${{ runner.arch }}
diff --git a/.github/workflows/Invalidations.yml b/.github/workflows/Invalidations.yml
index e4ce0dc8..6f33a462 100644
--- a/.github/workflows/Invalidations.yml
+++ b/.github/workflows/Invalidations.yml
@@ -20,7 +20,7 @@ jobs:
     if: github.base_ref == github.event.repository.default_branch
     runs-on: ubuntu-latest
     steps:
-    - uses: julia-actions/setup-julia@v1
+    - uses: julia-actions/setup-julia@v2
       with:
         version: '1'
     - uses: actions/checkout@v4
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index a75d6423..bf59dbfb 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -48,7 +48,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
diff --git a/.github/workflows/downstream.yml b/.github/workflows/downstream.yml
index 454f4777..30fde31f 100644
--- a/.github/workflows/downstream.yml
+++ b/.github/workflows/downstream.yml
@@ -50,7 +50,7 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.julia-version }}
           arch: x64

From 4c5018636650629046ca4a01b32b13d3807ec448 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 15 Apr 2024 18:46:25 +0530
Subject: [PATCH 19/35] Remove single-argument map tests on nightly (#358)

* Remove single-argument map tests on nightly

* Add comment on method removal

Co-authored-by: Frames White <oxinabox@ucc.asn.au>

---------

Co-authored-by: Frames White <oxinabox@ucc.asn.au>
---
 test/runtests.jl | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 380c15f8..02992733 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1132,8 +1132,13 @@ end
     @test_throws DimensionMismatch map(+, x2', x2)
 
     # Issue https://github.com/JuliaArrays/FillArrays.jl/issues/179
-    @test map(() -> "ok") == "ok"  # was MethodError: reducing over an empty collection is not allowed
-    @test mapreduce(() -> "ok", *) == "ok"
+    if VERSION < v"1.11.0-"  # In 1.11, 1-arg map & mapreduce was removed
+        @test map(() -> "ok") == "ok"  # was MethodError: reducing over an empty collection is not allowed
+        @test mapreduce(() -> "ok", *) == "ok"
+    else
+        @test_throws "no method matching map" map(() -> "ok")
+        @test_throws "no method matching map" mapreduce(() -> "ok", *)
+    end
 end
 
 @testset "mapreduce" begin

From 7177c59ffbfa16e51de16bd0d5bf5e3b82839016 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 15 Apr 2024 23:33:44 +0530
Subject: [PATCH 20/35] ignore versioned manifest file (#359)

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index 28d0ac95..bfd7878c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,4 +4,5 @@
 deps/deps.jl
 .DS_Store
 Manifest.toml
+Manifest-v*.*.toml
 docs/build

From 8509b6b1ed24570daefec317337669c02466603a Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 22 Apr 2024 18:47:22 +0530
Subject: [PATCH 21/35] Bugfix in matrix multiplication involving adj/trans
 (#360)

* Bugfix in matrix multiplication

* _mulfill implementation for all matrix orderings

* Loop in copyfirstrow!

* Revert some unnecessary changes
---
 src/fillalgebra.jl | 92 ++++++++++++++++++++++++++++++----------------
 test/runtests.jl   | 10 +++++
 2 files changed, 70 insertions(+), 32 deletions(-)

diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl
index ea454fdf..a887eca2 100644
--- a/src/fillalgebra.jl
+++ b/src/fillalgebra.jl
@@ -212,7 +212,8 @@ for (T, f) in ((:Adjoint, :adjoint), (:Transpose, :transpose))
     end
 end
 
-function mul!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractFillMatrix, alpha::Number, beta::Number)
+# unnecessary indirection, added for ambiguity resolution
+function _mulfill!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractFillMatrix, alpha, beta)
     check_matmul_sizes(C, A, B)
     ABα = getindex_value(A) * getindex_value(B) * alpha * size(B,1)
     if iszero(beta)
@@ -220,7 +221,12 @@ function mul!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractFillMatrix, a
     else
         C .= ABα .+ C .* beta
     end
-    C
+    return C
+end
+
+function mul!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractFillMatrix, alpha::Number, beta::Number)
+    _mulfill!(C, A, B, alpha, beta)
+    return C
 end
 
 function copyfirstcol!(C)
@@ -229,50 +235,72 @@ function copyfirstcol!(C)
     end
     return C
 end
-function copyfirstcol!(C::Union{Adjoint, Transpose})
-    # in this case, we copy the first row of the parent to others
-    Cp = parent(C)
-    for colind in axes(Cp, 2)
-        Cp[2:end, colind] .= Cp[1, colind]
+
+_firstcol(C::AbstractMatrix) = first(eachcol(C))
+
+function copyfirstrow!(C)
+    # C[begin+1:end, ind] .= permutedims(_firstrow(C))
+    # we loop here as the aliasing check isn't smart enough to
+    # detect that the two sides don't alias, and ends up materializing the RHS
+    for (ind, v) in pairs(_firstrow(C))
+        C[begin+1:end, ind] .= Ref(v)
     end
     return C
 end
+_firstrow(C::AbstractMatrix) = first(eachrow(C))
 
-_firstcol(C::AbstractMatrix) = view(C, :, 1)
-_firstcol(C::Union{Adjoint, Transpose}) = view(parent(C), 1, :)
-
-function _mulfill!(C, A, B::AbstractFillMatrix, alpha, beta)
+function _mulfill!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractFillMatrix, alpha, beta)
+    check_matmul_sizes(C, A, B)
+    iszero(size(B,2)) && return C # no columns in B and C, empty matrix
+    if iszero(beta)
+        # the mat-vec product sums along the rows of A
+        mul!(_firstcol(C), A, _firstcol(B), alpha, beta)
+        copyfirstcol!(C)
+    else
+        # the mat-vec product sums along the rows of A, which produces the first column of ABα
+        # allocate a temporary column vector to store the result
+        v = A * (_firstcol(B) * alpha)
+        C .= v .+ C .* beta
+    end
+    return C
+end
+function _mulfill!(C::AbstractMatrix, A::AbstractFillMatrix, B::AbstractMatrix, alpha, beta)
     check_matmul_sizes(C, A, B)
-    if iszero(size(B,2))
-        return rmul!(C, beta)
+    iszero(size(A,1)) && return C # no rows in A and C, empty matrix
+    Aval = getindex_value(A)
+    if iszero(beta)
+        Crow = _firstrow(C)
+        # sum along the columns of B
+        Crow .= Ref(Aval) .* sum.(eachcol(B)) .* alpha
+        copyfirstrow!(C)
+    else
+        # sum along the columns of B, and allocate the result.
+        # This is the first row of ABα
+        ABα_row = Ref(Aval) .* sum.(eachcol(B)) .* alpha
+        C .= permutedims(ABα_row) .+ C .* beta
     end
-    mul!(_firstcol(C), A, view(B, :, 1), alpha, beta)
-    copyfirstcol!(C)
-    C
+    return C
 end
 
 function mul!(C::StridedMatrix, A::StridedMatrix, B::AbstractFillMatrix, alpha::Number, beta::Number)
     _mulfill!(C, A, B, alpha, beta)
+    return C
 end
-
-for T in (:Adjoint, :Transpose)
-    @eval function mul!(C::StridedMatrix, A::$T{<:Any, <:StridedMatrix}, B::AbstractFillMatrix, alpha::Number, beta::Number)
-        _mulfill!(C, A, B, alpha, beta)
-    end
-end
-
 function mul!(C::StridedMatrix, A::AbstractFillMatrix, B::StridedMatrix, alpha::Number, beta::Number)
-    check_matmul_sizes(C, A, B)
-    for (colC, colB) in zip(eachcol(C), eachcol(B))
-        mul!(colC, A, colB, alpha, beta)
-    end
-    C
+    _mulfill!(C, A, B, alpha, beta)
+    return C
 end
 
-for (T, f) in ((:Adjoint, :adjoint), (:Transpose, :transpose))
-    @eval function mul!(C::StridedMatrix, A::AbstractFillMatrix, B::$T{<:Any, <:StridedMatrix}, alpha::Number, beta::Number)
-        _mulfill!($f(C), $f(B), $f(A), alpha, beta)
-        C
+for T in (:Adjoint, :Transpose)
+    @eval begin
+        function mul!(C::StridedMatrix, A::$T{<:Any, <:StridedMatrix}, B::AbstractFillMatrix, alpha::Number, beta::Number)
+            _mulfill!(C, A, B, alpha, beta)
+            return C
+        end
+        function mul!(C::StridedMatrix, A::AbstractFillMatrix, B::$T{<:Any, <:StridedMatrix}, alpha::Number, beta::Number)
+            _mulfill!(C, A, B, alpha, beta)
+            return C
+        end
     end
 end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 02992733..34c77108 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1716,6 +1716,16 @@ end
         @test transpose(A)*fillmat ≈ transpose(A)*Array(fillmat)
         @test adjoint(A)*fillvec ≈ adjoint(A)*Array(fillvec)
         @test adjoint(A)*fillmat ≈ adjoint(A)*Array(fillmat)
+
+        # inplace C = F * B' * alpha + C * beta
+        F = Fill(fv, m, k)
+        A = Array(F)
+        B = rand(T, n, k)
+        C = rand(T, m, n)
+        @testset for f in (adjoint, transpose)
+            @test mul!(copy(C), F, f(B)) ≈ mul!(copy(C), A, f(B))
+            @test mul!(copy(C), F, f(B), 1.0, 2.0) ≈ mul!(copy(C), A, f(B), 1.0, 2.0)
+        end
     end
 
     @testset "non-commutative" begin

From b0ee65f0d05516dcde1c21c911297372d82352b5 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 23 Apr 2024 11:27:26 +0530
Subject: [PATCH 22/35] Bump version to v1.10.1 (#361)

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d980a8a5..e6fff2d1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.10.0"
+version = "1.10.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 873437171f2320e2cf4e99ef44682cbfbae40962 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 25 Apr 2024 20:01:31 +0530
Subject: [PATCH 23/35] Fix OneElement multiplication with array elements
 (#335)

* Fix OneElement multiplication with array elements

* Fix matmul for array elements in OneElMat * StridedMat
---
 src/oneelement.jl | 67 +++++++++++++++---------------------
 test/runtests.jl  | 87 +++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 115 insertions(+), 39 deletions(-)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index f1e2fd8e..28bc2c4b 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -145,12 +145,12 @@ function mul!(C::AbstractVector, A::OneElementMatrix, B::OneElementVector, alpha
 end
 
 @inline function __mul!(y, A::AbstractMatrix, x::OneElement, alpha, beta)
-    αx = alpha * x.val
+    xα = Ref(x.val * alpha)
     ind1 = x.ind[1]
     if iszero(beta)
-        y .= αx .* view(A, :, ind1)
+        y .= view(A, :, ind1) .* xα
     else
-        y .= αx .* view(A, :, ind1) .+ beta .* y
+        y .= view(A, :, ind1) .* xα .+ y .* beta
     end
     return y
 end
@@ -171,13 +171,14 @@ function _mul!(C::AbstractMatrix, A::AbstractMatrix, B::OneElementMatrix, alpha,
         mul!(C, A, Zeros{eltype(B)}(axes(B)), alpha, beta)
         return C
     end
+    nzrow, nzcol = B.ind
     if iszero(beta)
-        C .= zero(eltype(C))
+        C .= Ref(zero(eltype(C)))
     else
-        view(C, :, 1:B.ind[2]-1) .*= beta
-        view(C, :, B.ind[2]+1:size(C,2)) .*= beta
+        view(C, :, 1:nzcol-1) .*= beta
+        view(C, :, nzcol+1:size(C,2)) .*= beta
     end
-    y = view(C, :, B.ind[2])
+    y = view(C, :, nzcol)
     __mul!(y, A, B, alpha, beta)
     C
 end
@@ -187,17 +188,14 @@ function _mul!(C::AbstractMatrix, A::Diagonal, B::OneElementMatrix, alpha, beta)
         mul!(C, A, Zeros{eltype(B)}(axes(B)), alpha, beta)
         return C
     end
-    if iszero(beta)
-        C .= zero(eltype(C))
-    else
-        view(C, :, 1:B.ind[2]-1) .*= beta
-        view(C, :, B.ind[2]+1:size(C,2)) .*= beta
-    end
-    ABα = A * B * alpha
     nzrow, nzcol = B.ind
+    ABα = A * B * alpha
     if iszero(beta)
-        C[B.ind...] = ABα[B.ind...]
+        C .= Ref(zero(eltype(C)))
+        C[nzrow, nzcol] = ABα[nzrow, nzcol]
     else
+        view(C, :, 1:nzcol-1) .*= beta
+        view(C, :, nzcol+1:size(C,2)) .*= beta
         y = view(C, :, nzcol)
         y .= view(ABα, :, nzcol) .+ y .* beta
     end
@@ -210,19 +208,16 @@ function _mul!(C::AbstractMatrix, A::OneElementMatrix, B::AbstractMatrix, alpha,
         mul!(C, Zeros{eltype(A)}(axes(A)), B, alpha, beta)
         return C
     end
-    if iszero(beta)
-        C .= zero(eltype(C))
-    else
-        view(C, 1:A.ind[1]-1, :) .*= beta
-        view(C, A.ind[1]+1:size(C,1), :) .*= beta
-    end
-    y = view(C, A.ind[1], :)
-    ind2 = A.ind[2]
+    nzrow, nzcol = A.ind
+    y = view(C, nzrow, :)
     Aval = A.val
     if iszero(beta)
-        y .= Aval .* view(B, ind2, :) .* alpha
+        C .= Ref(zero(eltype(C)))
+        y .= Ref(Aval) .* view(B, nzcol, :) .* alpha
     else
-        y .= Aval .* view(B, ind2, :) .* alpha .+ y .* beta
+        view(C, 1:nzrow-1, :) .*= beta
+        view(C, nzrow+1:size(C,1), :) .*= beta
+        y .= Ref(Aval) .* view(B, nzcol, :) .* alpha .+ y .* beta
     end
     C
 end
@@ -232,17 +227,14 @@ function _mul!(C::AbstractMatrix, A::OneElementMatrix, B::Diagonal, alpha, beta)
         mul!(C, Zeros{eltype(A)}(axes(A)), B, alpha, beta)
         return C
     end
-    if iszero(beta)
-        C .= zero(eltype(C))
-    else
-        view(C, 1:A.ind[1]-1, :) .*= beta
-        view(C, A.ind[1]+1:size(C,1), :) .*= beta
-    end
-    ABα = A * B * alpha
     nzrow, nzcol = A.ind
+    ABα = A * B * alpha
     if iszero(beta)
-        C[A.ind...] = ABα[A.ind...]
+        C .= Ref(zero(eltype(C)))
+        C[nzrow, nzcol] = ABα[nzrow, nzcol]
     else
+        view(C, 1:nzrow-1, :) .*= beta
+        view(C, nzrow+1:size(C,1), :) .*= beta
         y = view(C, nzrow, :)
         y .= view(ABα, nzrow, :) .+ y .* beta
     end
@@ -256,16 +248,13 @@ function _mul!(C::AbstractVector, A::OneElementMatrix, B::AbstractVector, alpha,
         return C
     end
     nzrow, nzcol = A.ind
-    if iszero(beta)
-        C .= zero(eltype(C))
-    else
-        view(C, 1:nzrow-1) .*= beta
-        view(C, nzrow+1:size(C,1)) .*= beta
-    end
     Aval = A.val
     if iszero(beta)
+        C .= Ref(zero(eltype(C)))
         C[nzrow] = Aval * B[nzcol] * alpha
     else
+        view(C, 1:nzrow-1) .*= beta
+        view(C, nzrow+1:size(C,1)) .*= beta
         C[nzrow] = Aval * B[nzcol] * alpha + C[nzrow] * beta
     end
     C
diff --git a/test/runtests.jl b/test/runtests.jl
index 34c77108..fe5fafa3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2318,6 +2318,93 @@ end
                 @test mul!(C, O, D, 2, 2) == 2 * O * D .+ 2
             end
         end
+        @testset "array elements" begin
+            A = [SMatrix{2,3}(1:6)*(i+j) for i in 1:3, j in 1:2]
+            @testset "StridedMatrix * OneElementMatrix" begin
+                B = OneElement(SMatrix{3,2}(1:6), (size(A,2),2), (size(A,2),4))
+                C = [SMatrix{2,2}(1:4) for i in axes(A,1), j in axes(B,2)]
+                @test mul!(copy(C), A, B) == A * B
+                @test mul!(copy(C), A, B, 2, 2) == 2 * A * B + 2 * C
+            end
+            @testset "StridedMatrix * OneElementVector" begin
+                B = OneElement(SMatrix{3,2}(1:6), (size(A,2),), (size(A,2),))
+                C = [SMatrix{2,2}(1:4) for i in axes(A,1)]
+                @test mul!(copy(C), A, B) == A * B
+                @test mul!(copy(C), A, B, 2, 2) == 2 * A * B + 2 * C
+            end
+
+            A = OneElement(SMatrix{3,2}(1:6), (3,2), (5,4))
+            @testset "OneElementMatrix * StridedMatrix" begin
+                B = [SMatrix{2,3}(1:6)*(i+j) for i in axes(A,2), j in 1:2]
+                C = [SMatrix{3,3}(1:9) for i in axes(A,1), j in axes(B,2)]
+                @test mul!(copy(C), A, B) == A * B
+                @test mul!(copy(C), A, B, 2, 2) == 2 * A * B + 2 * C
+            end
+            @testset "OneElementMatrix * StridedVector" begin
+                B = [SMatrix{2,3}(1:6)*i for i in axes(A,2)]
+                C = [SMatrix{3,3}(1:9) for i in axes(A,1)]
+                @test mul!(copy(C), A, B) == A * B
+                @test mul!(copy(C), A, B, 2, 2) == 2 * A * B + 2 * C
+            end
+            @testset "OneElementMatrix * OneElementMatrix" begin
+                B = OneElement(SMatrix{2,3}(1:6), (2,4), (size(A,2), 3))
+                C = [SMatrix{3,3}(1:9) for i in axes(A,1), j in axes(B,2)]
+                @test mul!(copy(C), A, B) == A * B
+                @test mul!(copy(C), A, B, 2, 2) == 2 * A * B + 2 * C
+            end
+            @testset "OneElementMatrix * OneElementVector" begin
+                B = OneElement(SMatrix{2,3}(1:6), 2, size(A,2))
+                C = [SMatrix{3,3}(1:9) for i in axes(A,1)]
+                @test mul!(copy(C), A, B) == A * B
+                @test mul!(copy(C), A, B, 2, 2) == 2 * A * B + 2 * C
+            end
+        end
+        @testset "non-commutative" begin
+            A = OneElement(quat(rand(4)...), (2,3), (3,4))
+            for (B,C) in (
+                        # OneElementMatrix * OneElementVector
+                        (OneElement(quat(rand(4)...), 3, size(A,2)),
+                            [quat(rand(4)...) for i in axes(A,1)]),
+
+                        # OneElementMatrix * OneElementMatrix
+                        (OneElement(quat(rand(4)...), (3,2), (size(A,2), 4)),
+                            [quat(rand(4)...) for i in axes(A,1), j in 1:4]),
+                        )
+                @test mul!(copy(C), A, B) ≈ A * B
+                α, β = quat(0,0,1,0), quat(1,0,1,0)
+                @test mul!(copy(C), A, B, α, β) ≈ mul!(copy(C), A, Array(B), α, β) ≈ A * B * α + C * β
+            end
+
+            A = [quat(rand(4)...)*(i+j) for i in 1:2, j in 1:3]
+            for (B,C) in (
+                        # StridedMatrix * OneElementVector
+                        (OneElement(quat(rand(4)...), 1, size(A,2)),
+                            [quat(rand(4)...) for i in axes(A,1)]),
+
+                        # StridedMatrix * OneElementMatrix
+                        (OneElement(quat(rand(4)...), (2,2), (size(A,2), 4)),
+                            [quat(rand(4)...) for i in axes(A,1), j in 1:4]),
+                        )
+                @test mul!(copy(C), A, B) ≈ A * B
+                α, β = quat(0,0,1,0), quat(1,0,1,0)
+                @test mul!(copy(C), A, B, α, β) ≈ mul!(copy(C), A, Array(B), α, β) ≈ A * B * α + C * β
+            end
+
+            A = OneElement(quat(rand(4)...), (2,2), (3, 4))
+            for (B,C) in (
+                        # OneElementMatrix * StridedMatrix
+                        ([quat(rand(4)...) for i in axes(A,2), j in 1:3],
+                            [quat(rand(4)...) for i in axes(A,1), j in 1:3]),
+
+                        # OneElementMatrix * StridedVector
+                        ([quat(rand(4)...) for i in axes(A,2)],
+                            [quat(rand(4)...) for i in axes(A,1)]),
+                        )
+                @test mul!(copy(C), A, B) ≈ A * B
+                α, β = quat(0,0,1,0), quat(1,0,1,0)
+                @test mul!(copy(C), A, B, α, β) ≈ mul!(copy(C), A, Array(B), α, β) ≈ A * B * α + C * β
+            end
+        end
     end
 
     @testset "multiplication/division by a number" begin

From 545fa12ec7057e68d48690a73e8b9da3d3d969b2 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 25 Apr 2024 20:01:45 +0530
Subject: [PATCH 24/35] Fix adjtransvec multiplication with
 `AbstractZerosMatrix` (#341)

* Fix adjtransvec multiplication with zerosmatrix

* Add three-term mul tests
---
 src/fillalgebra.jl |  5 ++--
 test/runtests.jl   | 64 ++++++++++++++++++++++++++++++----------------
 2 files changed, 44 insertions(+), 25 deletions(-)

diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl
index a887eca2..928c04a7 100644
--- a/src/fillalgebra.jl
+++ b/src/fillalgebra.jl
@@ -102,9 +102,8 @@ for MT in (:(AbstractMatrix{T}), :(Transpose{<:Any, <:AbstractMatrix{T}}), :(Adj
             :(AbstractTriangular{T}))
     @eval *(a::$MT, b::AbstractZerosVector) where {T} = mult_zeros(a, b)
 end
-for MT in (:(Transpose{<:Any, <:AbstractVector}), :(Adjoint{<:Any, <:AbstractVector}))
-    @eval *(a::$MT, b::AbstractZerosMatrix) = mult_zeros(a, b)
-end
+*(a::Transpose{<:Any, <:AbstractVector}, b::AbstractZerosMatrix) = transpose(transpose(b) * parent(a))
+*(a::Adjoint{<:Any, <:AbstractVector}, b::AbstractZerosMatrix) = adjoint(adjoint(b) * parent(a))
 *(a::AbstractZerosMatrix, b::AbstractVector) = mult_zeros(a, b)
 
 function lmul_diag(a::Diagonal, b)
diff --git a/test/runtests.jl b/test/runtests.jl
index fe5fafa3..e7a5eda6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1633,8 +1633,8 @@ end
             @test transpose(A) * Zeros(mA) ≡ Zeros(nA)
             @test A' * Zeros(mA) ≡ Zeros(nA)
 
-            @test transpose(a) * Zeros(la, 3) ≡ Zeros(1,3)
-            @test a' * Zeros(la,3) ≡ Zeros(1,3)
+            @test transpose(a) * Zeros(la, 3) ≡ transpose(Zeros(3))
+            @test a' * Zeros(la,3) ≡ adjoint(Zeros(3))
 
             @test Zeros(la)' * Transpose(Adjoint(a)) == 0.0
 
@@ -1701,30 +1701,50 @@ end
     @test (1:5)'E == (1.0:5)'
     @test E*E ≡ E
 
+    n  = 10
+    k  = 12
+    m  = 15
     for T in (Float64, ComplexF64)
-        fv = T == Float64 ? Float64(1.6) : ComplexF64(1.6, 1.3)
-        n  = 10
-        k  = 12
-        m  = 15
-        fillvec = Fill(fv, k)
-        fillmat = Fill(fv, k, m)
-        A  = rand(ComplexF64, n, k)
-        @test A*fillvec ≈ A*Array(fillvec)
-        @test A*fillmat ≈ A*Array(fillmat)
-        A  = rand(ComplexF64, k, n)
-        @test transpose(A)*fillvec ≈ transpose(A)*Array(fillvec)
-        @test transpose(A)*fillmat ≈ transpose(A)*Array(fillmat)
-        @test adjoint(A)*fillvec ≈ adjoint(A)*Array(fillvec)
-        @test adjoint(A)*fillmat ≈ adjoint(A)*Array(fillmat)
-
-        # inplace C = F * B' * alpha + C * beta
+        Ank  = rand(T, n, k)
+        Akn = rand(T, k, n)
+        Ak = rand(T, k)
+        onesm = ones(m)
+        zerosm = zeros(m)
+
+        fv = T == Float64 ? T(1.6) : T(1.6, 1.3)
+
+        for (fillvec, fillmat) in ((Fill(fv, k), Fill(fv, k, m)),
+                                    (Ones(T, k), Ones(T, k, m)),
+                                    (Zeros(T, k), Zeros(T, k, m)))
+
+            Afillvec = Array(fillvec)
+            Afillmat = Array(fillmat)
+            @test Ank * fillvec ≈ Ank * Afillvec
+            @test Ank * fillmat ≈ Ank * Afillmat
+
+            for A  in (Akn, Ak)
+                @test transpose(A)*fillvec ≈ transpose(A)*Afillvec
+                AtF = transpose(A)*fillmat
+                AtM = transpose(A)*Afillmat
+                @test AtF ≈ AtM
+                @test AtF * Ones(m) ≈ AtM * onesm
+                @test AtF * Zeros(m) ≈ AtM * zerosm
+                @test adjoint(A)*fillvec ≈ adjoint(A)*Afillvec
+                AadjF = adjoint(A)*fillmat
+                AadjM = adjoint(A)*Afillmat
+                @test AadjF ≈ AadjM
+                @test AadjF * Ones(m) ≈ AadjM * onesm
+                @test AadjF * Zeros(m) ≈ AadjM * zerosm
+            end
+        end
+
+        # inplace C = F * A' * alpha + C * beta
         F = Fill(fv, m, k)
-        A = Array(F)
-        B = rand(T, n, k)
+        M = Array(F)
         C = rand(T, m, n)
         @testset for f in (adjoint, transpose)
-            @test mul!(copy(C), F, f(B)) ≈ mul!(copy(C), A, f(B))
-            @test mul!(copy(C), F, f(B), 1.0, 2.0) ≈ mul!(copy(C), A, f(B), 1.0, 2.0)
+            @test mul!(copy(C), F, f(Ank)) ≈ mul!(copy(C), M, f(Ank))
+            @test mul!(copy(C), F, f(Ank), 1.0, 2.0) ≈ mul!(copy(C), M, f(Ank), 1.0, 2.0)
         end
     end
 

From e2ec19d9b5ad6ae2eba95972d040bc37c098aa64 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Fri, 26 Apr 2024 00:07:32 +0530
Subject: [PATCH 25/35] Bump version to v1.10.2 (#362)

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index e6fff2d1..caea01ad 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.10.1"
+version = "1.10.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

From 10afd8ed6ca4aa69ffa8f77c2486c2b699cdd8c5 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sat, 27 Apr 2024 20:15:29 +0530
Subject: [PATCH 26/35] Import deps transitively in SparseArrays extension
 (#307)

* Import deps transitively in SparseArray extension

* exclude Base

* Move comment
---
 ext/FillArraysSparseArraysExt.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/ext/FillArraysSparseArraysExt.jl b/ext/FillArraysSparseArraysExt.jl
index 6a82f246..a33b89d8 100644
--- a/ext/FillArraysSparseArraysExt.jl
+++ b/ext/FillArraysSparseArraysExt.jl
@@ -1,10 +1,12 @@
 module FillArraysSparseArraysExt
 
 using SparseArrays
-import Base: convert, kron
 using FillArrays
 using FillArrays: RectDiagonalFill, RectOrDiagonalFill, ZerosVector, ZerosMatrix, getindex_value
-using LinearAlgebra
+import Base: convert, kron
+# Specifying the full namespace is necessary because of https://github.com/JuliaLang/julia/issues/48533
+# See https://github.com/JuliaStats/LogExpFunctions.jl/pull/63
+using FillArrays.LinearAlgebra
 
 ##################
 ## Sparse arrays

From 45113c3b41dcd3bb9c322415c12bb8db57325d14 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 28 Apr 2024 08:07:20 +0530
Subject: [PATCH 27/35] Function to query the non-zero index of a OneElement
 (#345)

* Function to query the non-zero index of a OneElement

* Return a CartesianIndex
---
 src/oneelement.jl | 48 +++++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl  |  3 +++
 2 files changed, 51 insertions(+)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index 28bc2c4b..b731071f 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -47,6 +47,54 @@ Base.@propagate_inbounds function Base.getindex(A::OneElement{T,N}, kj::Vararg{I
     ifelse(kj == A.ind, A.val, zero(T))
 end
 
+"""
+    nzind(A::OneElement{T,N}) -> CartesianIndex{N}
+
+Return the index where `A` contains a non-zero value.
+
+!!! note
+    The indices are not guaranteed to lie within the valid index bounds for `A`,
+    and if `FillArrays.nzind(A) ∉ CartesianIndices(A)` then `all(iszero, A)`.
+    On the other hand, if `FillArrays.nzind(A) in CartesianIndices(A)` then
+    `A[FillArrays.nzind(A)] == FillArrays.getindex_value(A)`
+
+# Examples
+```jldoctest
+julia> A = OneElement(2, (1,2), (2,2))
+2×2 OneElement{Int64, 2, Tuple{Int64, Int64}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
+ ⋅  2
+ ⋅  ⋅
+
+julia> FillArrays.nzind(A)
+CartesianIndex(1, 2)
+
+julia> A[FillArrays.nzind(A)]
+2
+```
+"""
+nzind(f::OneElement) = CartesianIndex(f.ind)
+
+"""
+    getindex_value(A::OneElement)
+
+Return the only non-zero value stored in `A`.
+
+!!! note
+    If the index at which the value is stored doesn't lie within the valid indices of `A`, then
+    this returns `zero(eltype(A))`.
+
+# Examples
+```jldoctest
+julia> A = OneElement(2, 3)
+3-element OneElement{Int64, 1, Tuple{Int64}, Tuple{Base.OneTo{Int64}}}:
+ ⋅
+ 1
+ ⋅
+
+julia> FillArrays.getindex_value(A)
+1
+```
+"""
 getindex_value(A::OneElement) = all(in.(A.ind, axes(A))) ? A.val : zero(eltype(A))
 
 Base.AbstractArray{T,N}(A::OneElement{<:Any,N}) where {T,N} = OneElement(T(A.val), A.ind, A.axes)
diff --git a/test/runtests.jl b/test/runtests.jl
index e7a5eda6..2fd3c5ca 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2114,11 +2114,13 @@ end
 
 @testset "OneElement" begin
     A = OneElement(2, (), ())
+    @test FillArrays.nzind(A) == CartesianIndex()
     @test A == Fill(2, ())
     @test A[] === 2
 
     e₁ = OneElement(2, 5)
     @test e₁ == [0,1,0,0,0]
+    @test FillArrays.nzind(e₁) == CartesianIndex(2)
     @test_throws BoundsError e₁[6]
 
     f₁ = AbstractArray{Float64}(e₁)
@@ -2138,6 +2140,7 @@ end
 
     V = OneElement(2, (2,3), (3,4))
     @test V == [0 0 0 0; 0 0 2 0; 0 0 0 0]
+    @test FillArrays.nzind(V) == CartesianIndex(2,3)
 
     Vf = AbstractArray{Float64}(V)
     @test Vf isa OneElement{Float64,2}

From dfe0e62452e8c4281f8b8cd362d3cf3c1c522ecf Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sat, 4 May 2024 21:38:24 +0530
Subject: [PATCH 28/35] Specialize isassigned for AbstractFill (#351)

* Specialize isassigned for AbstractFill

* Define method only for Integer indices

* isassigned for OneElement
---
 src/FillArrays.jl |  4 ++++
 src/oneelement.jl |  5 +++++
 test/runtests.jl  | 17 +++++++++++++++++
 3 files changed, 26 insertions(+)

diff --git a/src/FillArrays.jl b/src/FillArrays.jl
index 581e8252..d8db9f5d 100644
--- a/src/FillArrays.jl
+++ b/src/FillArrays.jl
@@ -33,6 +33,10 @@ const AbstractFillVecOrMat{T} = Union{AbstractFillVector{T},AbstractFillMatrix{T
 
 ==(a::AbstractFill, b::AbstractFill) = axes(a) == axes(b) && getindex_value(a) == getindex_value(b)
 
+@inline function Base.isassigned(F::AbstractFill, i::Integer...)
+    @boundscheck checkbounds(Bool, F, to_indices(F, i)...) || return false
+    return true
+end
 
 @inline function _fill_getindex(F::AbstractFill, kj::Integer...)
     @boundscheck checkbounds(F, kj...)
diff --git a/src/oneelement.jl b/src/oneelement.jl
index b731071f..f6204abb 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -97,6 +97,11 @@ julia> FillArrays.getindex_value(A)
 """
 getindex_value(A::OneElement) = all(in.(A.ind, axes(A))) ? A.val : zero(eltype(A))
 
+@inline function Base.isassigned(F::OneElement, i::Integer...)
+    @boundscheck checkbounds(Bool, F, to_indices(F, i)...) || return false
+    return true
+end
+
 Base.AbstractArray{T,N}(A::OneElement{<:Any,N}) where {T,N} = OneElement(T(A.val), A.ind, A.axes)
 
 Base.replace_in_print_matrix(o::OneElementVector, k::Integer, j::Integer, s::AbstractString) =
diff --git a/test/runtests.jl b/test/runtests.jl
index 2fd3c5ca..8af5ed45 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -261,6 +261,15 @@ end
     end
 end
 
+@testset "isassigned" begin
+    for f in (Fill("", 3, 4), Zeros(3,4), Ones(3,4))
+        @test !isassigned(f, 0, 0)
+        @test isassigned(f, 2, 2)
+        @test !isassigned(f, 10, 10)
+        @test_throws ArgumentError isassigned(f, true)
+    end
+end
+
 @testset "indexing" begin
     A = Fill(3.0,5)
     @test A[1:3] ≡ Fill(3.0,3)
@@ -2184,6 +2193,14 @@ end
         @test adjoint(A) == OneElement(3, (1,2), (1,4))
     end
 
+    @testset "isassigned" begin
+        f = OneElement(2, (3,3), (4,4))
+        @test !isassigned(f, 0, 0)
+        @test isassigned(f, 2, 2)
+        @test !isassigned(f, 10, 10)
+        @test_throws ArgumentError isassigned(f, true)
+    end
+
     @testset "matmul" begin
         A = reshape(Float64[1:9;], 3, 3)
         v = reshape(Float64[1:3;], 3)

From 39df7ec50dac0170e7ef6dd82c7fd4b8a9f7f7cd Mon Sep 17 00:00:00 2001
From: Sheehan Olver <solver@mac.com>
Date: Sat, 4 May 2024 18:18:01 +0100
Subject: [PATCH 29/35] Overload accumulate and make types of cumsum consistent
 with Vector (#363)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Support accumulate(±, ::AbstractFill)

* Overload accumulate and make types of cumsum consistent with Vector
---
 Project.toml         |  2 +-
 src/FillArrays.jl    | 28 +++++++++++++++++++----
 src/fillbroadcast.jl |  2 ++
 test/runtests.jl     | 54 +++++++++++++++++++++++++++++++-------------
 4 files changed, 64 insertions(+), 22 deletions(-)

diff --git a/Project.toml b/Project.toml
index caea01ad..dd70d301 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "FillArrays"
 uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
-version = "1.10.2"
+version = "1.11"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/src/FillArrays.jl b/src/FillArrays.jl
index d8db9f5d..4671053a 100644
--- a/src/FillArrays.jl
+++ b/src/FillArrays.jl
@@ -7,7 +7,7 @@ import Base: size, getindex, setindex!, IndexStyle, checkbounds, convert,
     any, all, axes, isone, iszero, iterate, unique, allunique, permutedims, inv,
     copy, vec, setindex!, count, ==, reshape, map, zero,
     show, view, in, mapreduce, one, reverse, promote_op, promote_rule, repeat,
-    parent, similar, issorted
+    parent, similar, issorted, add_sum, accumulate, OneTo
 
 import LinearAlgebra: rank, svdvals!, tril, triu, tril!, triu!, diag, transpose, adjoint, fill!,
     dot, norm2, norm1, normInf, normMinusInf, normp, lmul!, rmul!, diagzero, AdjointAbsVec, TransposeAbsVec,
@@ -576,14 +576,32 @@ sum(x::AbstractZeros) = getindex_value(x)
 
 # needed to support infinite case
 steprangelen(st...) = StepRangeLen(st...)
-cumsum(x::AbstractFill{<:Any,1}) = steprangelen(getindex_value(x), getindex_value(x), length(x))
+function cumsum(x::AbstractFill{T,1}) where T
+    V = promote_op(add_sum, T, T)
+    steprangelen(convert(V,getindex_value(x)), getindex_value(x), length(x))
+end
 
-cumsum(x::AbstractZerosVector) = x
-cumsum(x::AbstractZerosVector{Bool}) = x
-cumsum(x::AbstractOnesVector{II}) where II<:Integer = convert(AbstractVector{II}, oneto(length(x)))
+cumsum(x::AbstractZerosVector{T}) where T = _range_convert(AbstractVector{promote_op(add_sum, T, T)}, x)
+cumsum(x::AbstractZerosVector{Bool}) = _range_convert(AbstractVector{Int}, x)
+cumsum(x::AbstractOnesVector{T}) where T<:Integer = _range_convert(AbstractVector{promote_op(add_sum, T, T)}, oneto(length(x)))
 cumsum(x::AbstractOnesVector{Bool}) = oneto(length(x))
 
 
+for op in (:+, :-)
+    @eval begin
+        function accumulate(::typeof($op), x::AbstractFill{T,1}) where T
+            V = promote_op($op, T, T)
+            steprangelen(convert(V,getindex_value(x)), $op(getindex_value(x)), length(x))
+        end
+
+        accumulate(::typeof($op), x::AbstractZerosVector{T}) where T = _range_convert(AbstractVector{promote_op($op, T, T)}, x)
+        accumulate(::typeof($op), x::AbstractZerosVector{Bool}) = _range_convert(AbstractVector{Int}, x)
+    end
+end
+
+accumulate(::typeof(+), x::AbstractOnesVector{T}) where T<:Integer = _range_convert(AbstractVector{promote_op(+, T, T)}, oneto(length(x)))
+accumulate(::typeof(+), x::AbstractOnesVector{Bool}) = oneto(length(x))
+
 #########
 # Diff
 #########
diff --git a/src/fillbroadcast.jl b/src/fillbroadcast.jl
index d286418a..2b5ea59c 100644
--- a/src/fillbroadcast.jl
+++ b/src/fillbroadcast.jl
@@ -177,7 +177,9 @@ end
 # special case due to missing converts for ranges
 _range_convert(::Type{AbstractVector{T}}, a::AbstractRange{T}) where T = a
 _range_convert(::Type{AbstractVector{T}}, a::AbstractUnitRange) where T = convert(T,first(a)):convert(T,last(a))
+_range_convert(::Type{AbstractVector{T}}, a::OneTo) where T = OneTo(convert(T, a.stop))
 _range_convert(::Type{AbstractVector{T}}, a::AbstractRange) where T = convert(T,first(a)):step(a):convert(T,last(a))
+_range_convert(::Type{AbstractVector{T}}, a::ZerosVector) where T = ZerosVector{T}(length(a))
 
 
 # TODO: replacing with the following will support more general broadcasting.
diff --git a/test/runtests.jl b/test/runtests.jl
index 8af5ed45..2e2aba22 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -826,31 +826,53 @@ end
     @test_throws MethodError sort!(Fill(im, 2))
 end
 
-@testset "Cumsum and diff" begin
-    @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)) ≡ StepRangeLen(3,3,10)
-
-    @test sum(Ones(10)) ≡ 10.0
-    @test sum(x -> x + 1, Ones(10)) ≡ 20.0
-    @test cumsum(Ones(10)) ≡ StepRangeLen(1.0, 1.0, 10)
+@testset "Cumsum, accumulate and diff" begin
+    @test @inferred(sum(Fill(3,10))) ≡ 30
+    @test @inferred(reduce(+, Fill(3,10))) ≡ 30
+    @test @inferred(sum(x -> x + 1, Fill(3,10))) ≡ 40
+    @test @inferred(cumsum(Fill(3,10))) ≡ @inferred(accumulate(+, Fill(3,10))) ≡ StepRangeLen(3,3,10)
+    @test @inferred(accumulate(-, Fill(3,10))) ≡ StepRangeLen(3,-3,10)
+
+    @test @inferred(sum(Ones(10))) ≡ 10.0
+    @test @inferred(sum(x -> x + 1, Ones(10))) ≡ 20.0
+    @test @inferred(cumsum(Ones(10))) ≡ @inferred(accumulate(+, Ones(10))) ≡ StepRangeLen(1.0, 1.0, 10)
+    @test @inferred(accumulate(-, Ones(10))) ≡ StepRangeLen(1.0,-1.0,10)
 
     @test sum(Ones{Int}(10)) ≡ 10
     @test sum(x -> x + 1, Ones{Int}(10)) ≡ 20
-    @test cumsum(Ones{Int}(10)) ≡ Base.OneTo(10)
+    @test cumsum(Ones{Int}(10)) ≡ accumulate(+,Ones{Int}(10)) ≡ Base.OneTo(10)
+    @test accumulate(-, Ones{Int}(10)) ≡ StepRangeLen(1,-1,10)
 
     @test sum(Zeros(10)) ≡ 0.0
     @test sum(x -> x + 1, Zeros(10)) ≡ 10.0
-    @test cumsum(Zeros(10)) ≡ Zeros(10)
+    @test cumsum(Zeros(10)) ≡ accumulate(+,Zeros(10)) ≡ accumulate(-,Zeros(10)) ≡ Zeros(10)
 
     @test sum(Zeros{Int}(10)) ≡ 0
     @test sum(x -> x + 1, Zeros{Int}(10)) ≡ 10
-    @test cumsum(Zeros{Int}(10)) ≡ Zeros{Int}(10)
-
-    @test cumsum(Zeros{Bool}(10)) ≡ Zeros{Bool}(10)
-    @test cumsum(Ones{Bool}(10)) ≡ Base.OneTo{Int}(10)
-    @test cumsum(Fill(true,10)) ≡ StepRangeLen(true, true, 10)
+    @test cumsum(Zeros{Int}(10)) ≡ accumulate(+,Zeros{Int}(10)) ≡ accumulate(-,Zeros{Int}(10)) ≡ Zeros{Int}(10)
+
+    # we want cumsum of fills to match the types of the standard cusum
+    @test all(cumsum(Zeros{Bool}(10)) .≡ cumsum(zeros(Bool,10)))
+    @test all(accumulate(+, Zeros{Bool}(10)) .≡ accumulate(+, zeros(Bool,10)) .≡ accumulate(-, zeros(Bool,10)))
+    @test cumsum(Zeros{Bool}(10)) ≡ accumulate(+, Zeros{Bool}(10)) ≡ accumulate(-, Zeros{Bool}(10)) ≡ Zeros{Int}(10)
+    @test cumsum(Ones{Bool}(10)) ≡ accumulate(+, Ones{Bool}(10)) ≡ Base.OneTo{Int}(10)
+    @test all(cumsum(Fill(true,10)) .≡ cumsum(fill(true,10)))
+    @test cumsum(Fill(true,10)) ≡ StepRangeLen(1, true, 10)
+
+    @test all(cumsum(Zeros{UInt8}(10)) .≡ cumsum(zeros(UInt8,10)))
+    @test all(accumulate(+, Zeros{UInt8}(10)) .≡ accumulate(+, zeros(UInt8,10)))
+    @test cumsum(Zeros{UInt8}(10)) ≡ Zeros{UInt64}(10)
+    @test accumulate(+, Zeros{UInt8}(10)) ≡ accumulate(-, Zeros{UInt8}(10)) ≡ Zeros{UInt8}(10)
+    
+    @test all(cumsum(Ones{UInt8}(10)) .≡ cumsum(ones(UInt8,10)))
+    @test all(accumulate(+, Ones{UInt8}(10)) .≡ accumulate(+, ones(UInt8,10)))
+    @test cumsum(Ones{UInt8}(10)) ≡ Base.OneTo(UInt64(10))
+    @test accumulate(+, Ones{UInt8}(10)) ≡ Base.OneTo(UInt8(10))
+    
+    @test all(cumsum(Fill(UInt8(2),10)) .≡ cumsum(fill(UInt8(2),10)))
+    @test all(accumulate(+,  Fill(UInt8(2))) .≡ accumulate(+, fill(UInt8(2))))
+    @test cumsum(Fill(UInt8(2),10)) ≡ StepRangeLen(UInt64(2), UInt8(2), 10)
+    @test accumulate(+, Fill(UInt8(2),10)) ≡ StepRangeLen(UInt8(2), UInt8(2), 10)
 
     @test diff(Fill(1,10)) ≡ Zeros{Int}(9)
     @test diff(Ones{Float64}(10)) ≡ Zeros{Float64}(9)

From 5d88772fd0a886901a1ecc7caeaa020a34c04feb Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Mon, 6 May 2024 11:15:11 +0530
Subject: [PATCH 30/35] Bump julia-actions/cache from 1 to 2 (#364)

Bumps [julia-actions/cache](https://github.com/julia-actions/cache) from 1 to 2.
- [Release notes](https://github.com/julia-actions/cache/releases)
- [Commits](https://github.com/julia-actions/cache/compare/v1...v2)

---
updated-dependencies:
- dependency-name: julia-actions/cache
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index bf59dbfb..18c5f197 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -52,7 +52,7 @@ jobs:
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1

From 4f8a966e931208a8a5aa56c909cc5e579044e421 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 9 May 2024 10:17:06 +0530
Subject: [PATCH 31/35] Reshape for a OneElement (#365)

* Reshape for a OneElement

* Tests for 0D
---
 src/oneelement.jl | 10 ++++++++++
 test/runtests.jl  | 12 ++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index f6204abb..c72bcd2d 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -372,6 +372,16 @@ function broadcasted(::DefaultArrayStyle{N}, ::typeof(\), x::Number, r::OneEleme
     OneElement(x \ r.val, r.ind, axes(r))
 end
 
+# reshape
+
+function Base.reshape(A::OneElement, shape::Tuple{Vararg{Int}})
+    prod(shape) == length(A) || throw(DimensionMismatch("new dimension $shape must be consistent with array size $(length(A))"))
+    # we use the fact that the linear index of the non-zero value is preserved
+    oldlinind = LinearIndices(A)[A.ind...]
+    newcartind = CartesianIndices(shape)[oldlinind]
+    OneElement(A.val, Tuple(newcartind), shape)
+end
+
 # show
 _maybesize(t::Tuple{Base.OneTo{Int}, Vararg{Base.OneTo{Int}}}) = size.(t,1)
 _maybesize(t) = t
diff --git a/test/runtests.jl b/test/runtests.jl
index 2e2aba22..d4b2bb95 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2215,6 +2215,18 @@ end
         @test adjoint(A) == OneElement(3, (1,2), (1,4))
     end
 
+    @testset "reshape" begin
+        for O in (OneElement(2, (2,3), (4,5)), OneElement(2, (2,), (20,)),
+                    OneElement(2, (1,2,2), (2,2,5)))
+            A = Array(O)
+            for shp in ((2,5,2), (5,1,4), (20,), (2,2,5,1,1))
+                @test reshape(O, shp) == reshape(A, shp)
+            end
+        end
+        O = OneElement(2, (), ())
+        @test reshape(O, ()) === O
+    end
+
     @testset "isassigned" begin
         f = OneElement(2, (3,3), (4,4))
         @test !isassigned(f, 0, 0)

From 52b7de2dcae4df240455c9b747c07d85cebf3f85 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Fri, 28 Jun 2024 18:44:57 +0530
Subject: [PATCH 32/35] Use "lts" version specifier in CI script (#366)

---
 .github/workflows/ci.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 18c5f197..8e55678d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -38,7 +38,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - 'lts'
           - '1'
         os:
           - ubuntu-latest
@@ -48,7 +48,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@v2.2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

From 4bfbeecdab1c1f40959b3b2e2b268816e1c75b9d Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 9 Jul 2024 01:14:30 +0530
Subject: [PATCH 33/35] tril/triu for ZerosMatrix and OneElementMatrix (#334)

---
 src/fillalgebra.jl |  3 +++
 src/oneelement.jl  | 13 +++++++++++++
 test/runtests.jl   | 19 +++++++++++++++++++
 3 files changed, 35 insertions(+)

diff --git a/src/fillalgebra.jl b/src/fillalgebra.jl
index 928c04a7..c17e46cd 100644
--- a/src/fillalgebra.jl
+++ b/src/fillalgebra.jl
@@ -551,3 +551,6 @@ end
 function LinearAlgebra.istril(A::AbstractFillMatrix, k::Integer = 0)
     iszero(A) || k >= size(A,2)-1
 end
+
+triu(A::AbstractZerosMatrix, k::Integer=0) = A
+tril(A::AbstractZerosMatrix, k::Integer=0) = A
diff --git a/src/oneelement.jl b/src/oneelement.jl
index c72bcd2d..7cffbee4 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -349,7 +349,20 @@ end
 adjoint(A::OneElementMatrix) = OneElement(adjoint(A.val), reverse(A.ind), reverse(A.axes))
 transpose(A::OneElementMatrix) = OneElement(transpose(A.val), reverse(A.ind), reverse(A.axes))
 
+# tril/triu
+
+function tril(A::OneElementMatrix, k::Integer=0)
+    nzband = A.ind[2] - A.ind[1]
+    OneElement(nzband > k ? zero(A.val) : A.val, A.ind, axes(A))
+end
+
+function triu(A::OneElementMatrix, k::Integer=0)
+    nzband = A.ind[2] - A.ind[1]
+    OneElement(nzband < k ? zero(A.val) : A.val, A.ind, axes(A))
+end
+
 # broadcast
+
 function broadcasted(::DefaultArrayStyle{N}, ::typeof(conj), r::OneElement{<:Any,N}) where {N}
     OneElement(conj(r.val), r.ind, axes(r))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index d4b2bb95..1e066665 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2501,6 +2501,17 @@ end
         @test A / (2 + 3.0im) === OneElement(val / (2 + 3.0im), (2,2), (3,4)) == B / (2 + 3.0im)
     end
 
+
+    @testset "tril/triu" begin
+        for A in (OneElement(3, (4,2), (4,5)), OneElement(3, (2,3), (4,5)), OneElement(3, (3,3), (4,5)))
+            B = Array(A)
+            for k in -5:5
+                @test tril(A, k) == tril(B, k)
+                @test triu(A, k) == triu(B, k)
+            end
+        end
+    end
+
     @testset "broadcasting" begin
         for v in (OneElement(2, 3, 4), OneElement(2im, (1,2), (3,4)))
             w = Array(v)
@@ -2688,3 +2699,11 @@ end
         end
     end
 end
+
+@testset "triu/tril for Zeros" begin
+    Z = Zeros(3, 4)
+    @test triu(Z) === Z
+    @test tril(Z) === Z
+    @test triu(Z, 2) === Z
+    @test tril(Z, 2) === Z
+end

From 1896e307489eed06cb09cb67c4d5abdfc48529ee Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Tue, 9 Jul 2024 14:07:47 +0530
Subject: [PATCH 34/35] zero for OneElement (#332)

* zero for OneElement

* Specialize iszero

* Specialize isone

* Return a OneElement instead of a Zero
---
 src/oneelement.jl | 11 +++++++++++
 test/runtests.jl  | 22 ++++++++++++++++++++++
 2 files changed, 33 insertions(+)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index 7cffbee4..7fc63594 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -115,6 +115,17 @@ Base.@propagate_inbounds function Base.setindex(A::AbstractZeros{T,N}, v, kj::Va
     OneElement(convert(T, v), kj, axes(A))
 end
 
+zero(A::OneElement) = OneElement(zero(A.val), A.ind, A.axes)
+
+iszero(A::OneElement) = iszero(getindex_value(A))
+
+function isone(A::OneElementMatrix)
+    lenA = length(A)
+    lenA == 0 && return true
+    lenA > 1 && return false
+    isone(getindex_value(A))
+end
+
 *(x::OneElement, b::Number) = OneElement(x.val * b, x.ind, x.axes)
 *(b::Number, x::OneElement) = OneElement(b * x.val, x.ind, x.axes)
 /(x::OneElement, b::Number) = OneElement(x.val / b, x.ind, x.axes)
diff --git a/test/runtests.jl b/test/runtests.jl
index 1e066665..08bee6ac 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2501,6 +2501,28 @@ end
         @test A / (2 + 3.0im) === OneElement(val / (2 + 3.0im), (2,2), (3,4)) == B / (2 + 3.0im)
     end
 
+    @testset "zero/iszero" begin
+        v = OneElement(10, 3, 4)
+        @test v + zero(v) == v
+        @test typeof(zero(v)) == typeof(v)
+
+        @test !iszero(v)
+        @test iszero(OneElement(0, 3, 4))
+        @test iszero(OneElement(0, 5, 4))
+        @test iszero(OneElement(3, (2,2), (0,0)))
+        @test iszero(OneElement(3, (2,2), (1,2)))
+
+        v = OneElement(SMatrix{2,2}(1:4), 3, 4)
+        @test v + zero(v) == v
+        @test typeof(zero(v)) == typeof(v)
+    end
+
+    @testset "isone" begin
+        @test isone(OneElement(3, (0,0), (0,0)))
+        @test isone(OneElement(1, (1,1), (1,1)))
+        @test !isone(OneElement(2, (1,1), (1,1)))
+        @test !isone(OneElement(1, (2,2), (2,2)))
+    end
 
     @testset "tril/triu" begin
         for A in (OneElement(3, (4,2), (4,5)), OneElement(3, (2,3), (4,5)), OneElement(3, (3,3), (4,5)))

From a7ddcca0ac7b5a76be09864e6a3fed5fffca820f Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 21 Jul 2024 12:47:41 +0530
Subject: [PATCH 35/35] Specialize isbanded for OneElement (#333)

---
 src/oneelement.jl |  5 +++++
 test/runtests.jl  | 31 +++++++++++++++++++++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/src/oneelement.jl b/src/oneelement.jl
index 7fc63594..f3325547 100644
--- a/src/oneelement.jl
+++ b/src/oneelement.jl
@@ -360,6 +360,11 @@ end
 adjoint(A::OneElementMatrix) = OneElement(adjoint(A.val), reverse(A.ind), reverse(A.axes))
 transpose(A::OneElementMatrix) = OneElement(transpose(A.val), reverse(A.ind), reverse(A.axes))
 
+# isbanded
+function LinearAlgebra.isbanded(A::OneElementMatrix, kl::Integer, ku::Integer)
+    iszero(getindex_value(A)) || kl <= A.ind[2] - A.ind[1] <= ku
+end
+
 # tril/triu
 
 function tril(A::OneElementMatrix, k::Integer=0)
diff --git a/test/runtests.jl b/test/runtests.jl
index 08bee6ac..e8e79a1d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2501,6 +2501,37 @@ end
         @test A / (2 + 3.0im) === OneElement(val / (2 + 3.0im), (2,2), (3,4)) == B / (2 + 3.0im)
     end
 
+
+    @testset "isbanded" begin
+        A = OneElement(3, (2,3), (4,5))
+        @test !isdiag(A)
+        @test istriu(A)
+        @test !istril(A)
+        @test LinearAlgebra.isbanded(A, 1, 2)
+        @test LinearAlgebra.isbanded(A, -1, 2)
+        @test !LinearAlgebra.isbanded(A, 2, 2)
+
+        A = OneElement(3, (4,2), (4,5))
+        @test !isdiag(A)
+        @test !istriu(A)
+        @test istril(A)
+        @test LinearAlgebra.isbanded(A, -2, -2)
+        @test LinearAlgebra.isbanded(A, -2, 2)
+        @test !LinearAlgebra.isbanded(A, 2, 2)
+
+        for A in (OneElement(3, (2,2), (4,5)), OneElement(3, (1,1), (1,2)), OneElement(3, (8,7), (2,1)))
+            @test isdiag(A)
+            @test istriu(A)
+            @test istril(A)
+        end
+
+        for A in (OneElement(0, (2,3), (4,5)), OneElement(0, (3,2), (4,5)))
+            @test isdiag(A)
+            @test istriu(A)
+            @test istril(A)
+        end
+    end
+
     @testset "zero/iszero" begin
         v = OneElement(10, 3, 4)
         @test v + zero(v) == v