From 73b82885a8ab10763cd48e1067a6b651342fd890 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Wed, 9 Sep 2026 09:06:55 -0400
Subject: [PATCH 1/7] Return sparse results from reductions along a dimension
 of a sparse matrix
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Fixes #43. `sum(A; dims)` and the other dimensional reductions of a
`SparseMatrixCSC` returned a dense `Matrix`, which for a hypersparse
matrix costs O(size) time and memory for a result with a handful of
entries. They now return a `SparseMatrixCSC` of the reduced shape that
stores an entry only for the rows or columns that store one themselves,
unless the reduction of a structurally empty slice is nonzero, in which
case the result is fully stored. Results without a `zero`, such as the
tuples of `extrema`, stay dense.

`reducedim_initarray` provides a structurally empty destination when the
initial value is zero and a fully stored one otherwise. A new
`_mapreducedim!` for a sparse destination reduces a fully stored one as
the dense array its values form, fills an empty one without touching the
slices that store nothing, and sends anything in between through the
element-wise kernel. Row reductions build the `1 x n` result column by
column. Column reductions use the result's value vector as a dense
workspace compressed in place when there are enough stored entries, and
otherwise sort the stored entries by row so that only rows storing
something are visited.

Measured on nightly against main, min of 9: `sum(A; dims=2)` for a
10^7 x 150 matrix with 151 entries goes from 6.5 ms and 156 MiB to
2.3 µs and 16 KiB; for 10^4 x 10^4 at 1e-3 and 1e-2 density the
reductions are within noise, with the result's extra index vectors as
the only added allocation.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01V6EdE4F3CCGE3vxr8gQYKf
---
 src/sparsematrix.jl      | 133 +++++++++++++++++++++++++++++++++++++--
 test/higherorderfns.jl   |   5 +-
 test/sparsematrix_ops.jl |  74 ++++++++++++++++++++++
 3 files changed, 206 insertions(+), 6 deletions(-)

diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index bc6ddc26..e8ef8d2e 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -2558,11 +2558,20 @@ Base.isequal(A::Transpose{<:Any,<:SparseMatrixCSCMaybeAdjOrTrans}, B::SparseMatr
 
 ## Reductions
 
-# In general, output of sparse matrix reductions will not be sparse,
-# and computing reductions along columns into SparseMatrixCSC is
-# non-trivial, so use Arrays for output. Array element type is given by `R`.
-function Base.reducedim_initarray(A::AbstractSparseMatrixCSC, region, v0, ::Type{R}) where {R}
-    fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
+# Reductions along a dimension of a sparse matrix return a sparse matrix (issue #43): the
+# result only stores an entry for the rows or columns that store one themselves, unless
+# the reduction of a structurally empty slice is nonzero. The initial array is therefore
+# structurally empty when the initial value is zero, and fully stored otherwise.
+function Base.reducedim_initarray(A::AbstractSparseMatrixCSC{<:Any,Ti}, region, v0, ::Type{R}) where {R,Ti}
+    m, n = Base.to_shape(Base.reduced_indices(A, region))
+    if !applicable(zero, R)
+        # no structural zero for the result, e.g. the tuples of `extrema`: reduce densely
+        return fill!(Array{R}(undef, m, n), v0)
+    elseif isequal(v0, zero(R))
+        return spzeros(R, Ti, m, n)
+    else
+        return SparseMatrixCSC(m, n, Ti[1 + m*j for j in 0:n], repeat(Ti.(1:m), n), fill!(Vector{R}(undef, m*n), v0))
+    end
 end
 
 # General mapreduce
@@ -2628,6 +2637,120 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas
     end
 end
 
+# Reduction of a sparse matrix into a sparse destination. A fully stored destination is
+# reduced into as the dense array its stored values form; a structurally empty one is
+# filled without touching the slices that store nothing, so the cost stays proportional to
+# the stored entries plus the length of the result; anything in between is rare and goes
+# through the element-wise kernel below.
+function Base._mapreducedim!(f::F, op::G, R::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC{T}) where {F,G,T}
+    require_one_based_indexing(A, R)
+    Base.check_reducedims(R, A)
+    isempty(A) && return R
+    if nnz(R) == length(R)
+        Base._mapreducedim!(f, op, reshape(view(nonzeros(R), 1:nnz(R)), size(R)), A)
+    elseif nnz(R) != 0
+        invoke(Base._mapreducedim!, Tuple{F,G,AbstractArray,AbstractSparseMatrixCSC{T}}, f, op, R, A)
+    elseif size(R) == (1, 1)
+        R[1, 1] = op(zero(eltype(R)), mapreduce(f, op, A))
+    elseif size(R, 1) == 1
+        _mapreducerows_sparse!(f, op, R, A)
+    elseif size(R, 2) == 1
+        _mapreducecols_sparse!(f, op, R, A)
+    else
+        # reduction over a dimension beyond 2: `R` has the shape of `A`
+        copyto!(R, op.(zero(eltype(R)), f.(A)))
+    end
+    return R
+end
+
+# `R` is a structurally empty `1 x n` sparse matrix: its columns are built in order
+function _mapreducerows_sparse!(f, op, R::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC{T}) where T
+    colptr = getcolptr(A)
+    nzval = nonzeros(A)
+    m, n = size(A)
+    z = zero(eltype(R))
+    # the reduction of a column that stores nothing, stored only when it is nonzero
+    zempty = m == 0 ? z : _mapreducezeros(f, op, T, m, z)
+    store_empty = !isequal(zempty, z)
+    Rcolptr, Rrowval, Rnzval = getcolptr(R), rowvals(R), nonzeros(R)
+    nstored = store_empty ? n : count(col -> colptr[col+1] > colptr[col], 1:n)
+    resize!(Rrowval, nstored)
+    fill!(Rrowval, 1)
+    resize!(Rnzval, nstored)
+    k = 0
+    @inbounds for col in 1:n
+        rng = colptr[col]:colptr[col+1]-1
+        if isempty(rng)
+            store_empty || (Rcolptr[col+1] = k + 1; continue)
+            v = zempty
+        else
+            r = z
+            @simd for j in rng
+                r = op(r, f(nzval[j]))
+            end
+            v = _mapreducezeros(f, op, T, m - length(rng), r)
+        end
+        k += 1
+        Rnzval[k] = v
+        Rcolptr[col+1] = k + 1
+    end
+    return R
+end
+
+# `R` is a structurally empty `m x 1` sparse matrix. With enough stored entries its value
+# vector serves as a dense workspace of length `m` that is then compressed in place; a
+# hypersparse `A` instead has its stored entries sorted by row so that only the rows storing
+# something are ever visited.
+function _mapreducecols_sparse!(f, op, R::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC{T}) where T
+    m, n = size(A)
+    z = zero(eltype(R))
+    zempty = n == 0 ? z : _mapreducezeros(f, op, T, n, z)
+    store_empty = !isequal(zempty, z)
+    Rcolptr, Rrowval, Rnzval = getcolptr(R), rowvals(R), nonzeros(R)
+    resize!(Rrowval, 0)
+    resize!(Rnzval, 0)
+    if store_empty || 8 * nnz(A) >= m
+        resize!(Rnzval, m)
+        fill!(Rnzval, z)
+        _mapreducecols!(f, op, Rnzval, A)
+        resize!(Rrowval, m)
+        if store_empty
+            Rrowval .= 1:m
+        else
+            k = 0
+            @inbounds for i in 1:m
+                w = Rnzval[i]
+                if !isequal(w, z)
+                    k += 1
+                    Rrowval[k] = i
+                    Rnzval[k] = w
+                end
+            end
+            resize!(Rrowval, k)
+            resize!(Rnzval, k)
+        end
+    else
+        rows = view(rowvals(A), 1:nnz(A))
+        vals = view(nonzeros(A), 1:nnz(A))
+        perm = sortperm(rows; alg=Base.Sort.DEFAULT_STABLE)   # keeps each row's entries in column order
+        s = 1
+        @inbounds while s <= length(perm)
+            row = rows[perm[s]]
+            r = op(z, f(vals[perm[s]]))
+            t = s + 1
+            while t <= length(perm) && rows[perm[t]] == row
+                r = op(r, f(vals[perm[t]]))
+                t += 1
+            end
+            push!(Rrowval, row)
+            push!(Rnzval, _mapreducezeros(f, op, T, n - (t - s), r))
+            s = t
+        end
+    end
+    Rcolptr[2] = length(Rnzval) + 1
+    return R
+end
+
 # General mapreducedim
 function _mapreducerows!(f, op, R::AbstractArray, A::AbstractSparseMatrixCSC{T}) where T
     require_one_based_indexing(A, R)
diff --git a/test/higherorderfns.jl b/test/higherorderfns.jl
index 66483d72..a70ca758 100644
--- a/test/higherorderfns.jl
+++ b/test/higherorderfns.jl
@@ -702,7 +702,10 @@ end
 end
 
 @testset "Issue #27836" begin
-    @test minimum(sparse([1, 2], [1, 2], ones(Int32, 2)), dims = 1) isa Matrix
+    # the result keeps the reduced eltype (and is sparse, see #43)
+    r = minimum(sparse([1, 2], [1, 2], ones(Int32, 2)), dims = 1)
+    @test r isa SparseMatrixCSC{Int32}
+    @test r == [0 0]
 end
 
 @testset "Issue #30118" begin
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index 717cba9c..d844b5f0 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -614,6 +614,80 @@ end
     @test B ≈ mapreduce(identity, +, Matrix(A), dims=2)
 end
 
+@testset "reductions along a dimension are sparse (#43)" begin
+    # every reduction with a `zero` for its result returns a sparse matrix of the reduced
+    # shape that agrees with the dense result
+    reductions = (
+        (X; dims) -> sum(X; dims), (X; dims) -> prod(X; dims),
+        (X; dims) -> maximum(X; dims), (X; dims) -> minimum(X; dims),
+        (X; dims) -> sum(abs2, X; dims), (X; dims) -> maximum(abs, X; dims),
+        (X; dims) -> count(x -> x > 0.5, X; dims),
+        (X; dims) -> any(x -> x > 0.5, X; dims), (X; dims) -> all(x -> x >= 0, X; dims),
+        (X; dims) -> mapreduce(x -> x + 1, +, X; dims),   # f(0) != 0: dense result
+        (X; dims) -> prod(x -> x + 1, X; dims),
+        (X; dims) -> sum(X; dims, init = 2.5),
+    )
+    @testset "size = ($m, $n), density = $d" for (m, n) in ((6, 5), (1, 1), (1, 9), (9, 1), (30, 20)),
+                                                 d in (0.0, 0.2, 1.0)
+        A = sprand(m, n, d)
+        M = Matrix(A)
+        for dims in (1, 2, (1, 2), 3), f in reductions
+            rs = f(A; dims)
+            rd = f(M; dims)
+            @test rs isa SparseMatrixCSC
+            @test size(rs) == size(rd)
+            @test Matrix(rs) ≈ rd
+        end
+    end
+    # only rows and columns that store something get an entry, unless the reduction of a
+    # structurally empty slice is nonzero
+    A = sparse([1, 2], [1, 1], [-1.0, 1.0], 4, 3)
+    @test nnz(sum(A; dims = 1)) == 1   # a stored, cancelled zero
+    @test nnz(sum(A; dims = 2)) == 2
+    @test nnz(prod(A; dims = 2)) == 4 && iszero(prod(A; dims = 2))
+    @test nnz(mapreduce(x -> x + 1, +, A; dims = 2)) == 4
+    @test Matrix(sum(A; dims = 2)) == sum(Matrix(A); dims = 2)
+    # stored and negative zeros survive as they do densely
+    A = sparse([1, 3], [2, 2], [0.0, -0.0], 4, 3)
+    @test isequal(Matrix(sum(A; dims = 1)), sum(Matrix(A); dims = 1))
+    @test isequal(Matrix(sum(A; dims = 2)), sum(Matrix(A); dims = 2))
+    # reductions over empty dimensions
+    for (m, n) in ((0, 4), (4, 0), (0, 0)), dims in (1, 2)
+        A = spzeros(m, n)
+        @test sum(A; dims) isa SparseMatrixCSC
+        @test size(sum(A; dims)) == size(sum(Matrix(A); dims))
+        @test Matrix(prod(A; dims)) == prod(Matrix(A); dims)
+    end
+    # results without a zero, such as the tuples of `extrema`, stay dense
+    A = sprand(5, 4, 0.5)
+    @test extrema(A; dims = 1) isa Matrix
+    @test extrema(A; dims = 1) == extrema(Matrix(A); dims = 1)
+    # the index type is kept
+    A = SparseMatrixCSC{Float64,Int32}(sprand(7, 4, 0.4))
+    @test sum(A; dims = 2) isa SparseMatrixCSC{Float64,Int32}
+    @test sum(A; dims = 1) isa SparseMatrixCSC{Float64,Int32}
+    # hypersparse: only the rows that store something are visited, so reducing along the
+    # rows of a tall matrix costs no more than a few hundred bytes beyond the result
+    A = sparse([5, 10^6, 5], [1, 2, 3], [1.0, 2.0, 3.0], 10^6, 3)
+    r = sum(A; dims = 2)
+    @test nnz(r) == 2 && r[5] == 4.0 && r[10^6] == 2.0
+    @test nnz(sum(A; dims = 1)) == 3
+    sum(A; dims = 2)
+    @test (@allocated sum(A; dims = 2)) < 2^12
+    # destinations: `sum!` resets its destination first, as for dense, while
+    # `mapreducedim!` folds into whatever it already stores
+    A = sprand(8, 6, 0.4); M = Matrix(A)
+    @test sum!(zeros(8, 1), A) ≈ sum(M; dims = 2)
+    @test sum!(spzeros(8, 1), A) ≈ sum(M; dims = 2)
+    @test sum!(spzeros(1, 6), A) ≈ sum(M; dims = 1)
+    R = sparse([2], [1], [1.0], 8, 1)   # partially stored
+    @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 2) .+ [0; 1; 0; 0; 0; 0; 0; 0]
+    R = sparse(ones(8, 1))              # fully stored
+    @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 2) .+ 1
+    R = sparse(ones(1, 6))
+    @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 1) .+ 1
+end
+
 @testset "oneunit of sparse matrix" begin
     A = sparse([Second(0) Second(0); Second(0) Second(0)])
     @test oneunit(sprand(2, 2, 0.5)) isa SparseMatrixCSC{Float64}

From 86a9f55e1ecdceba023257586c9f682430cf55c0 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <virals@gmail.com>
Date: Thu, 10 Sep 2026 09:44:53 +0000
Subject: [PATCH 2/7] Reduce views of a column range through the sparse kernels

Fixes #377.

`sum(view(A, :, j:k); dims)` and the other reductions of a column-range
view went through Base's element-wise fallback, indexing the parent once
per element. The reduction kernels only need the column pointers, row
indices and values, which a `SparseMatrixCSCView` already exposes off
the parent's storage, so they now accept `SparseMatrixCSCUnion`. The two
fully-stored fast paths index through the column pointers rather than
assuming the values start at one, and the hypersparse column reduction
sorts the view's stored range. `nnz` of such a view is now O(1).

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DzbyurS9gH2pdqWJrWaPxD
---
 src/sparsematrix.jl      | 69 ++++++++++++++++++++++------------------
 test/sparsematrix_ops.jl | 10 ++++++
 2 files changed, 48 insertions(+), 31 deletions(-)

diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index e8ef8d2e..e707ae12 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -254,6 +254,11 @@ getnzval( S::SparseMatrixCSCColumnSubset) = nonzeros(parent(S))
 getnzval( S::UpperTriangular{<:Any,<:SparseMatrixCSCOrView}) = nonzeros(S.data)
 getnzval( S::LowerTriangular{<:Any,<:SparseMatrixCSCOrView}) = nonzeros(S.data)
 nzvalview(S::AbstractSparseMatrixCSC) = view(nonzeros(S), 1:nnz(S))
+nzvalview(S::SparseMatrixCSCView) = view(nonzeros(S), _storedrange(S))
+# the stored entries of a column-range view are a contiguous range of the parent's
+_storedrange(S::AbstractSparseMatrixCSC) = 1:nnz(S)
+_storedrange(S::SparseMatrixCSCView) = (colptr = getcolptr(S); Int(colptr[1]):Int(colptr[end]) - 1)
+widelength(S::SparseMatrixCSCView) = prod(Int64.(size(S)))
 
 """
     nnz(A)
@@ -278,9 +283,10 @@ nnz(S::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = nnz(parent(S))
 nnz(S::UpperTriangular{<:Any,<:SparseMatrixCSCOrView}) = nnz1(S)
 nnz(S::LowerTriangular{<:Any,<:SparseMatrixCSCOrView}) = nnz1(S)
 nnz(S::SparseMatrixCSCColumnSubset) = nnz1(S)
+nnz(S::SparseMatrixCSCView) = length(_storedrange(S))
 nnz1(S) = @inbounds sum(length.(nzrange.(Ref(S), axes(S, 2))))
 
-function Base._simple_count(pred, S::AbstractSparseMatrixCSC, init::T) where T
+function Base._simple_count(pred, S::SparseMatrixCSCUnion, init::T) where T
     init + T(count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S)))
 end
 
@@ -2561,8 +2567,9 @@ Base.isequal(A::Transpose{<:Any,<:SparseMatrixCSCMaybeAdjOrTrans}, B::SparseMatr
 # Reductions along a dimension of a sparse matrix return a sparse matrix (issue #43): the
 # result only stores an entry for the rows or columns that store one themselves, unless
 # the reduction of a structurally empty slice is nonzero. The initial array is therefore
-# structurally empty when the initial value is zero, and fully stored otherwise.
-function Base.reducedim_initarray(A::AbstractSparseMatrixCSC{<:Any,Ti}, region, v0, ::Type{R}) where {R,Ti}
+# structurally empty when the initial value is zero, and fully stored otherwise. A view of
+# a column range reduces through the same kernels, off the parent's storage (issue #377).
+function Base.reducedim_initarray(A::SparseMatrixCSCUnion{<:Any,Ti}, region, v0, ::Type{R}) where {R,Ti}
     m, n = Base.to_shape(Base.reduced_indices(A, region))
     if !applicable(zero, R)
         # no structural zero for the result, e.g. the tuples of `extrema`: reduce densely
@@ -2594,7 +2601,7 @@ function _mapreducezeros(f::F, op::G, ::Type{T}, nzeros::Integer, v0) where {F,G
     v
 end
 
-function Base._mapreduce(f::F, op::G, ::Base.IndexCartesian, A::AbstractSparseMatrixCSC{T}) where {F,G,T}
+function Base._mapreduce(f::F, op::G, ::Base.IndexCartesian, A::SparseMatrixCSCUnion{T}) where {F,G,T}
     z = nnz(A)
     n = widelength(A)
     if z == 0
@@ -2619,12 +2626,12 @@ _mapreducezeros(f::Base.ExtremaMap, op::typeof(Base._extrema_rf), ::Type{T}, nze
     nzeros == 0 ? v0 : op(v0, f(zero(T)))
 
 # Specialized mapreduce for any and all
-Base._any(f, A::AbstractSparseMatrixCSC, ::Colon) =
+Base._any(f, A::SparseMatrixCSCUnion, ::Colon) =
     iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)
-Base._all(f, A::AbstractSparseMatrixCSC, ::Colon) =
+Base._all(f, A::SparseMatrixCSCUnion, ::Colon) =
     iszero(widelength(A)) ? true  : Base._mapreduce(f, &, IndexCartesian(), A)
 
-function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::AbstractSparseMatrixCSC{T}) where {F,T}
+function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::SparseMatrixCSCUnion{T}) where {F,T}
     nnzA = nnz(A)
     nzeros = widelength(A) - nnzA
     if nzeros == 0
@@ -2642,14 +2649,14 @@ end
 # filled without touching the slices that store nothing, so the cost stays proportional to
 # the stored entries plus the length of the result; anything in between is rare and goes
 # through the element-wise kernel below.
-function Base._mapreducedim!(f::F, op::G, R::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC{T}) where {F,G,T}
+function Base._mapreducedim!(f::F, op::G, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where {F,G,T}
     require_one_based_indexing(A, R)
     Base.check_reducedims(R, A)
     isempty(A) && return R
     if nnz(R) == length(R)
         Base._mapreducedim!(f, op, reshape(view(nonzeros(R), 1:nnz(R)), size(R)), A)
     elseif nnz(R) != 0
-        invoke(Base._mapreducedim!, Tuple{F,G,AbstractArray,AbstractSparseMatrixCSC{T}}, f, op, R, A)
+        invoke(Base._mapreducedim!, Tuple{F,G,AbstractArray,SparseMatrixCSCUnion{T}}, f, op, R, A)
     elseif size(R) == (1, 1)
         R[1, 1] = op(zero(eltype(R)), mapreduce(f, op, A))
     elseif size(R, 1) == 1
@@ -2664,7 +2671,7 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractSparseMatrixCSC, A::Abstrac
 end
 
 # `R` is a structurally empty `1 x n` sparse matrix: its columns are built in order
-function _mapreducerows_sparse!(f, op, R::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC{T}) where T
+function _mapreducerows_sparse!(f, op, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
     colptr = getcolptr(A)
     nzval = nonzeros(A)
     m, n = size(A)
@@ -2701,7 +2708,7 @@ end
 # vector serves as a dense workspace of length `m` that is then compressed in place; a
 # hypersparse `A` instead has its stored entries sorted by row so that only the rows storing
 # something are ever visited.
-function _mapreducecols_sparse!(f, op, R::AbstractSparseMatrixCSC, A::AbstractSparseMatrixCSC{T}) where T
+function _mapreducecols_sparse!(f, op, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
     m, n = size(A)
     z = zero(eltype(R))
     zempty = n == 0 ? z : _mapreducezeros(f, op, T, n, z)
@@ -2730,8 +2737,8 @@ function _mapreducecols_sparse!(f, op, R::AbstractSparseMatrixCSC, A::AbstractSp
             resize!(Rnzval, k)
         end
     else
-        rows = view(rowvals(A), 1:nnz(A))
-        vals = view(nonzeros(A), 1:nnz(A))
+        rows = view(rowvals(A), _storedrange(A))
+        vals = view(nonzeros(A), _storedrange(A))
         perm = sortperm(rows; alg=Base.Sort.DEFAULT_STABLE)   # keeps each row's entries in column order
         s = 1
         @inbounds while s <= length(perm)
@@ -2752,7 +2759,7 @@ function _mapreducecols_sparse!(f, op, R::AbstractSparseMatrixCSC, A::AbstractSp
 end
 
 # General mapreducedim
-function _mapreducerows!(f, op, R::AbstractArray, A::AbstractSparseMatrixCSC{T}) where T
+function _mapreducerows!(f, op, R::AbstractArray, A::SparseMatrixCSCUnion{T}) where T
     require_one_based_indexing(A, R)
     colptr = getcolptr(A)
     rowval = rowvals(A)
@@ -2768,7 +2775,7 @@ function _mapreducerows!(f, op, R::AbstractArray, A::AbstractSparseMatrixCSC{T})
     R
 end
 
-function _mapreducecols!(f, op, R::AbstractArray, A::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
+function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSCUnion{Tv,Ti}) where {Tv,Ti}
     require_one_based_indexing(A, R)
     colptr = getcolptr(A)
     rowval = rowvals(A)
@@ -2788,7 +2795,7 @@ function _mapreducecols!(f, op, R::AbstractArray, A::AbstractSparseMatrixCSC{Tv,
     R
 end
 
-function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::AbstractSparseMatrixCSC{T}) where {F,G,T}
+function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCUnion{T}) where {F,G,T}
     require_one_based_indexing(A, R)
     lsiz = Base.check_reducedims(R,A)
     isempty(A) && return R
@@ -2806,17 +2813,17 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::AbstractSparseMat
         # Reduction along a dimension > 2
         # Compute op(R, f(A))
         m, n = size(A)
+        colptr = getcolptr(A)
+        rowval = rowvals(A)
         nzval = nonzeros(A)
-        if length(nzval) == m*n
+        if nnz(A) == m*n
             # No zeros, so don't compute f(0) since it might throw
-            for col in axes(A,2)
-                @simd for row in axes(A,1)
-                    @inbounds R[row, col] = op(R[row, col], f(nzval[(col-1)*m+row]))
+            @inbounds for col in axes(A,2)
+                @simd for j = colptr[col]:colptr[col+1]-1
+                    R[rowval[j], col] = op(R[rowval[j], col], f(nzval[j]))
                 end
             end
         else
-            colptr = getcolptr(A)
-            rowval = rowvals(A)
             zeroval = f(zero(T))
             @inbounds for col in axes(A,2)
                 lastrow = 0
@@ -2839,20 +2846,20 @@ end
 
 # Specialized mapreducedim for + cols to avoid allocating a
 # temporary array when f(0) == 0
-function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::AbstractSparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
+function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCUnion{Tv,Ti}) where {Tv,Ti}
     require_one_based_indexing(A, R)
+    colptr = getcolptr(A)
+    rowval = rowvals(A)
     nzval = nonzeros(A)
     m, n = size(A)
-    if length(nzval) == m*n
+    if nnz(A) == m*n
         # No zeros, so don't compute f(0) since it might throw
-        for col in axes(A,2)
-            @simd for row in axes(A,1)
-                @inbounds R[row, 1] = op(R[row, 1], f(nzval[(col-1)*m+row]))
+        @inbounds for col in axes(A,2)
+            @simd for j = colptr[col]:colptr[col+1]-1
+                R[rowval[j], 1] = op(R[rowval[j], 1], f(nzval[j]))
             end
         end
     else
-        colptr = getcolptr(A)
-        rowval = rowvals(A)
         zeroval = f(zero(Tv))
         if isequal(zeroval, zero(Tv))
             # Case where f(0) == 0
@@ -2881,7 +2888,7 @@ end
 
 # any(pred, A, dims = 1) => mapreduce(pred, |, A, dims = 1)
 # final argument `post` is to allow post-mapping each columnar mapreduce
-function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::AbstractSparseMatrixCSC{Tv},
+function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::SparseMatrixCSCUnion{Tv},
                          post::F = identity) where {P, F, Tv}
     nzval = nonzeros(A)
     colptr = getcolptr(A)
@@ -2911,7 +2918,7 @@ function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::Abstr
 end
 # all(pred, A, dims = 1) => mapreduce(pred, &, A, dims = 1) == .!mapreduce(!pred, |, A, dims = 1)
 _mapreducerows!(pred::P, ::typeof(&), R::AbstractMatrix{Bool},
-                A::AbstractSparseMatrixCSC) where {P} = _mapreducerows!(!pred, |, R, A, !)
+                A::SparseMatrixCSCUnion) where {P} = _mapreducerows!(!pred, |, R, A, !)
 
 # findmax/min and argmax/min methods
 # find first zero value in sparse matrix - return linear index in full matrix
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index d844b5f0..9b8cec3f 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -631,12 +631,18 @@ end
                                                  d in (0.0, 0.2, 1.0)
         A = sprand(m, n, d)
         M = Matrix(A)
+        cols = (n + 1) ÷ 2:n   # a view of a column range reduces like its copy (#377)
+        V = view(A, :, cols)
+        C = A[:, cols]
         for dims in (1, 2, (1, 2), 3), f in reductions
             rs = f(A; dims)
             rd = f(M; dims)
             @test rs isa SparseMatrixCSC
             @test size(rs) == size(rd)
             @test Matrix(rs) ≈ rd
+            rv = f(V; dims)
+            rc = f(C; dims)
+            @test typeof(rv) == typeof(rc) && nnz(rv) == nnz(rc) && isequal(rv, rc)
         end
     end
     # only rows and columns that store something get an entry, unless the reduction of a
@@ -674,6 +680,10 @@ end
     @test nnz(sum(A; dims = 1)) == 3
     sum(A; dims = 2)
     @test (@allocated sum(A; dims = 2)) < 2^12
+    # a column-range view goes through the sparse kernels, not the element-wise fallback (#377)
+    V = view(A, :, 2:3)
+    @test (@which Base._mapreducedim!(identity, +, spzeros(10^6, 1), V)).module == SparseArrays
+    @test (@which Base._mapreduce(identity, +, IndexCartesian(), V)).module == SparseArrays
     # destinations: `sum!` resets its destination first, as for dense, while
     # `mapreducedim!` folds into whatever it already stores
     A = sprand(8, 6, 0.4); M = Matrix(A)

From 466dacdb8476a9dc958db051718fc9ca0331c7f0 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <virals@gmail.com>
Date: Fri, 11 Sep 2026 16:11:45 +0000
Subject: [PATCH 3/7] Make the sparse result of reductions along a dimension
 opt-in

Reductions along a dimension of a sparse matrix return a dense `Matrix` again, as
before this PR and as for dense input: the result has one dimension fewer and is
usually dense, and downstream code expects it dense. The sparse result is now
opt-in by reducing into a sparse destination, `sum!(spzeros(size(A, 1), 1), A)` or
`Base.mapreducedim!` and the other in-place reductions, which keeps the hypersparse
kernels and their cost proportional to the stored entries plus the length of the
result. A destination that stores only zeros, as a reused `sum!` destination does
after its `fill!`, folds like an empty one so that reuse stays on the fast path.
Column-range views keep reducing off the parent's storage and now return the same
dense result as their copy, where Base's `similar` gave them a sparse one.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DGaTC2P39YrE2aAU5es55x
---
 docs/src/index.md        |  6 +++
 src/sparsematrix.jl      | 34 ++++++---------
 test/higherorderfns.jl   |  5 +--
 test/sparsematrix_ops.jl | 92 +++++++++++++++-------------------------
 4 files changed, 55 insertions(+), 82 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 132eeb95..3388b629 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -226,6 +226,12 @@ In many cases it may be better to convert the sparse matrix into `(I,J,V)` forma
 manipulate the values or the structure in the dense vectors `(I,J,V)`, and then reconstruct
 the sparse matrix.
 
+Reductions along a dimension, such as `sum(S; dims = 2)`, return a dense `Matrix`, as for dense
+input. To keep the result sparse, reduce into a sparse destination: `sum!(spzeros(size(S, 1), 1), S)`
+stores an entry only for the rows of `S` that store one, at a cost proportional to the number of
+stored entries rather than to the number of rows. `Base.mapreducedim!` and the other in-place
+reductions accept a sparse destination in the same way.
+
 ### [Broadcasting and `map`](@id man-sparse-broadcast)
 
 [`broadcast`](@ref) (including dot syntax such as `A .* B`) and [`map`](@ref) over sparse vectors
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index e707ae12..c3f5c1c7 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -2564,21 +2564,12 @@ Base.isequal(A::Transpose{<:Any,<:SparseMatrixCSCMaybeAdjOrTrans}, B::SparseMatr
 
 ## Reductions
 
-# Reductions along a dimension of a sparse matrix return a sparse matrix (issue #43): the
-# result only stores an entry for the rows or columns that store one themselves, unless
-# the reduction of a structurally empty slice is nonzero. The initial array is therefore
-# structurally empty when the initial value is zero, and fully stored otherwise. A view of
-# a column range reduces through the same kernels, off the parent's storage (issue #377).
-function Base.reducedim_initarray(A::SparseMatrixCSCUnion{<:Any,Ti}, region, v0, ::Type{R}) where {R,Ti}
-    m, n = Base.to_shape(Base.reduced_indices(A, region))
-    if !applicable(zero, R)
-        # no structural zero for the result, e.g. the tuples of `extrema`: reduce densely
-        return fill!(Array{R}(undef, m, n), v0)
-    elseif isequal(v0, zero(R))
-        return spzeros(R, Ti, m, n)
-    else
-        return SparseMatrixCSC(m, n, Ti[1 + m*j for j in 0:n], repeat(Ti.(1:m), n), fill!(Vector{R}(undef, m*n), v0))
-    end
+# Reductions along a dimension return a dense `Array`, as for dense input. A sparse result
+# (issue #43) is opt-in by reducing into a sparse destination, e.g. `sum!(spzeros(size(A, 1), 1), A)`,
+# see `_mapreducedim!` below. Covers column-range views so they reduce like their copy (#377),
+# where Base's `similar` would otherwise give a sparse result.
+function Base.reducedim_initarray(A::SparseMatrixCSCUnion, region, v0, ::Type{R}) where {R}
+    fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
 end
 
 # General mapreduce
@@ -2644,18 +2635,19 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas
     end
 end
 
-# Reduction of a sparse matrix into a sparse destination. A fully stored destination is
-# reduced into as the dense array its stored values form; a structurally empty one is
-# filled without touching the slices that store nothing, so the cost stays proportional to
-# the stored entries plus the length of the result; anything in between is rare and goes
-# through the element-wise kernel below.
+# Reduction into a sparse destination, the opt-in for a sparse result. A fully stored `R` is
+# reduced into as the dense array its values form; an empty one, e.g. `spzeros(m, 1)`, gets an
+# entry only for the rows or columns of `A` that store one (all of them if an empty slice
+# reduces to something nonzero, as for `f(0) != 0`) in time proportional to nnz(A) + length(R);
+# a partially stored `R` is rare and goes through the element-wise kernel below.
 function Base._mapreducedim!(f::F, op::G, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where {F,G,T}
     require_one_based_indexing(A, R)
     Base.check_reducedims(R, A)
     isempty(A) && return R
     if nnz(R) == length(R)
         Base._mapreducedim!(f, op, reshape(view(nonzeros(R), 1:nnz(R)), size(R)), A)
-    elseif nnz(R) != 0
+    elseif nnz(R) != 0 && !all(isequal(zero(eltype(R))), nzvalview(R))
+        # stored zeros only, e.g. a reused `sum!` destination after its `fill!`, fold like an empty one
         invoke(Base._mapreducedim!, Tuple{F,G,AbstractArray,SparseMatrixCSCUnion{T}}, f, op, R, A)
     elseif size(R) == (1, 1)
         R[1, 1] = op(zero(eltype(R)), mapreduce(f, op, A))
diff --git a/test/higherorderfns.jl b/test/higherorderfns.jl
index a70ca758..66483d72 100644
--- a/test/higherorderfns.jl
+++ b/test/higherorderfns.jl
@@ -702,10 +702,7 @@ end
 end
 
 @testset "Issue #27836" begin
-    # the result keeps the reduced eltype (and is sparse, see #43)
-    r = minimum(sparse([1, 2], [1, 2], ones(Int32, 2)), dims = 1)
-    @test r isa SparseMatrixCSC{Int32}
-    @test r == [0 0]
+    @test minimum(sparse([1, 2], [1, 2], ones(Int32, 2)), dims = 1) isa Matrix
 end
 
 @testset "Issue #30118" begin
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index 9b8cec3f..05e25f25 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -614,82 +614,60 @@ end
     @test B ≈ mapreduce(identity, +, Matrix(A), dims=2)
 end
 
-@testset "reductions along a dimension are sparse (#43)" begin
-    # every reduction with a `zero` for its result returns a sparse matrix of the reduced
-    # shape that agrees with the dense result
-    reductions = (
-        (X; dims) -> sum(X; dims), (X; dims) -> prod(X; dims),
-        (X; dims) -> maximum(X; dims), (X; dims) -> minimum(X; dims),
-        (X; dims) -> sum(abs2, X; dims), (X; dims) -> maximum(abs, X; dims),
-        (X; dims) -> count(x -> x > 0.5, X; dims),
-        (X; dims) -> any(x -> x > 0.5, X; dims), (X; dims) -> all(x -> x >= 0, X; dims),
-        (X; dims) -> mapreduce(x -> x + 1, +, X; dims),   # f(0) != 0: dense result
-        (X; dims) -> prod(x -> x + 1, X; dims),
-        (X; dims) -> sum(X; dims, init = 2.5),
+@testset "reductions along a dimension: dense by default, sparse into a sparse destination (#43), column views (#377)" begin
+    reductions = (   # (f, op); the last two have f(0) != 0
+        (identity, +), (identity, *), (identity, max), (identity, min), (abs2, +),
+        (x -> x > 0.5, +), (x -> x > 0.5, |), (x -> x >= 0, &), (x -> x + 1, +), (x -> x + 1, *),
     )
     @testset "size = ($m, $n), density = $d" for (m, n) in ((6, 5), (1, 1), (1, 9), (9, 1), (30, 20)),
                                                  d in (0.0, 0.2, 1.0)
         A = sprand(m, n, d)
         M = Matrix(A)
-        cols = (n + 1) ÷ 2:n   # a view of a column range reduces like its copy (#377)
-        V = view(A, :, cols)
-        C = A[:, cols]
-        for dims in (1, 2, (1, 2), 3), f in reductions
-            rs = f(A; dims)
-            rd = f(M; dims)
-            @test rs isa SparseMatrixCSC
-            @test size(rs) == size(rd)
-            @test Matrix(rs) ≈ rd
-            rv = f(V; dims)
-            rc = f(C; dims)
-            @test typeof(rv) == typeof(rc) && nnz(rv) == nnz(rc) && isequal(rv, rc)
+        V = view(A, :, (n + 1) ÷ 2:n)   # a view of a column range reduces like its copy (#377)
+        C = A[:, (n + 1) ÷ 2:n]
+        @test nnz(V) == nnz(C)
+        for dims in (1, 2, (1, 2), 3), (f, op) in reductions
+            rd = mapreduce(f, op, M; dims)
+            r = mapreduce(f, op, A; dims)
+            @test r isa Matrix && r ≈ rd
+            rv, rc = mapreduce(f, op, V; dims), mapreduce(f, op, C; dims)
+            @test typeof(rv) == typeof(rc) && isequal(rv, rc)
+            # opt-in: reduce into a sparse destination
+            T = eltype(rd)
+            rs = Base.mapreducedim!(f, op, spzeros(T, size(rd)...), A)
+            @test rs isa SparseMatrixCSC{T} && rs ≈ Base.mapreducedim!(f, op, zeros(T, size(rd)...), M)
         end
     end
-    # only rows and columns that store something get an entry, unless the reduction of a
-    # structurally empty slice is nonzero
+    # only rows and columns that store something get an entry, unless a structurally empty
+    # slice reduces to something nonzero
     A = sparse([1, 2], [1, 1], [-1.0, 1.0], 4, 3)
-    @test nnz(sum(A; dims = 1)) == 1   # a stored, cancelled zero
-    @test nnz(sum(A; dims = 2)) == 2
-    @test nnz(prod(A; dims = 2)) == 4 && iszero(prod(A; dims = 2))
-    @test nnz(mapreduce(x -> x + 1, +, A; dims = 2)) == 4
-    @test Matrix(sum(A; dims = 2)) == sum(Matrix(A); dims = 2)
-    # stored and negative zeros survive as they do densely
-    A = sparse([1, 3], [2, 2], [0.0, -0.0], 4, 3)
-    @test isequal(Matrix(sum(A; dims = 1)), sum(Matrix(A); dims = 1))
-    @test isequal(Matrix(sum(A; dims = 2)), sum(Matrix(A); dims = 2))
-    # reductions over empty dimensions
+    @test nnz(sum!(spzeros(1, 3), A)) == 1   # a stored, cancelled zero
+    @test nnz(sum!(spzeros(4, 1), A)) == 2
+    @test nnz(Base.mapreducedim!(x -> x + 1, +, spzeros(4, 1), A)) == 4
+    @test sum!(spzeros(4, 1), A) == sum(Matrix(A); dims = 2)
     for (m, n) in ((0, 4), (4, 0), (0, 0)), dims in (1, 2)
         A = spzeros(m, n)
-        @test sum(A; dims) isa SparseMatrixCSC
-        @test size(sum(A; dims)) == size(sum(Matrix(A); dims))
-        @test Matrix(prod(A; dims)) == prod(Matrix(A); dims)
+        @test sum(A; dims) == sum(Matrix(A); dims)
+        @test sum!(spzeros(size(sum(A; dims))...), A) == sum(Matrix(A); dims)
     end
-    # results without a zero, such as the tuples of `extrema`, stay dense
-    A = sprand(5, 4, 0.5)
-    @test extrema(A; dims = 1) isa Matrix
-    @test extrema(A; dims = 1) == extrema(Matrix(A); dims = 1)
-    # the index type is kept
-    A = SparseMatrixCSC{Float64,Int32}(sprand(7, 4, 0.4))
-    @test sum(A; dims = 2) isa SparseMatrixCSC{Float64,Int32}
-    @test sum(A; dims = 1) isa SparseMatrixCSC{Float64,Int32}
-    # hypersparse: only the rows that store something are visited, so reducing along the
-    # rows of a tall matrix costs no more than a few hundred bytes beyond the result
+    # hypersparse: only the rows that store something are visited
     A = sparse([5, 10^6, 5], [1, 2, 3], [1.0, 2.0, 3.0], 10^6, 3)
-    r = sum(A; dims = 2)
+    r = sum!(spzeros(10^6, 1), A)
     @test nnz(r) == 2 && r[5] == 4.0 && r[10^6] == 2.0
-    @test nnz(sum(A; dims = 1)) == 3
-    sum(A; dims = 2)
-    @test (@allocated sum(A; dims = 2)) < 2^12
+    @test nnz(sum!(spzeros(1, 3), A)) == 3
+    R = spzeros(10^6, 1)
+    sum!(R, A)
+    @test (@allocated sum!(R, A)) < 2^12
     # a column-range view goes through the sparse kernels, not the element-wise fallback (#377)
     V = view(A, :, 2:3)
-    @test (@which Base._mapreducedim!(identity, +, spzeros(10^6, 1), V)).module == SparseArrays
+    @test (@which Base._mapreducedim!(identity, +, zeros(10^6, 1), V)).module == SparseArrays
     @test (@which Base._mapreduce(identity, +, IndexCartesian(), V)).module == SparseArrays
-    # destinations: `sum!` resets its destination first, as for dense, while
-    # `mapreducedim!` folds into whatever it already stores
+    # `sum!` resets its destination first, as for dense; `mapreducedim!` folds into whatever
+    # it already stores
     A = sprand(8, 6, 0.4); M = Matrix(A)
-    @test sum!(zeros(8, 1), A) ≈ sum(M; dims = 2)
     @test sum!(spzeros(8, 1), A) ≈ sum(M; dims = 2)
     @test sum!(spzeros(1, 6), A) ≈ sum(M; dims = 1)
+    @test maximum!(spzeros(8, 1), A) == maximum(M; dims = 2)
     R = sparse([2], [1], [1.0], 8, 1)   # partially stored
     @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 2) .+ [0; 1; 0; 0; 0; 0; 0; 0]
     R = sparse(ones(8, 1))              # fully stored

From e82a89d2c3949db413b57cb6674f42b28f506cd0 Mon Sep 17 00:00:00 2001
From: Viral Shah <viral@mayin.org>
Date: Wed, 16 Sep 2026 17:07:03 -0400
Subject: [PATCH 4/7] Opt in to a sparse reduction result with a keyword
 instead of a sparse destination

`sum(A; dims = 2, sparse = true)` and the other reductions along a dimension
now return a `SparseMatrixCSC`; reducing into a sparse destination is no
longer special-cased. Base forwards unknown keywords from `sum`, `prod`,
`maximum`, `minimum` and `extrema` to `mapreduce`, so one `mapreduce` method
on the sparse types carries the keyword; `any`, `all` and `count` do not
forward it and get their own methods.

Without a destination the kernels can no longer rely on `sum!` and friends
having initialized it, so each slice is now seeded the way Base seeds the
dense result: with `init` when given, otherwise with `mapreduce_first` of its
first stored entry, and the unstored entries folded in afterwards. The
element type and the value of a slice with nothing to reduce come from Base's
`reducedim_init` on a stand-in, so they match the dense result exactly,
including the widening of small integers and the error for `maximum` over an
empty axis.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SgWK99c6dYwt3gxH44MxCs
---
 docs/src/index.md        |   8 +-
 src/sparsematrix.jl      | 174 ++++++++++++++++++++++++---------------
 test/sparsematrix_ops.jl |  93 ++++++++++++++-------
 3 files changed, 174 insertions(+), 101 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 3388b629..f18309cc 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -227,10 +227,10 @@ manipulate the values or the structure in the dense vectors `(I,J,V)`, and then
 the sparse matrix.
 
 Reductions along a dimension, such as `sum(S; dims = 2)`, return a dense `Matrix`, as for dense
-input. To keep the result sparse, reduce into a sparse destination: `sum!(spzeros(size(S, 1), 1), S)`
-stores an entry only for the rows of `S` that store one, at a cost proportional to the number of
-stored entries rather than to the number of rows. `Base.mapreducedim!` and the other in-place
-reductions accept a sparse destination in the same way.
+input. To keep the result sparse, pass `sparse = true`: `sum(S; dims = 2, sparse = true)` stores an
+entry only for the rows of `S` that store one, at a cost proportional to the number of stored
+entries rather than to the number of rows. `prod`, `maximum`, `minimum`, `count`, `any`, `all` and
+`mapreduce` accept the keyword in the same way.
 
 ### [Broadcasting and `map`](@id man-sparse-broadcast)
 
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index c3f5c1c7..303ae4c7 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -2565,9 +2565,9 @@ Base.isequal(A::Transpose{<:Any,<:SparseMatrixCSCMaybeAdjOrTrans}, B::SparseMatr
 ## Reductions
 
 # Reductions along a dimension return a dense `Array`, as for dense input. A sparse result
-# (issue #43) is opt-in by reducing into a sparse destination, e.g. `sum!(spzeros(size(A, 1), 1), A)`,
-# see `_mapreducedim!` below. Covers column-range views so they reduce like their copy (#377),
-# where Base's `similar` would otherwise give a sparse result.
+# (issue #43) is opt-in with the keyword `sparse = true`, see `_mapreduce_dim_sparse` below.
+# Covers column-range views so they reduce like their copy (#377), where Base's `similar`
+# would otherwise give a sparse result.
 function Base.reducedim_initarray(A::SparseMatrixCSCUnion, region, v0, ::Type{R}) where {R}
     fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
 end
@@ -2635,44 +2635,81 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas
     end
 end
 
-# Reduction into a sparse destination, the opt-in for a sparse result. A fully stored `R` is
-# reduced into as the dense array its values form; an empty one, e.g. `spzeros(m, 1)`, gets an
-# entry only for the rows or columns of `A` that store one (all of them if an empty slice
-# reduces to something nonzero, as for `f(0) != 0`) in time proportional to nnz(A) + length(R);
-# a partially stored `R` is rare and goes through the element-wise kernel below.
-function Base._mapreducedim!(f::F, op::G, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where {F,G,T}
-    require_one_based_indexing(A, R)
-    Base.check_reducedims(R, A)
-    isempty(A) && return R
-    if nnz(R) == length(R)
-        Base._mapreducedim!(f, op, reshape(view(nonzeros(R), 1:nnz(R)), size(R)), A)
-    elseif nnz(R) != 0 && !all(isequal(zero(eltype(R))), nzvalview(R))
-        # stored zeros only, e.g. a reused `sum!` destination after its `fill!`, fold like an empty one
-        invoke(Base._mapreducedim!, Tuple{F,G,AbstractArray,SparseMatrixCSCUnion{T}}, f, op, R, A)
-    elseif size(R) == (1, 1)
-        R[1, 1] = op(zero(eltype(R)), mapreduce(f, op, A))
-    elseif size(R, 1) == 1
-        _mapreducerows_sparse!(f, op, R, A)
-    elseif size(R, 2) == 1
-        _mapreducecols_sparse!(f, op, R, A)
+# The `sparse` keyword, the opt-in for a sparse result (issue #43): `sum(A; dims = 2, sparse = true)`
+# and the other reductions along a dimension return a `SparseMatrixCSC` that stores an entry
+# only for the rows or columns of `A` that store one (all of them when a slice that stores
+# nothing reduces to something nonzero, as for `f(0) != 0`), in time proportional to
+# nnz(A) + length(result) rather than to length(A). Base's `sum`, `prod`, `maximum`, `minimum`
+# and `extrema` forward unknown keywords to `mapreduce`, so this method serves them all;
+# `any`, `all` and `count` do not and get their own methods below.
+function Base.mapreduce(f, op, A::SparseMatrixCSCUnion; dims=:, init=Base._InitialValue(), sparse::Bool=false)
+    sparse || return Base._mapreduce_dim(f, op, init, A, dims)
+    dims === (:) && throw(ArgumentError("a sparse result needs a reduction along a dimension, pass `dims`"))
+    return _mapreduce_dim_sparse(f, op, init, A, dims)
+end
+for (fname, _fname, op) in ((:any, :_any, :(Base.or_any)), (:all, :_all, :(Base.and_all)))
+    @eval begin
+        Base.$fname(A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) = Base.$fname(identity, A; dims, sparse)
+        Base.$fname(f::Function, A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) =
+            sparse ? mapreduce(f, $op, A; dims, sparse) : Base.$_fname(f, A, dims)
+    end
+end
+Base.count(A::SparseMatrixCSCUnion; dims=:, init=0, sparse::Bool=false) = count(identity, A; dims, init, sparse)
+Base.count(f, A::SparseMatrixCSCUnion; dims=:, init=0, sparse::Bool=false) =
+    sparse ? mapreduce(Base._bool(f), Base.add_sum, A; dims, init, sparse) : Base._count(f, A, dims, init)
+
+# The slices are reduced the way Base's dense result is: seeded with `init` when given and
+# otherwise with `mapreduce_first` of their first element, with the entries a slice does not
+# store folded in through `_mapreducezeros`.
+_seed(f, op, ::Base._InitialValue, x) = Base.mapreduce_first(f, op, x)
+_seed(f, op, init, x) = op(init, f(x))
+# element type of the dense result, taken from Base's initialization of a 1 x 1 stand-in
+_reduced_eltype(f, op, ::Base._InitialValue, ::Type{T}) where T =
+    eltype(Base.reducedim_init(f, op, fill(zero(T), 1, 1), 1))
+_reduced_eltype(f, op, init, ::Type{T}) where T = typeof(init)
+# the reduction of a slice with no entries at all, as Base initializes the dense result
+_reduced_empty(f, op, ::Base._InitialValue, ::Type{T}) where T =
+    Base.reducedim_init(f, op, Matrix{T}(undef, 0, 1), 1)[1]
+_reduced_empty(f, op, init, ::Type{T}) where T = init
+# the reduction of a slice of length `len` that stores nothing
+_reduce_unstored(f, op, init, ::Type{T}, len) where T =
+    len == 0 ? _reduced_empty(f, op, init, T) : _mapreducezeros(f, op, T, len - 1, _seed(f, op, init, zero(T)))
+
+function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCUnion{T,Ti}, dims) where {T,Ti}
+    m, n = size(A)
+    rm, rn = map(length, Base.reduced_indices(A, dims))   # also validates `dims`
+    Tr = _reduced_eltype(f, op, init, T)
+    if rm == rn == 1
+        R = spzeros(Tr, Ti, 1, 1)
+        v = isempty(A) ? _reduced_empty(f, op, init, T) : _seed(identity, op, init, mapreduce(f, op, A))
+        if nnz(A) > 0 || !isequal(v, zero(Tr))
+            push!(rowvals(R), 1)
+            push!(nonzeros(R), v)
+            getcolptr(R)[2] = 2
+        end
+        return R
+    elseif rm == 1
+        return _mapreducerows_sparse!(f, op, init, spzeros(Tr, Ti, 1, n), A)
+    elseif rn == 1
+        return _mapreducecols_sparse!(f, op, init, spzeros(Tr, Ti, m, 1), A)
     else
-        # reduction over a dimension beyond 2: `R` has the shape of `A`
-        copyto!(R, op.(zero(eltype(R)), f.(A)))
+        # a dimension beyond 2: every entry is a slice of its own
+        return convert(SparseMatrixCSC{Tr,Ti}, map(x -> _seed(f, op, init, x), A))
     end
-    return R
 end
 
 # `R` is a structurally empty `1 x n` sparse matrix: its columns are built in order
-function _mapreducerows_sparse!(f, op, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
+function _mapreducerows_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
     colptr = getcolptr(A)
     nzval = nonzeros(A)
     m, n = size(A)
     z = zero(eltype(R))
-    # the reduction of a column that stores nothing, stored only when it is nonzero
-    zempty = m == 0 ? z : _mapreducezeros(f, op, T, m, z)
-    store_empty = !isequal(zempty, z)
+    # the reduction of a column that stores nothing; when every column is full it is not
+    # needed and f(0) is not evaluated, since it might throw
+    zunstored = nnz(A) == m*n && m > 0 ? z : _reduce_unstored(f, op, init, T, m)
+    store_unstored = !isequal(zunstored, z)
     Rcolptr, Rrowval, Rnzval = getcolptr(R), rowvals(R), nonzeros(R)
-    nstored = store_empty ? n : count(col -> colptr[col+1] > colptr[col], 1:n)
+    nstored = store_unstored ? n : count(col -> colptr[col+1] > colptr[col], 1:n)
     resize!(Rrowval, nstored)
     fill!(Rrowval, 1)
     resize!(Rnzval, nstored)
@@ -2680,11 +2717,11 @@ function _mapreducerows_sparse!(f, op, R::AbstractSparseMatrixCSC, A::SparseMatr
     @inbounds for col in 1:n
         rng = colptr[col]:colptr[col+1]-1
         if isempty(rng)
-            store_empty || (Rcolptr[col+1] = k + 1; continue)
-            v = zempty
+            store_unstored || (Rcolptr[col+1] = k + 1; continue)
+            v = zunstored
         else
-            r = z
-            @simd for j in rng
+            r = _seed(f, op, init, nzval[first(rng)])
+            @simd for j in first(rng)+1:last(rng)
                 r = op(r, f(nzval[j]))
             end
             v = _mapreducezeros(f, op, T, m - length(rng), r)
@@ -2696,48 +2733,51 @@ function _mapreducerows_sparse!(f, op, R::AbstractSparseMatrixCSC, A::SparseMatr
     return R
 end
 
-# `R` is a structurally empty `m x 1` sparse matrix. With enough stored entries its value
-# vector serves as a dense workspace of length `m` that is then compressed in place; a
-# hypersparse `A` instead has its stored entries sorted by row so that only the rows storing
-# something are ever visited.
-function _mapreducecols_sparse!(f, op, R::AbstractSparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
+# `R` is a structurally empty `m x 1` sparse matrix. With enough stored entries the rows are
+# reduced in a dense workspace that is then compressed into `R`; a hypersparse `A` instead has
+# its stored entries sorted by row so that only the rows storing something are ever visited.
+function _mapreducecols_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
     m, n = size(A)
-    z = zero(eltype(R))
-    zempty = n == 0 ? z : _mapreducezeros(f, op, T, n, z)
-    store_empty = !isequal(zempty, z)
+    Tr = eltype(R)
+    z = zero(Tr)
+    rows = view(rowvals(A), _storedrange(A))
+    vals = view(nonzeros(A), _storedrange(A))
+    nz = length(rows)
+    zunstored = nz == m*n && n > 0 ? z : _reduce_unstored(f, op, init, T, n)
+    store_unstored = !isequal(zunstored, z)
     Rcolptr, Rrowval, Rnzval = getcolptr(R), rowvals(R), nonzeros(R)
-    resize!(Rrowval, 0)
-    resize!(Rnzval, 0)
-    if store_empty || 8 * nnz(A) >= m
-        resize!(Rnzval, m)
-        fill!(Rnzval, z)
-        _mapreducecols!(f, op, Rnzval, A)
+    if store_unstored || 8 * nz >= m
+        W = Vector{Tr}(undef, m)
+        cnt = zeros(Int, m)   # stored entries seen in each row
+        @inbounds for j in eachindex(rows, vals)
+            i = rows[j]
+            W[i] = cnt[i] == 0 ? _seed(f, op, init, vals[j]) : op(W[i], f(vals[j]))
+            cnt[i] += 1
+        end
         resize!(Rrowval, m)
-        if store_empty
-            Rrowval .= 1:m
-        else
-            k = 0
-            @inbounds for i in 1:m
-                w = Rnzval[i]
-                if !isequal(w, z)
-                    k += 1
-                    Rrowval[k] = i
-                    Rnzval[k] = w
-                end
+        resize!(Rnzval, m)
+        k = 0
+        @inbounds for i in 1:m
+            if cnt[i] > 0
+                k += 1
+                Rrowval[k] = i
+                Rnzval[k] = _mapreducezeros(f, op, T, n - cnt[i], W[i])
+            elseif store_unstored
+                k += 1
+                Rrowval[k] = i
+                Rnzval[k] = zunstored
             end
-            resize!(Rrowval, k)
-            resize!(Rnzval, k)
         end
+        resize!(Rrowval, k)
+        resize!(Rnzval, k)
     else
-        rows = view(rowvals(A), _storedrange(A))
-        vals = view(nonzeros(A), _storedrange(A))
         perm = sortperm(rows; alg=Base.Sort.DEFAULT_STABLE)   # keeps each row's entries in column order
         s = 1
-        @inbounds while s <= length(perm)
+        @inbounds while s <= nz
             row = rows[perm[s]]
-            r = op(z, f(vals[perm[s]]))
+            r = _seed(f, op, init, vals[perm[s]])
             t = s + 1
-            while t <= length(perm) && rows[perm[t]] == row
+            while t <= nz && rows[perm[t]] == row
                 r = op(r, f(vals[perm[t]]))
                 t += 1
             end
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index 05e25f25..ec8655ed 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -614,14 +614,14 @@ end
     @test B ≈ mapreduce(identity, +, Matrix(A), dims=2)
 end
 
-@testset "reductions along a dimension: dense by default, sparse into a sparse destination (#43), column views (#377)" begin
+@testset "reductions along a dimension: dense by default, sparse with `sparse = true` (#43), column views (#377)" begin
     reductions = (   # (f, op); the last two have f(0) != 0
         (identity, +), (identity, *), (identity, max), (identity, min), (abs2, +),
         (x -> x > 0.5, +), (x -> x > 0.5, |), (x -> x >= 0, &), (x -> x + 1, +), (x -> x + 1, *),
     )
     @testset "size = ($m, $n), density = $d" for (m, n) in ((6, 5), (1, 1), (1, 9), (9, 1), (30, 20)),
                                                  d in (0.0, 0.2, 1.0)
-        A = sprand(m, n, d)
+        A = sparse(sprand(m, n, d) .- 0.5)   # negative entries, so that max and min do not see 0 as a bound
         M = Matrix(A)
         V = view(A, :, (n + 1) ÷ 2:n)   # a view of a column range reduces like its copy (#377)
         C = A[:, (n + 1) ÷ 2:n]
@@ -632,48 +632,81 @@ end
             @test r isa Matrix && r ≈ rd
             rv, rc = mapreduce(f, op, V; dims), mapreduce(f, op, C; dims)
             @test typeof(rv) == typeof(rc) && isequal(rv, rc)
-            # opt-in: reduce into a sparse destination
+            # opt-in: the sparse result has the element type and values of the dense one
             T = eltype(rd)
-            rs = Base.mapreducedim!(f, op, spzeros(T, size(rd)...), A)
-            @test rs isa SparseMatrixCSC{T} && rs ≈ Base.mapreducedim!(f, op, zeros(T, size(rd)...), M)
+            rs = mapreduce(f, op, A; dims, sparse = true)
+            @test rs isa SparseMatrixCSC{T} && rs ≈ rd
+            rvs = mapreduce(f, op, V; dims, sparse = true)
+            @test rvs isa SparseMatrixCSC{T} && rvs ≈ mapreduce(f, op, Matrix(C); dims)
+        end
+        for dims in (1, 2)
+            @test sum(A; dims, sparse = true) ≈ sum(M; dims)
+            @test sum(abs, V; dims, sparse = true) ≈ sum(abs, Matrix(C); dims)
+            @test prod(A; dims, sparse = true) ≈ prod(M; dims)
+            @test maximum(A; dims, sparse = true) == maximum(M; dims)
+            @test minimum(abs2, A; dims, sparse = true) == minimum(abs2, M; dims)
+            @test sum(A; dims, init = 2.5, sparse = true) ≈ sum(M; dims, init = 2.5)
+            @test mapreduce(abs, (x, y) -> x + y, A; dims, init = 1.5, sparse = true) ≈
+                  mapreduce(abs, (x, y) -> x + y, M; dims, init = 1.5)
+            @test count(>(0), A; dims, sparse = true) == count(>(0), M; dims)
+            @test count(A .> 0; dims, sparse = true) == count(M .> 0; dims)
+            @test count(A .> 0; dims, init = 3, sparse = true) == count(M .> 0; dims, init = 3)
+            @test any(>(0), A; dims, sparse = true) == any(>(0), M; dims)
+            @test any(A .> 0; dims, sparse = true) == any(M .> 0; dims)
+            @test all(<(0.4), A; dims, sparse = true) == all(<(0.4), M; dims)
+            @test all(A .< 0.4; dims, sparse = true) == all(M .< 0.4; dims)
+            for r in (count(>(0), A; dims, sparse = true), any(A .> 0; dims, sparse = true), all(A .< 0.4; dims, sparse = true))
+                @test r isa SparseMatrixCSC
+            end
+            # the default result and the scalar reductions are unchanged
+            @test sum(A; dims) isa Matrix{Float64} && count(A .> 0; dims) isa Matrix{Int} && any(A .> 0; dims) isa Matrix{Bool}
         end
+        @test sum(A) ≈ sum(M) && count(>(0), A) == count(>(0), M) && any(A .> 0) == any(M .> 0) && all(A .< 0.4) == all(M .< 0.4)
+        @test_throws ArgumentError sum(A; sparse = true)
     end
-    # only rows and columns that store something get an entry, unless a structurally empty
-    # slice reduces to something nonzero
+    # only rows and columns that store something get an entry, unless a slice that stores
+    # nothing reduces to something nonzero
     A = sparse([1, 2], [1, 1], [-1.0, 1.0], 4, 3)
-    @test nnz(sum!(spzeros(1, 3), A)) == 1   # a stored, cancelled zero
-    @test nnz(sum!(spzeros(4, 1), A)) == 2
-    @test nnz(Base.mapreducedim!(x -> x + 1, +, spzeros(4, 1), A)) == 4
-    @test sum!(spzeros(4, 1), A) == sum(Matrix(A); dims = 2)
-    for (m, n) in ((0, 4), (4, 0), (0, 0)), dims in (1, 2)
+    @test nnz(sum(A; dims = 1, sparse = true)) == 1   # a stored, cancelled zero
+    @test nnz(sum(A; dims = 2, sparse = true)) == 2
+    @test nnz(sum(x -> x + 1, A; dims = 2, sparse = true)) == 4
+    @test nnz(sum(A; dims = 2, init = 1.0, sparse = true)) == 4
+    @test nnz(prod(A; dims = 1, sparse = true)) == 1   # the product of an unstored column is 0
+    @test sum(A; dims = 2, sparse = true) == sum(Matrix(A); dims = 2)
+    # the element type is that of the dense result
+    @test sum(sparse(Int8[1 2; 3 4]); dims = 1, sparse = true) isa SparseMatrixCSC{Int}
+    @test sum(sparse([true false]); dims = 2, sparse = true) isa SparseMatrixCSC{Int}
+    @test maximum(sparse(Int8[1 2; 3 4]); dims = 1, sparse = true) isa SparseMatrixCSC{Int8}
+    @test sum(sparse(Int8[1 2; 3 4]); dims = 1, init = Int8(1), sparse = true) isa SparseMatrixCSC{Int8}
+    # empty dimensions
+    for (m, n) in ((0, 4), (4, 0), (0, 0)), dims in (1, 2, (1, 2))
         A = spzeros(m, n)
         @test sum(A; dims) == sum(Matrix(A); dims)
-        @test sum!(spzeros(size(sum(A; dims))...), A) == sum(Matrix(A); dims)
+        @test sum(A; dims, sparse = true) == sum(Matrix(A); dims)
+        @test prod(A; dims, sparse = true) == prod(Matrix(A); dims)
+        @test sum(x -> x + 1, A; dims, sparse = true) == sum(x -> x + 1, Matrix(A); dims)
+        @test all(A .> 0; dims, sparse = true) == all(Matrix(A) .> 0; dims)
+        md = try maximum(Matrix(A); dims) catch err; err end   # throws over an empty axis
+        if md isa ArgumentError
+            @test_throws ArgumentError maximum(A; dims, sparse = true)
+        else
+            @test maximum(A; dims, sparse = true) == md
+        end
     end
+    @test_throws ArgumentError sum(spzeros(3, 3); dims = 0, sparse = true)
     # hypersparse: only the rows that store something are visited
     A = sparse([5, 10^6, 5], [1, 2, 3], [1.0, 2.0, 3.0], 10^6, 3)
-    r = sum!(spzeros(10^6, 1), A)
+    r = sum(A; dims = 2, sparse = true)
     @test nnz(r) == 2 && r[5] == 4.0 && r[10^6] == 2.0
-    @test nnz(sum!(spzeros(1, 3), A)) == 3
-    R = spzeros(10^6, 1)
-    sum!(R, A)
-    @test (@allocated sum!(R, A)) < 2^12
+    @test maximum(A; dims = 2, sparse = true) == maximum(Matrix(A); dims = 2)
+    @test nnz(sum(A; dims = 1, sparse = true)) == 3
+    sum(A; dims = 2, sparse = true)
+    @test (@allocated sum(A; dims = 2, sparse = true)) < 2^12
     # a column-range view goes through the sparse kernels, not the element-wise fallback (#377)
     V = view(A, :, 2:3)
     @test (@which Base._mapreducedim!(identity, +, zeros(10^6, 1), V)).module == SparseArrays
     @test (@which Base._mapreduce(identity, +, IndexCartesian(), V)).module == SparseArrays
-    # `sum!` resets its destination first, as for dense; `mapreducedim!` folds into whatever
-    # it already stores
-    A = sprand(8, 6, 0.4); M = Matrix(A)
-    @test sum!(spzeros(8, 1), A) ≈ sum(M; dims = 2)
-    @test sum!(spzeros(1, 6), A) ≈ sum(M; dims = 1)
-    @test maximum!(spzeros(8, 1), A) == maximum(M; dims = 2)
-    R = sparse([2], [1], [1.0], 8, 1)   # partially stored
-    @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 2) .+ [0; 1; 0; 0; 0; 0; 0; 0]
-    R = sparse(ones(8, 1))              # fully stored
-    @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 2) .+ 1
-    R = sparse(ones(1, 6))
-    @test Base.mapreducedim!(identity, +, R, A) ≈ sum(M; dims = 1) .+ 1
+    @test nnz(sum(V; dims = 2, sparse = true)) == 2
 end
 
 @testset "oneunit of sparse matrix" begin

From abe8e9487d5a69814c67e480518537d4c21174db Mon Sep 17 00:00:00 2001
From: Viral Shah <viral@mayin.org>
Date: Thu, 17 Sep 2026 15:38:56 -0400
Subject: [PATCH 5/7] Reject element types without a zero, accept any callable
 in any/all, trim the tests

`extrema(A; dims, sparse = true)` failed with a MethodError on `zero(Tuple)`;
it now throws an ArgumentError. `any` and `all` with `sparse = true` accept
callables that are not `Function`s. The reduction grid in the tests drops
three pairs that exercise no new path, and gains complex, callable and
`extrema` cases.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/sparsematrix.jl      | 13 ++++++++-----
 test/sparsematrix_ops.jl | 16 +++++++++++++---
 2 files changed, 21 insertions(+), 8 deletions(-)

diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index 303ae4c7..26278078 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -2639,9 +2639,9 @@ end
 # and the other reductions along a dimension return a `SparseMatrixCSC` that stores an entry
 # only for the rows or columns of `A` that store one (all of them when a slice that stores
 # nothing reduces to something nonzero, as for `f(0) != 0`), in time proportional to
-# nnz(A) + length(result) rather than to length(A). Base's `sum`, `prod`, `maximum`, `minimum`
-# and `extrema` forward unknown keywords to `mapreduce`, so this method serves them all;
-# `any`, `all` and `count` do not and get their own methods below.
+# nnz(A) + length(result) rather than to length(A). Base's `sum`, `prod`, `maximum` and
+# `minimum` forward unknown keywords to `mapreduce`, so this method serves them all; `any`,
+# `all` and `count` do not and get their own methods below.
 function Base.mapreduce(f, op, A::SparseMatrixCSCUnion; dims=:, init=Base._InitialValue(), sparse::Bool=false)
     sparse || return Base._mapreduce_dim(f, op, init, A, dims)
     dims === (:) && throw(ArgumentError("a sparse result needs a reduction along a dimension, pass `dims`"))
@@ -2650,7 +2650,7 @@ end
 for (fname, _fname, op) in ((:any, :_any, :(Base.or_any)), (:all, :_all, :(Base.and_all)))
     @eval begin
         Base.$fname(A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) = Base.$fname(identity, A; dims, sparse)
-        Base.$fname(f::Function, A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) =
+        Base.$fname(f, A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) =
             sparse ? mapreduce(f, $op, A; dims, sparse) : Base.$_fname(f, A, dims)
     end
 end
@@ -2663,7 +2663,8 @@ Base.count(f, A::SparseMatrixCSCUnion; dims=:, init=0, sparse::Bool=false) =
 # store folded in through `_mapreducezeros`.
 _seed(f, op, ::Base._InitialValue, x) = Base.mapreduce_first(f, op, x)
 _seed(f, op, init, x) = op(init, f(x))
-# element type of the dense result, taken from Base's initialization of a 1 x 1 stand-in
+# element type of the dense result, taken from Base's initialization of a 1 x 1 stand-in,
+# which evaluates f(0) as Base does
 _reduced_eltype(f, op, ::Base._InitialValue, ::Type{T}) where T =
     eltype(Base.reducedim_init(f, op, fill(zero(T), 1, 1), 1))
 _reduced_eltype(f, op, init, ::Type{T}) where T = typeof(init)
@@ -2679,6 +2680,8 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCUnion{T,Ti}, dims)
     m, n = size(A)
     rm, rn = map(length, Base.reduced_indices(A, dims))   # also validates `dims`
     Tr = _reduced_eltype(f, op, init, T)
+    applicable(zero, Tr) || throw(ArgumentError("cannot store a sparse result of element type $Tr, " *
+        "which has no zero (as for `extrema`); reduce without `sparse = true`"))
     if rm == rn == 1
         R = spzeros(Tr, Ti, 1, 1)
         v = isempty(A) ? _reduced_empty(f, op, init, T) : _seed(identity, op, init, mapreduce(f, op, A))
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index ec8655ed..9269ba1c 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -615,9 +615,8 @@ end
 end
 
 @testset "reductions along a dimension: dense by default, sparse with `sparse = true` (#43), column views (#377)" begin
-    reductions = (   # (f, op); the last two have f(0) != 0
-        (identity, +), (identity, *), (identity, max), (identity, min), (abs2, +),
-        (x -> x > 0.5, +), (x -> x > 0.5, |), (x -> x >= 0, &), (x -> x + 1, +), (x -> x + 1, *),
+    reductions = (   # (f, op); the last one has f(0) != 0
+        (identity, +), (identity, *), (identity, max), (abs2, +), (x -> x > 0.5, |), (x -> x >= 0, &), (x -> x + 1, +),
     )
     @testset "size = ($m, $n), density = $d" for (m, n) in ((6, 5), (1, 1), (1, 9), (9, 1), (30, 20)),
                                                  d in (0.0, 0.2, 1.0)
@@ -664,6 +663,17 @@ end
         @test sum(A) ≈ sum(M) && count(>(0), A) == count(>(0), M) && any(A .> 0) == any(M .> 0) && all(A .< 0.4) == all(M .< 0.4)
         @test_throws ArgumentError sum(A; sparse = true)
     end
+    C = sprand(ComplexF64, 6, 5, 0.3)
+    MC, VC = Matrix(C), view(C, :, 2:5)
+    for dims in (1, 2)
+        @test sum(C; dims, sparse = true) isa SparseMatrixCSC{ComplexF64} && sum(C; dims, sparse = true) ≈ sum(MC; dims)
+        @test prod(abs2, C; dims, sparse = true) ≈ prod(abs2, MC; dims)
+        @test sum(VC; dims) isa Matrix{ComplexF64} && sum(VC; dims) == sum(Matrix(VC); dims)
+    end
+    struct Positive end   # a callable that is not a `Function`
+    (::Positive)(x) = x > 0
+    @test any(Positive(), C .|> real; dims = 1, sparse = true) == any(Positive(), real.(MC); dims = 1)
+    @test_throws ArgumentError extrema(C; dims = 1, sparse = true)   # a tuple has no zero
     # only rows and columns that store something get an entry, unless a slice that stores
     # nothing reduces to something nonzero
     A = sparse([1, 2], [1, 1], [-1.0, 1.0], 4, 3)

From d088871a346cda75ccd2ff1641038e45c91d58cb Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Sat, 19 Sep 2026 10:20:21 -0400
Subject: [PATCH 6/7] Reduce adjoints, column-subset views and sparse vectors
 through the sparse kernels

`sum(A'; dims)` and reductions of `view(A, :, cols)` got a fully stored sparse
destination from Base's `similar` and were filled through sparse `setindex!` or the
element-wise fallback. They now get the dense destination their copy gets. The
reduction kernels walk `nzrange` instead of the column pointers, so they take a view
of any column subset; an adjoint or transpose whose reduction LinearAlgebra does not
forward (a non-commutative `op`) is reduced as its parent along the other dimension.
The `sparse = true` keyword is accepted for all of these and for sparse vectors,
which return a `SparseVector`.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 docs/src/index.md        |   3 +-
 src/abstractsparse.jl    |   4 ++
 src/sparsematrix.jl      | 131 ++++++++++++++++++++++++---------------
 src/sparsevector.jl      |   3 -
 test/sparsematrix_ops.jl |  24 +++++++
 5 files changed, 110 insertions(+), 55 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index f18309cc..8177bfe1 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -230,7 +230,8 @@ Reductions along a dimension, such as `sum(S; dims = 2)`, return a dense `Matrix
 input. To keep the result sparse, pass `sparse = true`: `sum(S; dims = 2, sparse = true)` stores an
 entry only for the rows of `S` that store one, at a cost proportional to the number of stored
 entries rather than to the number of rows. `prod`, `maximum`, `minimum`, `count`, `any`, `all` and
-`mapreduce` accept the keyword in the same way.
+`mapreduce` accept the keyword in the same way, as do adjoints and transposes of sparse matrices,
+views of a subset of their columns, and sparse vectors, for which the result is a `SparseVector`.
 
 ### [Broadcasting and `map`](@id man-sparse-broadcast)
 
diff --git a/src/abstractsparse.jl b/src/abstractsparse.jl
index aaa2c935..b67386f0 100644
--- a/src/abstractsparse.jl
+++ b/src/abstractsparse.jl
@@ -75,6 +75,10 @@ const SparseVectorPartialView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseVector{Tv,T
 const SparseMatrixCSCMaybeAdjOrTrans = Union{AbstractSparseMatrixCSC, AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}
 const SparseVecOrMatMaybeAdjOrTrans = Union{SparseVecOrMat, AdjOrTrans{<:Any,<:SparseVecOrMat}}
 
+# arguments whose reductions along a dimension take the `sparse` keyword
+const SparseReducible = Union{SparseMatrixCSCOrColumnSubset, AdjOrTrans{<:Any,<:SparseMatrixCSCOrColumnSubset},
+                              SparseVectorOrView}
+
 # LinearAlgebra wrappers around CSC storage
 const SparseTriangular{Tv,Ti} = UpperOrLowerTriangular{Tv,<:SparseMatrixCSCOrView{Tv,Ti}}
 const SparseOrTri{Tv,Ti} = Union{SparseMatrixCSCOrView{Tv,Ti}, SparseTriangular{Tv,Ti}}
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index 26278078..6dfa9ed2 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -254,11 +254,18 @@ getnzval( S::SparseMatrixCSCColumnSubset) = nonzeros(parent(S))
 getnzval( S::UpperTriangular{<:Any,<:SparseMatrixCSCOrView}) = nonzeros(S.data)
 getnzval( S::LowerTriangular{<:Any,<:SparseMatrixCSCOrView}) = nonzeros(S.data)
 nzvalview(S::AbstractSparseMatrixCSC) = view(nonzeros(S), 1:nnz(S))
-nzvalview(S::SparseMatrixCSCView) = view(nonzeros(S), _storedrange(S))
-# the stored entries of a column-range view are a contiguous range of the parent's
-_storedrange(S::AbstractSparseMatrixCSC) = 1:nnz(S)
-_storedrange(S::SparseMatrixCSCView) = (colptr = getcolptr(S); Int(colptr[1]):Int(colptr[end]) - 1)
-widelength(S::SparseMatrixCSCView) = prod(Int64.(size(S)))
+nzvalview(S::SparseMatrixCSCColumnSubset) = view(nonzeros(S), _storedinds(S))
+# where the stored entries sit in the parent's storage: a contiguous range for a column range
+_storedinds(S::AbstractSparseMatrixCSC) = 1:nnz(S)
+_storedinds(S::SparseMatrixCSCView) = (colptr = getcolptr(S); Int(colptr[1]):Int(colptr[end]) - 1)
+function _storedinds(S::SparseMatrixCSCColumnSubset)
+    inds = Int[]
+    for col in axes(S, 2)
+        append!(inds, nzrange(S, col))
+    end
+    return inds
+end
+widelength(S::SparseMatrixCSCColumnSubset) = prod(Int64.(size(S)))
 
 """
     nnz(A)
@@ -283,10 +290,10 @@ nnz(S::AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}) = nnz(parent(S))
 nnz(S::UpperTriangular{<:Any,<:SparseMatrixCSCOrView}) = nnz1(S)
 nnz(S::LowerTriangular{<:Any,<:SparseMatrixCSCOrView}) = nnz1(S)
 nnz(S::SparseMatrixCSCColumnSubset) = nnz1(S)
-nnz(S::SparseMatrixCSCView) = length(_storedrange(S))
+nnz(S::SparseMatrixCSCView) = length(_storedinds(S))
 nnz1(S) = @inbounds sum(length.(nzrange.(Ref(S), axes(S, 2))))
 
-function Base._simple_count(pred, S::SparseMatrixCSCUnion, init::T) where T
+function Base._simple_count(pred, S::SparseMatrixCSCOrColumnSubset, init::T) where T
     init + T(count(pred, nzvalview(S)) + pred(zero(eltype(S)))*(prod(size(S)) - nnz(S)))
 end
 
@@ -2566,9 +2573,9 @@ Base.isequal(A::Transpose{<:Any,<:SparseMatrixCSCMaybeAdjOrTrans}, B::SparseMatr
 
 # Reductions along a dimension return a dense `Array`, as for dense input. A sparse result
 # (issue #43) is opt-in with the keyword `sparse = true`, see `_mapreduce_dim_sparse` below.
-# Covers column-range views so they reduce like their copy (#377), where Base's `similar`
+# Covers views and adjoints so they reduce like their copy (#377), where Base's `similar`
 # would otherwise give a sparse result.
-function Base.reducedim_initarray(A::SparseMatrixCSCUnion, region, v0, ::Type{R}) where {R}
+function Base.reducedim_initarray(A::SparseReducible, region, v0, ::Type{R}) where {R}
     fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
 end
 
@@ -2592,7 +2599,7 @@ function _mapreducezeros(f::F, op::G, ::Type{T}, nzeros::Integer, v0) where {F,G
     v
 end
 
-function Base._mapreduce(f::F, op::G, ::Base.IndexCartesian, A::SparseMatrixCSCUnion{T}) where {F,G,T}
+function Base._mapreduce(f::F, op::G, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
     z = nnz(A)
     n = widelength(A)
     if z == 0
@@ -2617,12 +2624,12 @@ _mapreducezeros(f::Base.ExtremaMap, op::typeof(Base._extrema_rf), ::Type{T}, nze
     nzeros == 0 ? v0 : op(v0, f(zero(T)))
 
 # Specialized mapreduce for any and all
-Base._any(f, A::SparseMatrixCSCUnion, ::Colon) =
+Base._any(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) =
     iszero(widelength(A)) ? false : Base._mapreduce(f, |, IndexCartesian(), A)
-Base._all(f, A::SparseMatrixCSCUnion, ::Colon) =
+Base._all(f, A::SparseMatrixCSCOrColumnSubset, ::Colon) =
     iszero(widelength(A)) ? true  : Base._mapreduce(f, &, IndexCartesian(), A)
 
-function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::SparseMatrixCSCUnion{T}) where {F,T}
+function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Base.IndexCartesian, A::SparseMatrixCSCOrColumnSubset{T}) where {F,T}
     nnzA = nnz(A)
     nzeros = widelength(A) - nnzA
     if nzeros == 0
@@ -2636,26 +2643,26 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas
 end
 
 # The `sparse` keyword, the opt-in for a sparse result (issue #43): `sum(A; dims = 2, sparse = true)`
-# and the other reductions along a dimension return a `SparseMatrixCSC` that stores an entry
+# and the other reductions along a dimension return a `SparseMatrixCSC` (a `SparseVector` for a sparse vector) that stores an entry
 # only for the rows or columns of `A` that store one (all of them when a slice that stores
 # nothing reduces to something nonzero, as for `f(0) != 0`), in time proportional to
 # nnz(A) + length(result) rather than to length(A). Base's `sum`, `prod`, `maximum` and
 # `minimum` forward unknown keywords to `mapreduce`, so this method serves them all; `any`,
 # `all` and `count` do not and get their own methods below.
-function Base.mapreduce(f, op, A::SparseMatrixCSCUnion; dims=:, init=Base._InitialValue(), sparse::Bool=false)
+function Base.mapreduce(f, op, A::SparseReducible; dims=:, init=Base._InitialValue(), sparse::Bool=false)
     sparse || return Base._mapreduce_dim(f, op, init, A, dims)
     dims === (:) && throw(ArgumentError("a sparse result needs a reduction along a dimension, pass `dims`"))
     return _mapreduce_dim_sparse(f, op, init, A, dims)
 end
 for (fname, _fname, op) in ((:any, :_any, :(Base.or_any)), (:all, :_all, :(Base.and_all)))
     @eval begin
-        Base.$fname(A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) = Base.$fname(identity, A; dims, sparse)
-        Base.$fname(f, A::SparseMatrixCSCUnion; dims=:, sparse::Bool=false) =
+        Base.$fname(A::SparseReducible; dims=:, sparse::Bool=false) = Base.$fname(identity, A; dims, sparse)
+        Base.$fname(f, A::SparseReducible; dims=:, sparse::Bool=false) =
             sparse ? mapreduce(f, $op, A; dims, sparse) : Base.$_fname(f, A, dims)
     end
 end
-Base.count(A::SparseMatrixCSCUnion; dims=:, init=0, sparse::Bool=false) = count(identity, A; dims, init, sparse)
-Base.count(f, A::SparseMatrixCSCUnion; dims=:, init=0, sparse::Bool=false) =
+Base.count(A::SparseReducible; dims=:, init=0, sparse::Bool=false) = count(identity, A; dims, init, sparse)
+Base.count(f, A::SparseReducible; dims=:, init=0, sparse::Bool=false) =
     sparse ? mapreduce(Base._bool(f), Base.add_sum, A; dims, init, sparse) : Base._count(f, A, dims, init)
 
 # The slices are reduced the way Base's dense result is: seeded with `init` when given and
@@ -2676,12 +2683,33 @@ _reduced_empty(f, op, init, ::Type{T}) where T = init
 _reduce_unstored(f, op, init, ::Type{T}, len) where T =
     len == 0 ? _reduced_empty(f, op, init, T) : _mapreducezeros(f, op, T, len - 1, _seed(f, op, init, zero(T)))
 
-function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCUnion{T,Ti}, dims) where {T,Ti}
-    m, n = size(A)
-    rm, rn = map(length, Base.reduced_indices(A, dims))   # also validates `dims`
+function _sparse_reduced_eltype(f, op, init, ::Type{T}) where T
     Tr = _reduced_eltype(f, op, init, T)
     applicable(zero, Tr) || throw(ArgumentError("cannot store a sparse result of element type $Tr, " *
         "which has no zero (as for `extrema`); reduce without `sparse = true`"))
+    return Tr
+end
+
+# a slice of `A'` is a slice of `A` along the other dimension
+_mapreduce_dim_sparse(f, op, init, A::Adjoint, dims) =
+    permutedims(_mapreduce_dim_sparse(f ∘ adjoint, op, init, parent(A), map(_switch_dim12, dims)), (2, 1))
+_mapreduce_dim_sparse(f, op, init, A::Transpose, dims) =
+    permutedims(_mapreduce_dim_sparse(f ∘ transpose, op, init, parent(A), map(_switch_dim12, dims)), (2, 1))
+_switch_dim12(d) = d == 1 ? 2 : d == 2 ? 1 : d
+
+function _mapreduce_dim_sparse(f, op, init, A::SparseVectorOrView{T,Ti}, dims) where {T,Ti}
+    Base.reduced_indices(A, dims)   # validates `dims`
+    Tr = _sparse_reduced_eltype(f, op, init, T)
+    # a dimension beyond 1: every entry is a slice of its own
+    1 in dims || return convert(SparseVector{Tr,Ti}, map(x -> _seed(f, op, init, x), A))
+    v = isempty(A) ? _reduced_empty(f, op, init, T) : _seed(identity, op, init, mapreduce(f, op, A))
+    return nnz(A) > 0 || !isequal(v, zero(Tr)) ? SparseVector(1, Ti[1], Tr[v]) : spzeros(Tr, Ti, 1)
+end
+
+function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,Ti}, dims) where {T,Ti}
+    m, n = size(A)
+    rm, rn = map(length, Base.reduced_indices(A, dims))   # also validates `dims`
+    Tr = _sparse_reduced_eltype(f, op, init, T)
     if rm == rn == 1
         R = spzeros(Tr, Ti, 1, 1)
         v = isempty(A) ? _reduced_empty(f, op, init, T) : _seed(identity, op, init, mapreduce(f, op, A))
@@ -2702,8 +2730,7 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCUnion{T,Ti}, dims)
 end
 
 # `R` is a structurally empty `1 x n` sparse matrix: its columns are built in order
-function _mapreducerows_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
-    colptr = getcolptr(A)
+function _mapreducerows_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCOrColumnSubset{T}) where T
     nzval = nonzeros(A)
     m, n = size(A)
     z = zero(eltype(R))
@@ -2712,13 +2739,13 @@ function _mapreducerows_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrix
     zunstored = nnz(A) == m*n && m > 0 ? z : _reduce_unstored(f, op, init, T, m)
     store_unstored = !isequal(zunstored, z)
     Rcolptr, Rrowval, Rnzval = getcolptr(R), rowvals(R), nonzeros(R)
-    nstored = store_unstored ? n : count(col -> colptr[col+1] > colptr[col], 1:n)
+    nstored = store_unstored ? n : count(col -> !isempty(nzrange(A, col)), 1:n)
     resize!(Rrowval, nstored)
     fill!(Rrowval, 1)
     resize!(Rnzval, nstored)
     k = 0
     @inbounds for col in 1:n
-        rng = colptr[col]:colptr[col+1]-1
+        rng = nzrange(A, col)
         if isempty(rng)
             store_unstored || (Rcolptr[col+1] = k + 1; continue)
             v = zunstored
@@ -2739,12 +2766,12 @@ end
 # `R` is a structurally empty `m x 1` sparse matrix. With enough stored entries the rows are
 # reduced in a dense workspace that is then compressed into `R`; a hypersparse `A` instead has
 # its stored entries sorted by row so that only the rows storing something are ever visited.
-function _mapreducecols_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCUnion{T}) where T
+function _mapreducecols_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrixCSCOrColumnSubset{T}) where T
     m, n = size(A)
     Tr = eltype(R)
     z = zero(Tr)
-    rows = view(rowvals(A), _storedrange(A))
-    vals = view(nonzeros(A), _storedrange(A))
+    rows = view(rowvals(A), _storedinds(A))
+    vals = view(nonzeros(A), _storedinds(A))
     nz = length(rows)
     zunstored = nz == m*n && n > 0 ? z : _reduce_unstored(f, op, init, T, n)
     store_unstored = !isequal(zunstored, z)
@@ -2794,31 +2821,29 @@ function _mapreducecols_sparse!(f, op, init, R::SparseMatrixCSC, A::SparseMatrix
 end
 
 # General mapreducedim
-function _mapreducerows!(f, op, R::AbstractArray, A::SparseMatrixCSCUnion{T}) where T
+function _mapreducerows!(f, op, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where T
     require_one_based_indexing(A, R)
-    colptr = getcolptr(A)
     rowval = rowvals(A)
     nzval = nonzeros(A)
     m, n = size(A)
     @inbounds for col in axes(A,2)
         r = R[1, col]
-        @simd for j = colptr[col]:colptr[col+1]-1
+        @simd for j in nzrange(A, col)
             r = op(r, f(nzval[j]))
         end
-        R[1, col] = _mapreducezeros(f, op, T, m-(colptr[col+1]-colptr[col]), r)
+        R[1, col] = _mapreducezeros(f, op, T, m-length(nzrange(A, col)), r)
     end
     R
 end
 
-function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSCUnion{Tv,Ti}) where {Tv,Ti}
+function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{Tv,Ti}) where {Tv,Ti}
     require_one_based_indexing(A, R)
-    colptr = getcolptr(A)
     rowval = rowvals(A)
     nzval = nonzeros(A)
     m, n = size(A)
     rownz = fill(convert(Ti, n), m)
     @inbounds for col in axes(A,2)
-        @simd for j = colptr[col]:colptr[col+1]-1
+        @simd for j in nzrange(A, col)
             row = rowval[j]
             R[row, 1] = op(R[row, 1], f(nzval[j]))
             rownz[row] -= 1
@@ -2830,7 +2855,7 @@ function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSCUnion{Tv,Ti}
     R
 end
 
-function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCUnion{T}) where {F,G,T}
+function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{T}) where {F,G,T}
     require_one_based_indexing(A, R)
     lsiz = Base.check_reducedims(R,A)
     isempty(A) && return R
@@ -2848,13 +2873,12 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCUn
         # Reduction along a dimension > 2
         # Compute op(R, f(A))
         m, n = size(A)
-        colptr = getcolptr(A)
         rowval = rowvals(A)
         nzval = nonzeros(A)
         if nnz(A) == m*n
             # No zeros, so don't compute f(0) since it might throw
             @inbounds for col in axes(A,2)
-                @simd for j = colptr[col]:colptr[col+1]-1
+                @simd for j in nzrange(A, col)
                     R[rowval[j], col] = op(R[rowval[j], col], f(nzval[j]))
                 end
             end
@@ -2862,7 +2886,7 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCUn
             zeroval = f(zero(T))
             @inbounds for col in axes(A,2)
                 lastrow = 0
-                for j = colptr[col]:colptr[col+1]-1
+                for j in nzrange(A, col)
                     row = rowval[j]
                     @simd for i = lastrow+1:row-1 # Zeros before this nonzero
                         R[i, col] = op(R[i, col], zeroval)
@@ -2879,18 +2903,24 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCUn
     R
 end
 
+# LinearAlgebra forwards the commutative reductions of an adjoint or transpose to its parent;
+# this covers the others
+Base._mapreducedim!(f, op, R::AbstractMatrix, A::Adjoint{<:Any,<:SparseMatrixCSCOrColumnSubset}) =
+    (Base._mapreducedim!(f ∘ adjoint, op, PermutedDimsArray(R, (2, 1)), parent(A)); R)
+Base._mapreducedim!(f, op, R::AbstractMatrix, A::Transpose{<:Any,<:SparseMatrixCSCOrColumnSubset}) =
+    (Base._mapreducedim!(f ∘ transpose, op, PermutedDimsArray(R, (2, 1)), parent(A)); R)
+
 # Specialized mapreducedim for + cols to avoid allocating a
 # temporary array when f(0) == 0
-function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCUnion{Tv,Ti}) where {Tv,Ti}
+function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCOrColumnSubset{Tv,Ti}) where {Tv,Ti}
     require_one_based_indexing(A, R)
-    colptr = getcolptr(A)
     rowval = rowvals(A)
     nzval = nonzeros(A)
     m, n = size(A)
     if nnz(A) == m*n
         # No zeros, so don't compute f(0) since it might throw
         @inbounds for col in axes(A,2)
-            @simd for j = colptr[col]:colptr[col+1]-1
+            @simd for j in nzrange(A, col)
                 R[rowval[j], 1] = op(R[rowval[j], 1], f(nzval[j]))
             end
         end
@@ -2899,7 +2929,7 @@ function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCU
         if isequal(zeroval, zero(Tv))
             # Case where f(0) == 0
             @inbounds for col in axes(A,2)
-                @simd for j = colptr[col]:colptr[col+1]-1
+                @simd for j in nzrange(A, col)
                     R[rowval[j], 1] += f(nzval[j])
                 end
             end
@@ -2907,7 +2937,7 @@ function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSCU
             # Case where f(0) != 0
             rownz = fill(convert(Ti, n), m)
             @inbounds for col in axes(A,2)
-                @simd for j = colptr[col]:colptr[col+1]-1
+                @simd for j in nzrange(A, col)
                     row = rowval[j]
                     R[row, 1] += f(nzval[j])
                     rownz[row] -= 1
@@ -2923,14 +2953,13 @@ end
 
 # any(pred, A, dims = 1) => mapreduce(pred, |, A, dims = 1)
 # final argument `post` is to allow post-mapping each columnar mapreduce
-function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::SparseMatrixCSCUnion{Tv},
+function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::SparseMatrixCSCOrColumnSubset{Tv},
                          post::F = identity) where {P, F, Tv}
     nzval = nonzeros(A)
-    colptr = getcolptr(A)
     m, n = size(A)
     @inbounds for ii in axes(A,2)
-        bi, ei = colptr[ii], colptr[ii+1]
-        len = ei - bi
+        rng = nzrange(A, ii)
+        len = length(rng)
         # An empty column is trivial
         if len == 0
             R[1, ii] = post(pred(zero(Tv)))
@@ -2943,7 +2972,7 @@ function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::Spars
         end
         # Otherwise reduce over the stored values
         r = false
-        for jj in bi:(ei - 1)
+        for jj in rng
             r = pred(nzval[jj])
             r && break
         end
@@ -2953,7 +2982,7 @@ function _mapreducerows!(pred::P, ::typeof(|), R::AbstractMatrix{Bool}, A::Spars
 end
 # all(pred, A, dims = 1) => mapreduce(pred, &, A, dims = 1) == .!mapreduce(!pred, |, A, dims = 1)
 _mapreducerows!(pred::P, ::typeof(&), R::AbstractMatrix{Bool},
-                A::SparseMatrixCSCUnion) where {P} = _mapreducerows!(!pred, |, R, A, !)
+                A::SparseMatrixCSCOrColumnSubset) where {P} = _mapreducerows!(!pred, |, R, A, !)
 
 # findmax/min and argmax/min methods
 # find first zero value in sparse matrix - return linear index in full matrix
diff --git a/src/sparsevector.jl b/src/sparsevector.jl
index 70c46899..e1a25863 100644
--- a/src/sparsevector.jl
+++ b/src/sparsevector.jl
@@ -1780,9 +1780,6 @@ for fun in (:+, :-)
 end
 
 ### Reduction
-Base.reducedim_initarray(A::SparseVectorOrView, region, v0, ::Type{R}) where {R} =
-    fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
-
 function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) where {F,G}
     T = eltype(A)
     isempty(A) && return Base.mapreduce_empty(f, op, T)
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index 9269ba1c..6f1e4fe2 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -717,6 +717,30 @@ end
     @test (@which Base._mapreducedim!(identity, +, zeros(10^6, 1), V)).module == SparseArrays
     @test (@which Base._mapreduce(identity, +, IndexCartesian(), V)).module == SparseArrays
     @test nnz(sum(V; dims = 2, sparse = true)) == 2
+    # adjoints, views of a column subset and sparse vectors reduce like their copy, calling `f`
+    # for the stored entries and once per slice rather than per element
+    A, C, v = sprand(60, 50, 0.05), sprand(ComplexF64, 60, 50, 0.05), sprand(60, 0.1)
+    S = view(A, :, [7, 2, 2, 15])
+    for X in (A', transpose(C), C', S, v), dims in (1, 2, (1, 2)),
+        (f, op) in ((abs2, +), (abs, max), (x -> abs(x) + 1, (x, y) -> x + y))   # LinearAlgebra does not forward the last
+        calls = Ref(0)
+        rd = mapreduce(f, op, Array(X); dims, init = 0.0)
+        r = mapreduce(x -> (calls[] += 1; f(x)), op, X; dims, init = 0.0)
+        @test r isa Array && r ≈ rd
+        @test calls[] <= nnz(X) + sum(size(X)) + 1
+        rs = mapreduce(f, op, X; dims, init = 0.0, sparse = true)
+        @test rs isa (X isa AbstractVector ? SparseVector{Float64} : SparseMatrixCSC{Float64}) && rs ≈ rd
+    end
+    for X in (A', S, v), dims in (1, 2)
+        M = Array(X)
+        @test sum(X; dims) isa Array && sum(X; dims) ≈ sum(M; dims)
+        @test prod(X; dims, sparse = true) ≈ prod(M; dims)
+        @test count(!iszero, X; dims, sparse = true) == count(!iszero, M; dims)
+        @test any(!iszero, X; dims, sparse = true) == any(!iszero, M; dims)
+        @test all(iszero, X; dims, sparse = true) == all(iszero, M; dims)
+    end
+    @test sum(S) ≈ sum(Matrix(S)) && prod(x -> x + 1, S) ≈ prod(x -> x + 1, Matrix(S))
+    @test nnz(sum(v; dims = 1, sparse = true)) == 1 && nnz(sum(spzeros(5); dims = 1, sparse = true)) == 0
 end
 
 @testset "oneunit of sparse matrix" begin

From 49692f7fb2131b212303b15223f9d98ac823c739 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Sat, 19 Sep 2026 11:34:44 -0400
Subject: [PATCH 7/7] Fix review findings in the sparse reductions and drop the
 `SparseReducible` alias

Reducing both dimensions of an adjoint folds a copy, so a non-commutative `op` sees
the adjoint's element order. The result's element type comes from a stand-in with the
input's element type, holding a stored value when the input has no zeros, so `Union`
element types keep the dense result's type and `f(0)` is not evaluated for a full
matrix. An empty column range outside the parent has no stored entries. A reduction
along a dimension beyond 2 maps a copy of a view rather than the view element by
element. The keyword methods are generated for the three existing argument types
instead of a new alias.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
---
 src/abstractsparse.jl    |  4 --
 src/sparsematrix.jl      | 96 ++++++++++++++++++++++++----------------
 src/sparsevector.jl      |  3 ++
 test/sparsematrix_ops.jl | 14 ++++++
 4 files changed, 75 insertions(+), 42 deletions(-)

diff --git a/src/abstractsparse.jl b/src/abstractsparse.jl
index b67386f0..aaa2c935 100644
--- a/src/abstractsparse.jl
+++ b/src/abstractsparse.jl
@@ -75,10 +75,6 @@ const SparseVectorPartialView{Tv,Ti} = SubArray{Tv,1,<:AbstractSparseVector{Tv,T
 const SparseMatrixCSCMaybeAdjOrTrans = Union{AbstractSparseMatrixCSC, AdjOrTrans{<:Any,<:AbstractSparseMatrixCSC}}
 const SparseVecOrMatMaybeAdjOrTrans = Union{SparseVecOrMat, AdjOrTrans{<:Any,<:SparseVecOrMat}}
 
-# arguments whose reductions along a dimension take the `sparse` keyword
-const SparseReducible = Union{SparseMatrixCSCOrColumnSubset, AdjOrTrans{<:Any,<:SparseMatrixCSCOrColumnSubset},
-                              SparseVectorOrView}
-
 # LinearAlgebra wrappers around CSC storage
 const SparseTriangular{Tv,Ti} = UpperOrLowerTriangular{Tv,<:SparseMatrixCSCOrView{Tv,Ti}}
 const SparseOrTri{Tv,Ti} = Union{SparseMatrixCSCOrView{Tv,Ti}, SparseTriangular{Tv,Ti}}
diff --git a/src/sparsematrix.jl b/src/sparsematrix.jl
index 6dfa9ed2..6b0aa3c1 100644
--- a/src/sparsematrix.jl
+++ b/src/sparsematrix.jl
@@ -257,7 +257,12 @@ nzvalview(S::AbstractSparseMatrixCSC) = view(nonzeros(S), 1:nnz(S))
 nzvalview(S::SparseMatrixCSCColumnSubset) = view(nonzeros(S), _storedinds(S))
 # where the stored entries sit in the parent's storage: a contiguous range for a column range
 _storedinds(S::AbstractSparseMatrixCSC) = 1:nnz(S)
-_storedinds(S::SparseMatrixCSCView) = (colptr = getcolptr(S); Int(colptr[1]):Int(colptr[end]) - 1)
+function _storedinds(S::SparseMatrixCSCView)
+    cols = S.indices[2]
+    isempty(cols) && return 1:0   # an empty range need not lie within the parent's columns
+    colptr = getcolptr(parent(S))
+    return Int(colptr[first(cols)]):Int(colptr[last(cols)+1]) - 1
+end
 function _storedinds(S::SparseMatrixCSCColumnSubset)
     inds = Int[]
     for col in axes(S, 2)
@@ -2575,7 +2580,8 @@ Base.isequal(A::Transpose{<:Any,<:SparseMatrixCSCMaybeAdjOrTrans}, B::SparseMatr
 # (issue #43) is opt-in with the keyword `sparse = true`, see `_mapreduce_dim_sparse` below.
 # Covers views and adjoints so they reduce like their copy (#377), where Base's `similar`
 # would otherwise give a sparse result.
-function Base.reducedim_initarray(A::SparseReducible, region, v0, ::Type{R}) where {R}
+function Base.reducedim_initarray(A::Union{SparseMatrixCSCOrColumnSubset,AdjOrTrans{<:Any,<:SparseMatrixCSCOrColumnSubset}},
+                                  region, v0, ::Type{R}) where {R}
     fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
 end
 
@@ -2643,38 +2649,45 @@ function Base._mapreduce(f::F, op::Union{typeof(Base.mul_prod),typeof(*)}, ::Bas
 end
 
 # The `sparse` keyword, the opt-in for a sparse result (issue #43): `sum(A; dims = 2, sparse = true)`
-# and the other reductions along a dimension return a `SparseMatrixCSC` (a `SparseVector` for a sparse vector) that stores an entry
-# only for the rows or columns of `A` that store one (all of them when a slice that stores
+# and the other reductions along a dimension return a `SparseMatrixCSC` (a `SparseVector` for a
+# sparse vector) that stores an entry only for the rows or columns of `A` that store one (all of
+# them when a slice that stores
 # nothing reduces to something nonzero, as for `f(0) != 0`), in time proportional to
 # nnz(A) + length(result) rather than to length(A). Base's `sum`, `prod`, `maximum` and
-# `minimum` forward unknown keywords to `mapreduce`, so this method serves them all; `any`,
+# `minimum` forward unknown keywords to `mapreduce`, so the `mapreduce` method serves them all; `any`,
 # `all` and `count` do not and get their own methods below.
-function Base.mapreduce(f, op, A::SparseReducible; dims=:, init=Base._InitialValue(), sparse::Bool=false)
-    sparse || return Base._mapreduce_dim(f, op, init, A, dims)
-    dims === (:) && throw(ArgumentError("a sparse result needs a reduction along a dimension, pass `dims`"))
-    return _mapreduce_dim_sparse(f, op, init, A, dims)
-end
-for (fname, _fname, op) in ((:any, :_any, :(Base.or_any)), (:all, :_all, :(Base.and_all)))
+for T in (:SparseMatrixCSCOrColumnSubset, :(AdjOrTrans{<:Any,<:SparseMatrixCSCOrColumnSubset}), :SparseVectorOrView)
+    @eval function Base.mapreduce(f, op, A::$T; dims=:, init=Base._InitialValue(), sparse::Bool=false)
+        sparse || return Base._mapreduce_dim(f, op, init, A, dims)
+        dims === (:) && throw(ArgumentError("a sparse result needs a reduction along a dimension, pass `dims`"))
+        return _mapreduce_dim_sparse(f, op, init, A, dims)
+    end
+    for (fname, _fname, op) in ((:any, :_any, :(Base.or_any)), (:all, :_all, :(Base.and_all)))
+        @eval begin
+            Base.$fname(A::$T; dims=:, sparse::Bool=false) = Base.$fname(identity, A; dims, sparse)
+            Base.$fname(f, A::$T; dims=:, sparse::Bool=false) =
+                sparse ? mapreduce(f, $op, A; dims, sparse) : Base.$_fname(f, A, dims)
+        end
+    end
     @eval begin
-        Base.$fname(A::SparseReducible; dims=:, sparse::Bool=false) = Base.$fname(identity, A; dims, sparse)
-        Base.$fname(f, A::SparseReducible; dims=:, sparse::Bool=false) =
-            sparse ? mapreduce(f, $op, A; dims, sparse) : Base.$_fname(f, A, dims)
+        Base.count(A::$T; dims=:, init=0, sparse::Bool=false) = count(identity, A; dims, init, sparse)
+        Base.count(f, A::$T; dims=:, init=0, sparse::Bool=false) =
+            sparse ? mapreduce(Base._bool(f), Base.add_sum, A; dims, init, sparse) : Base._count(f, A, dims, init)
     end
 end
-Base.count(A::SparseReducible; dims=:, init=0, sparse::Bool=false) = count(identity, A; dims, init, sparse)
-Base.count(f, A::SparseReducible; dims=:, init=0, sparse::Bool=false) =
-    sparse ? mapreduce(Base._bool(f), Base.add_sum, A; dims, init, sparse) : Base._count(f, A, dims, init)
 
 # The slices are reduced the way Base's dense result is: seeded with `init` when given and
 # otherwise with `mapreduce_first` of their first element, with the entries a slice does not
 # store folded in through `_mapreducezeros`.
 _seed(f, op, ::Base._InitialValue, x) = Base.mapreduce_first(f, op, x)
 _seed(f, op, init, x) = op(init, f(x))
-# element type of the dense result, taken from Base's initialization of a 1 x 1 stand-in,
-# which evaluates f(0) as Base does
-_reduced_eltype(f, op, ::Base._InitialValue, ::Type{T}) where T =
-    eltype(Base.reducedim_init(f, op, fill(zero(T), 1, 1), 1))
-_reduced_eltype(f, op, init, ::Type{T}) where T = typeof(init)
+# element type of the dense result, taken from Base's initialization of a 1 x 1 stand-in with
+# `A`'s element type. Base may map the stand-in's entry, so it is a zero only if `A` has one.
+_reduced_eltype(f, op, ::Base._InitialValue, A::AbstractArray{T}) where T =
+    eltype(Base.reducedim_init(f, op, fill!(Matrix{T}(undef, 1, 1), nnz(A) == length(A) > 0 ? _firststored(A) : zero(T)), 1))
+_reduced_eltype(f, op, init, A) = typeof(init)
+_firststored(A::AbstractVector) = first(nonzeros(A))
+_firststored(A::AbstractMatrix) = nonzeros(A)[first(nzrange(A, 1))]
 # the reduction of a slice with no entries at all, as Base initializes the dense result
 _reduced_empty(f, op, ::Base._InitialValue, ::Type{T}) where T =
     Base.reducedim_init(f, op, Matrix{T}(undef, 0, 1), 1)[1]
@@ -2683,25 +2696,27 @@ _reduced_empty(f, op, init, ::Type{T}) where T = init
 _reduce_unstored(f, op, init, ::Type{T}, len) where T =
     len == 0 ? _reduced_empty(f, op, init, T) : _mapreducezeros(f, op, T, len - 1, _seed(f, op, init, zero(T)))
 
-function _sparse_reduced_eltype(f, op, init, ::Type{T}) where T
-    Tr = _reduced_eltype(f, op, init, T)
+function _sparse_reduced_eltype(f, op, init, A)
+    Tr = _reduced_eltype(f, op, init, A)
     applicable(zero, Tr) || throw(ArgumentError("cannot store a sparse result of element type $Tr, " *
         "which has no zero (as for `extrema`); reduce without `sparse = true`"))
     return Tr
 end
 
-# a slice of `A'` is a slice of `A` along the other dimension
-_mapreduce_dim_sparse(f, op, init, A::Adjoint, dims) =
-    permutedims(_mapreduce_dim_sparse(f ∘ adjoint, op, init, parent(A), map(_switch_dim12, dims)), (2, 1))
-_mapreduce_dim_sparse(f, op, init, A::Transpose, dims) =
-    permutedims(_mapreduce_dim_sparse(f ∘ transpose, op, init, parent(A), map(_switch_dim12, dims)), (2, 1))
+# A slice of `A'` is a slice of `A` along the other dimension. Reducing both dimensions folds
+# the elements in their order, which the parent does not share.
+function _mapreduce_dim_sparse(f, op, init, A::AdjOrTrans, dims)
+    1 in dims && 2 in dims && return _mapreduce_dim_sparse(f, op, init, copy(A), dims)
+    g = A isa Adjoint ? adjoint : transpose
+    return permutedims(_mapreduce_dim_sparse(f ∘ g, op, init, parent(A), map(_switch_dim12, dims)), (2, 1))
+end
 _switch_dim12(d) = d == 1 ? 2 : d == 2 ? 1 : d
 
 function _mapreduce_dim_sparse(f, op, init, A::SparseVectorOrView{T,Ti}, dims) where {T,Ti}
     Base.reduced_indices(A, dims)   # validates `dims`
-    Tr = _sparse_reduced_eltype(f, op, init, T)
+    Tr = _sparse_reduced_eltype(f, op, init, A)
     # a dimension beyond 1: every entry is a slice of its own
-    1 in dims || return convert(SparseVector{Tr,Ti}, map(x -> _seed(f, op, init, x), A))
+    1 in dims || return convert(SparseVector{Tr,Ti}, map(x -> _seed(f, op, init, x), A isa SubArray ? copy(A) : A))
     v = isempty(A) ? _reduced_empty(f, op, init, T) : _seed(identity, op, init, mapreduce(f, op, A))
     return nnz(A) > 0 || !isequal(v, zero(Tr)) ? SparseVector(1, Ti[1], Tr[v]) : spzeros(Tr, Ti, 1)
 end
@@ -2709,7 +2724,7 @@ end
 function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,Ti}, dims) where {T,Ti}
     m, n = size(A)
     rm, rn = map(length, Base.reduced_indices(A, dims))   # also validates `dims`
-    Tr = _sparse_reduced_eltype(f, op, init, T)
+    Tr = _sparse_reduced_eltype(f, op, init, A)
     if rm == rn == 1
         R = spzeros(Tr, Ti, 1, 1)
         v = isempty(A) ? _reduced_empty(f, op, init, T) : _seed(identity, op, init, mapreduce(f, op, A))
@@ -2725,7 +2740,7 @@ function _mapreduce_dim_sparse(f, op, init, A::SparseMatrixCSCOrColumnSubset{T,T
         return _mapreducecols_sparse!(f, op, init, spzeros(Tr, Ti, m, 1), A)
     else
         # a dimension beyond 2: every entry is a slice of its own
-        return convert(SparseMatrixCSC{Tr,Ti}, map(x -> _seed(f, op, init, x), A))
+        return convert(SparseMatrixCSC{Tr,Ti}, map(x -> _seed(f, op, init, x), A isa SubArray ? copy(A) : A))
     end
 end
 
@@ -2904,11 +2919,16 @@ function Base._mapreducedim!(f::F, op::G, R::AbstractArray, A::SparseMatrixCSCOr
 end
 
 # LinearAlgebra forwards the commutative reductions of an adjoint or transpose to its parent;
-# this covers the others
-Base._mapreducedim!(f, op, R::AbstractMatrix, A::Adjoint{<:Any,<:SparseMatrixCSCOrColumnSubset}) =
-    (Base._mapreducedim!(f ∘ adjoint, op, PermutedDimsArray(R, (2, 1)), parent(A)); R)
-Base._mapreducedim!(f, op, R::AbstractMatrix, A::Transpose{<:Any,<:SparseMatrixCSCOrColumnSubset}) =
-    (Base._mapreducedim!(f ∘ transpose, op, PermutedDimsArray(R, (2, 1)), parent(A)); R)
+# this covers the others. Reducing both dimensions folds the elements in their order, which
+# the parent does not share.
+function Base._mapreducedim!(f, op, R::AbstractMatrix, A::AdjOrTrans{<:Any,<:SparseMatrixCSCOrColumnSubset})
+    if size(R, 1) == size(R, 2) == 1
+        Base._mapreducedim!(f, op, R, copy(A))
+    else
+        Base._mapreducedim!(f ∘ (A isa Adjoint ? adjoint : transpose), op, PermutedDimsArray(R, (2, 1)), parent(A))
+    end
+    return R
+end
 
 # Specialized mapreducedim for + cols to avoid allocating a
 # temporary array when f(0) == 0
diff --git a/src/sparsevector.jl b/src/sparsevector.jl
index e1a25863..70c46899 100644
--- a/src/sparsevector.jl
+++ b/src/sparsevector.jl
@@ -1780,6 +1780,9 @@ for fun in (:+, :-)
 end
 
 ### Reduction
+Base.reducedim_initarray(A::SparseVectorOrView, region, v0, ::Type{R}) where {R} =
+    fill!(Array{R}(undef, Base.to_shape(Base.reduced_indices(A, region))), v0)
+
 function Base._mapreduce(f::F, op::G, ::IndexCartesian, A::SparseVectorOrView) where {F,G}
     T = eltype(A)
     isempty(A) && return Base.mapreduce_empty(f, op, T)
diff --git a/test/sparsematrix_ops.jl b/test/sparsematrix_ops.jl
index 6f1e4fe2..026f2696 100644
--- a/test/sparsematrix_ops.jl
+++ b/test/sparsematrix_ops.jl
@@ -741,6 +741,20 @@ end
     end
     @test sum(S) ≈ sum(Matrix(S)) && prod(x -> x + 1, S) ≈ prod(x -> x + 1, Matrix(S))
     @test nnz(sum(v; dims = 1, sparse = true)) == 1 && nnz(sum(spzeros(5); dims = 1, sparse = true)) == 0
+    # reducing both dimensions of an adjoint keeps its element order for a non-commutative `op`
+    firstnz(x, y) = iszero(x) ? y : x
+    B = sparse([0 1; 2 0])
+    @test mapreduce(identity, firstnz, B'; dims = (1, 2), init = 0) == [1;;] == mapreduce(identity, firstnz, B'; dims = (1, 2), init = 0, sparse = true)
+    # the element type of the dense result for a `Union`, and no f(0) for a full matrix
+    @test sum(sparse(Union{Int,Float64}[1.5 2; 3 4]); dims = 1, sparse = true) == [4.5 6.0]
+    @test maximum(x -> 1 ÷ x, sparse([1 2; 3 4]); dims = 1, sparse = true) == [1 0]
+    # an empty column range outside the parent
+    V = view(spzeros(4, 5), :, 10:9)
+    @test nnz(V) == 0 && sum(V) == 0 && size(sum(V; dims = 1, sparse = true)) == (1, 0)
+    # a dimension beyond 2 maps the stored entries of a view only
+    calls = Ref(0)
+    @test mapreduce(x -> (calls[] += 1; x), +, view(A, :, [7, 2]); dims = 3, sparse = true) == A[:, [7, 2]]
+    @test calls[] <= nnz(A) + 1
 end
 
 @testset "oneunit of sparse matrix" begin