From f0ff81837db3a0775076374c227ded8d576d1e29 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 15 Oct 2026 17:03:48 +0200 Subject: [PATCH 1/6] Prepare contextual zeros for reductions and array products Add zero_for(T, inputs...) so scalar types whose additive identity depends on operand values can prepare their context before accumulation. The default retains type-based zeros without inspecting the inputs. Use the hook for sums, fused reductions, array output initialization, and sparse-to-dense conversion. Inspect stored sparse entries and preserve transpose/adjoint transformations. Give output entries independent zeros and reuse the existing array multiplication paths. Check reductions, output initialization, input ownership, inference, and allocation budgets. --- src/dispatch.jl | 24 +++++-- src/implementations/LinearAlgebra.jl | 22 +++++-- src/implementations/SparseArrays.jl | 20 ++++++ src/reduce.jl | 43 +++++++++++- test/test_matmul.jl | 53 +++++++++++++++ test/test_reduce.jl | 97 ++++++++++++++++++++++++++++ test/test_sparse_arrays.jl | 20 ++++++ 7 files changed, 264 insertions(+), 15 deletions(-) create mode 100644 test/test_reduce.jl diff --git a/src/dispatch.jl b/src/dispatch.jl index 9737343..2c01ccf 100644 --- a/src/dispatch.jl +++ b/src/dispatch.jl @@ -33,15 +33,25 @@ end function Base.sum( a::AbstractArray{T}; dims = :, - init = zero(promote_operation(+, T, T)), + init = nothing, ) where {T<:AbstractMutable} + accumulator = if isnothing(init) + zero_for(promote_operation(+, T, T), _zero_entries(a)) + elseif isempty(a) + copy_if_mutable(init) + else + # Base.sum owns its accumulator: include the seed's context, then + # add it without mutating the caller's value. + R = promote_operation(+, T, T) + operate(+, init, zero_for(R, _zero_entries(a), (init,))) + end if dims !== Colon() - # We cannot use `mapreduce` with `add!!` instead of `Base.add_mul` like + # We cannot use `mapreduce` with `add!!` instead of `Base.add_sum` like # `operate(sum, ...)` because the same instance given at `init` is used # at several places. - return mapreduce(identity, Base.add_sum, a; dims, init) + return mapreduce(identity, Base.add_sum, a; dims, init = accumulator) end - return operate(sum, a; init) + return operate(sum, a; init = accumulator) end # When doing `x'y` where the elements of `x` and/or `y` are arrays, redirecting @@ -120,7 +130,7 @@ function _mul!(output, A, B, α, β) # See SparseArrays/src/linalg.jl if !isone(β) if iszero(β) - operate!(zero, output) + _zero_for!(output, A, B) else LinearAlgebra.rmul!(output, scaling(β)) end @@ -129,7 +139,7 @@ function _mul!(output, A, B, α, β) end function _mul!(output, A, B, α) - operate!(zero, output) + _zero_for!(output, A, B) return operate!(add_mul, output, A, B, scaling(α)) end @@ -470,7 +480,7 @@ end function Base.Matrix(S::_SparseMat{T}) where {T<:AbstractMutable} U = promote_operation(+, promote_operation(zero, T), T) A = Matrix{U}(undef, size(S)...) - operate!(zero, A) + _zero_for!(A, S) return operate!(+, A, S) end diff --git a/src/implementations/LinearAlgebra.jl b/src/implementations/LinearAlgebra.jl index 0f81c4d..f9322f7 100644 --- a/src/implementations/LinearAlgebra.jl +++ b/src/implementations/LinearAlgebra.jl @@ -304,9 +304,9 @@ function _dim_check(C::AbstractVector, A::AbstractMatrix, B::AbstractVector) return end -function _dim_check(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix) +function _dim_check(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractMatrix) mB, nB = size(B) - mA, nA = size(A) + mA, nA = size(A, 1), size(A, 2) if mB != nA throw( DimensionMismatch( @@ -383,12 +383,20 @@ function operate!( return buffered_operate!(buffer, add_mul, C, A, B) end -function operate!(::typeof(zero), C::Union{Vector,Matrix}) +function _zero_for!(C::AbstractArray, inputs::Vararg{Any,N}) where {N} + isempty(C) && return C + # Prepare context once, then give each output entry its own zero. + z = zero_for(eltype(C), _zero_inputs(*, inputs...)...) # C may contain undefined values so we cannot call `zero!` + first_index = first(eachindex(C)) for i in eachindex(C) - @inbounds C[i] = zero(eltype(C)) + @inbounds C[i] = i == first_index ? z : zero(z) end - return + return C +end + +function operate!(::typeof(zero), C::AbstractArray) + return _zero_for!(C) end function operate_to!( @@ -397,7 +405,9 @@ function operate_to!( A::AbstractArray, B::AbstractArray, ) - operate!(zero, C) + _dim_check(C, A, B) + isempty(C) && return C + _zero_for!(C, A, B) return operate!(add_mul, C, A, B) end diff --git a/src/implementations/SparseArrays.jl b/src/implementations/SparseArrays.jl index bf36578..c2b5f15 100644 --- a/src/implementations/SparseArrays.jl +++ b/src/implementations/SparseArrays.jl @@ -9,6 +9,18 @@ const _SparseMat = SparseArrays.SparseMatrixCSC +function _zero_entries(input::SparseArrays.AbstractSparseArray) + return SparseArrays.nonzeros(input) +end +for (Wrapper, op) in + ((LinearAlgebra.Transpose, transpose), (LinearAlgebra.Adjoint, adjoint)) + @eval function _zero_entries( + input::$Wrapper{<:Any,<:SparseArrays.AbstractSparseArray}, + ) + return Iterators.map($op, SparseArrays.nonzeros(parent(input))) + end +end + function mutable_copy(v::SparseArrays.SparseVector) result = copy(v) vals = SparseArrays.nonzeros(result) @@ -39,6 +51,14 @@ function operate!(::typeof(zero), A::_SparseMat) return A end +function _zero_for!( + C::Union{SparseArrays.SparseVector,_SparseMat}, + inputs::Vararg{Any,N}, +) where {N} + # Sparse outputs are cleared structurally; there are no entries to fill. + return operate!(zero, C) +end + function promote_operation( op::Union{typeof(+),typeof(-)}, ::Type{<:SparseArrays.AbstractSparseArray{Tv,Ti,N}}, diff --git a/src/reduce.jl b/src/reduce.jl index fbd4c9a..f0b50c4 100644 --- a/src/reduce.jl +++ b/src/reduce.jl @@ -25,6 +25,27 @@ reduce_op(::typeof(add_dot)) = + neutral_element(::typeof(+), ::Type{T}) where {T} = zero(T) neutral_element(::typeof(+), ::Type{T}) where {T<:AbstractArray} = Zero() +""" + zero_for(::Type{T}, inputs...) + +Return a fresh additive identity for accumulating values from `inputs` with +promoted result type `T`. The inputs are iterables of scalar operands and must +not be mutated. The returned value can be mutated independently of the inputs +and of zeros returned by other calls. + +The default returns `zero(T)` without inspecting the inputs. +For array-valued results, it retains the symbolic [`Zero`](@ref) identity. +Types whose zeros need value information can specialize this function to scan +the inputs before constructing the zero. +For example, algebra elements can prepare a common parent for all operands. + +If empty inputs provide no context, a specialization can use `zero(T)` when +available. Otherwise it must report that no zero can be constructed. +""" +function zero_for(::Type{T}, inputs::Vararg{Any,N}) where {T,N} + return neutral_element(+, T) +end + map_op(::AddSubMul) = * map_op(::typeof(add_dot)) = LinearAlgebra.dot @@ -39,10 +60,25 @@ end _concrete_eltype(x) = isempty(x) ? eltype(x) : typeof(first(x)) +# Implicit sparse zeros add no context and need not be constructed or indexed. +_zero_entries(input) = input +_zero_inputs(op, args::Vararg{Any,N}) where {N} = map(_zero_entries, args) +# Prepare from the factors used by add_dot: adjoint can change their context. +function _zero_inputs(::typeof(add_dot), a, b) + return (Iterators.map(adjoint, _zero_entries(a)), _zero_entries(b)) +end + +""" + fused_map_reduce(op, args...) + +Reduce corresponding entries of `args` using the fused operation `op`. +The accumulator is prepared by [`zero_for`](@ref) from the promoted result type +and the inputs (with the left operands adjointed for `add_dot`). +""" function fused_map_reduce(op::F, args::Vararg{Any,N}) where {F<:Function,N} _check_same_length(args...) T = promote_map_reduce(op, _concrete_eltype.(args)...) - accumulator = neutral_element(reduce_op(op), T) + accumulator = zero_for(T, _zero_inputs(op, args...)...) buffer = buffer_for(op, T, eltype.(args)...) for I in zip(eachindex.(args)...) accumulator = @@ -54,7 +90,10 @@ end function operate( ::typeof(sum), a::AbstractArray; - init = zero(promote_operation(+, eltype(a), eltype(a))), + init = zero_for( + promote_operation(+, eltype(a), eltype(a)), + _zero_entries(a), + ), ) return mapreduce(identity, add!!, a; init) end diff --git a/test/test_matmul.jl b/test/test_matmul.jl index 2dce391..98d5a17 100644 --- a/test/test_matmul.jl +++ b/test/test_matmul.jl @@ -283,6 +283,59 @@ function test_matrix_matrix_product() return end +function test_product_zero_allocations() + for T in (Int, Float64, BigInt), matrix in (false, true) + A = T[1 2; 3 4] + B = matrix ? T[2 3; 4 5] : T[2, 3] + C = matrix ? Matrix{T}(undef, 2, 2) : Vector{T}(undef, 2) + @test (@inferred MA.operate_to!(C, *, A, B)) === C + expected = copy(C) + # The existing type-only initialization is the allocation baseline. + MA.operate!(zero, C) + MA.operate!(MA.add_mul, C, A, B) + budget = @allocated begin + MA.operate!(zero, C) + MA.operate!(MA.add_mul, C, A, B) + end + @test (@allocated MA.operate_to!(C, *, A, B)) <= budget + @test C == expected + bad = matrix ? B[1:1, :] : B[1:1] + @test_throws DimensionMismatch MA.operate_to!(C, *, A, bad) + @test C == expected + end + C = Matrix{BigInt}(undef, 2, 2) + @test MA.operate_to!( + C, + *, + zeros(Int, 2, 0), + Matrix{BigInt}(undef, 0, 2), + ) === C + MA.operate!(+, C[1], 1) + @test C == [1 0; 0 0] + return +end + +function test_product_array_outputs() + for T in (Int, BigInt) + A = T[1 2; 3 4] + b = T[2, 3] + for (left, right) in ((A, b), (A, hcat(b, b)), (b, transpose(b))) + expected = left * right + C = similar(expected) + @test (@inferred MA.operate_to!(C, *, left, right)) === C + @test C == expected + output_view = view(C, axes(C)...) + @test (@inferred MA.operate!(zero, output_view)) === output_view + @test all(iszero, C) + if T === BigInt + MA.operate!(+, C[1], 1) + @test all(iszero, C[2:end]) + end + end + end + return +end + function test_matrix_multiplication() X = ones(BigInt, 1, 1) M = ones(1, 1) diff --git a/test/test_reduce.jl b/test/test_reduce.jl new file mode 100644 index 0000000..05f036d --- /dev/null +++ b/test/test_reduce.jl @@ -0,0 +1,97 @@ +# Copyright (c) 2026 MutableArithmetics.jl contributors +# +# This Source Code Form is subject to the terms of the Mozilla Public License, +# v.2.0. If a copy of the MPL was not distributed with this file, You can obtain +# one at http://mozilla.org/MPL/2.0/. + +module TestReduce + +using Test +import LinearAlgebra +import MutableArithmetics as MA +import SparseArrays + +struct UnreadableInput end +function Base.iterate(::UnreadableInput) + return error("The type-based zero must not read inputs") +end + +@testset "Zeros from inputs" begin + for T in (Int, Float64, BigInt) + a = @inferred MA.zero_for(T, UnreadableInput()) + b = @inferred MA.zero_for(T, UnreadableInput()) + @test iszero(a) && iszero(b) + @test typeof(a) === T + if T === BigInt + @test a !== b + MA.operate!(+, a, 1) + @test iszero(b) + end + @test (@inferred MA.operate(sum, T[1, 2, 3])) == 6 + @test iszero(MA.operate(sum, T[])) + end + @test MA.zero_for(Matrix{Int}, UnreadableInput()) isa MA.Zero +end + +@testset "Sparse reductions preserve adjoints and inputs" begin + for values in ([1, 0, 2], BigInt[1, 0, 2], [1 + 2im, 0, 3 - im]) + a = SparseArrays.sparse(values) + b = [2, 3, 4] + @test (@inferred MA.operate(sum, a)) == sum(values) + @test (@inferred MA.operate(sum, transpose(a))) == sum(values) + @test (@inferred MA.operate(sum, adjoint(a))) == conj(sum(values)) + @test (@inferred MA.fused_map_reduce(MA.add_dot, a, b)) == + LinearAlgebra.dot(values, b) + @test a == values + end +end + +@testset "Fused reductions" begin + for T in (Int, BigInt), op in (MA.add_mul, MA.add_dot) + a, b = T[1, 2, 3], T[4, 5, 6] + originals = MA.mutable_copy.((a, b)) + @test (@inferred MA.fused_map_reduce(op, a, b)) == 32 + @test (@inferred MA.fused_map_reduce(MA.sub_mul, a, b)) == -32 + @test iszero(MA.fused_map_reduce(op, T[], T[])) + @test_throws DimensionMismatch MA.fused_map_reduce(op, a, b[1:1]) + @test (a, b) == originals + end + a, b = [1 + 2im, 3 - im], [2 - im, 1 + im] + @test (@inferred MA.fused_map_reduce(MA.add_dot, a, b)) == + LinearAlgebra.dot(a, b) + @test MA.fused_map_reduce(MA.add_mul, Matrix{Int}[], Float64[]) isa MA.Zero +end + +function test_allocations() + a, b = [1, 2, 3], [4, 5, 6] + @test (@inferred MA.fused_map_reduce(MA.add_mul, a, b)) == 32 + @test (@allocated MA.fused_map_reduce(MA.add_mul, a, b)) == 0 + @test (@inferred MA.fused_map_reduce(MA.add_dot, a, b)) == 32 + @test (@allocated MA.fused_map_reduce(MA.add_dot, a, b)) == 0 + @test (@inferred MA.operate(sum, a)) == 6 + @test (@allocated MA.operate(sum, a)) == 0 + @test (@inferred MA.zero_for(Int, UnreadableInput())) == 0 + @test (@allocated MA.zero_for(Int, UnreadableInput())) == 0 + return +end + +@testset "Fused reduction allocations" begin + test_allocations() +end + +function test_bigint_allocations() + a, b = BigInt[1, 2, 3], BigInt[4, 5, 6] + MA.fused_map_reduce(MA.add_mul, a, b) + MA.fused_map_reduce(MA.add_dot, a, b) + # Preserve the existing budget: result plus multiplication workspace. + budget = 3 * @allocated(BigInt(1)) + @test (@allocated MA.fused_map_reduce(MA.add_mul, a, b)) <= budget + @test (@allocated MA.fused_map_reduce(MA.add_dot, a, b)) <= budget + return +end + +@testset "BigInt reduction allocations" begin + test_bigint_allocations() +end + +end diff --git a/test/test_sparse_arrays.jl b/test/test_sparse_arrays.jl index 0f96b84..1391e2c 100644 --- a/test/test_sparse_arrays.jl +++ b/test/test_sparse_arrays.jl @@ -72,6 +72,26 @@ function test_spmatmul() return end +function test_operate_to_product() + for T in (Int, BigInt), + sparse_inputs in (false, true), + matrix in (false, true) + + a = T[1 0 2; 0 3 0] + b = matrix ? T[0 4; 5 0; 0 6] : T[4, 5, 6] + A, B = sparse_inputs ? SparseArrays.sparse.((a, b)) : (a, b) + C = SparseArrays.sparse(matrix ? T[7 8; 9 10] : T[7, 8]) + @test (@inferred MA.operate_to!(C, *, A, B)) === C + @test C == a * b + bad = matrix ? B[1:2, :] : B[1:2] + @test_throws DimensionMismatch MA.operate_to!(C, *, A, bad) + @test C == a * b + @test MA.operate_to!(C, *, A, zero(B)) === C + @test iszero(SparseArrays.nnz(C)) + end + return +end + function test_spmatmul_prefer_sort() Random.seed!(1234) m = n = 100 From 56a45a037a969688cae50bf0c43fd4fe3a1c79ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 15 Oct 2026 16:42:37 +0200 Subject: [PATCH 2/6] ae --- src/reduce.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/reduce.jl b/src/reduce.jl index f0b50c4..f002988 100644 --- a/src/reduce.jl +++ b/src/reduce.jl @@ -37,7 +37,8 @@ The default returns `zero(T)` without inspecting the inputs. For array-valued results, it retains the symbolic [`Zero`](@ref) identity. Types whose zeros need value information can specialize this function to scan the inputs before constructing the zero. -For example, algebra elements can prepare a common parent for all operands. +For example, `StarAlgebras.AlgebraElement`s can prepare a common parent for all +operands. If empty inputs provide no context, a specialization can use `zero(T)` when available. Otherwise it must report that no zero can be constructed. From d667a731138b385b841516f91462fecafe70faae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 15 Oct 2026 22:03:58 +0200 Subject: [PATCH 3/6] Reuse contextual sum accumulators and fix zero_for docs Add the seed into the prepared zero without allocating a second result. Preserve input ownership and include the seed type when choosing the accumulator. Measure seeded sums with a preallocated seed, and check inference and input preservation. Replace the broken documentation reference to the symbolic Zero identity. --- src/dispatch.jl | 4 ++-- src/reduce.jl | 2 +- test/test_basics.jl | 9 +++++++-- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/src/dispatch.jl b/src/dispatch.jl index 2c01ccf..84e4ef5 100644 --- a/src/dispatch.jl +++ b/src/dispatch.jl @@ -42,8 +42,8 @@ function Base.sum( else # Base.sum owns its accumulator: include the seed's context, then # add it without mutating the caller's value. - R = promote_operation(+, T, T) - operate(+, init, zero_for(R, _zero_entries(a), (init,))) + R = promote_operation(+, typeof(init), promote_operation(+, T, T)) + operate!!(+, zero_for(R, _zero_entries(a), (init,)), init) end if dims !== Colon() # We cannot use `mapreduce` with `add!!` instead of `Base.add_sum` like diff --git a/src/reduce.jl b/src/reduce.jl index f002988..19ecc2f 100644 --- a/src/reduce.jl +++ b/src/reduce.jl @@ -34,7 +34,7 @@ not be mutated. The returned value can be mutated independently of the inputs and of zeros returned by other calls. The default returns `zero(T)` without inspecting the inputs. -For array-valued results, it retains the symbolic [`Zero`](@ref) identity. +For array-valued results, it retains the symbolic `Zero` identity. Types whose zeros need value information can specialize this function to scan the inputs before constructing the zero. For example, `StarAlgebras.AlgebraElement`s can prepare a common parent for all diff --git a/test/test_basics.jl b/test/test_basics.jl index 2aa5af5..408b759 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -178,13 +178,18 @@ end function test_sum_with_init() x = convert(Vector{DummyBigInt}, 1:100) + init = DummyBigInt(7) + @test MA.isequal_canonical(@inferred(sum(x; init)), DummyBigInt(5057)) + @test MA.isequal_canonical(init, DummyBigInt(7)) + @test MA.isequal_canonical(x, convert(Vector{DummyBigInt}, 1:100)) # compilation @allocated sum(x) - @allocated sum(x; init = DummyBigInt(0)) + @allocated sum(x; init) @allocated non_mutable_sum_pr306(x) # now test actual allocations no_init = @allocated sum(x) - with_init = @allocated sum(x; init = DummyBigInt(0)) + # Both sums own an accumulator; exclude construction of the caller's seed. + with_init = @allocated sum(x; init) no_ma = @allocated non_mutable_sum_pr306(x) # There's an additional 16 bytes for kwarg version. Upper bound by 40 to be # safe between Julia versions From da3f147b50bfeab0a30bb348b37780105edfd506 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 16 Oct 2026 08:55:21 +0200 Subject: [PATCH 4/6] Prepare zeros directly from operand entries Remove _zero_inputs and its add_dot-specific adjoint mapping. Pass the original operand entries to zero_for in reductions and array output initialization, retaining sparse-entry filtering. The reduction and matrix-product test files pass all 390 assertions. Complex-polynomial dot products also preserve their result and inputs. --- src/implementations/LinearAlgebra.jl | 2 +- src/reduce.jl | 9 ++------- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/src/implementations/LinearAlgebra.jl b/src/implementations/LinearAlgebra.jl index f9322f7..c7d8691 100644 --- a/src/implementations/LinearAlgebra.jl +++ b/src/implementations/LinearAlgebra.jl @@ -386,7 +386,7 @@ end function _zero_for!(C::AbstractArray, inputs::Vararg{Any,N}) where {N} isempty(C) && return C # Prepare context once, then give each output entry its own zero. - z = zero_for(eltype(C), _zero_inputs(*, inputs...)...) + z = zero_for(eltype(C), map(_zero_entries, inputs)...) # C may contain undefined values so we cannot call `zero!` first_index = first(eachindex(C)) for i in eachindex(C) diff --git a/src/reduce.jl b/src/reduce.jl index 19ecc2f..54c6e8e 100644 --- a/src/reduce.jl +++ b/src/reduce.jl @@ -63,23 +63,18 @@ _concrete_eltype(x) = isempty(x) ? eltype(x) : typeof(first(x)) # Implicit sparse zeros add no context and need not be constructed or indexed. _zero_entries(input) = input -_zero_inputs(op, args::Vararg{Any,N}) where {N} = map(_zero_entries, args) -# Prepare from the factors used by add_dot: adjoint can change their context. -function _zero_inputs(::typeof(add_dot), a, b) - return (Iterators.map(adjoint, _zero_entries(a)), _zero_entries(b)) -end """ fused_map_reduce(op, args...) Reduce corresponding entries of `args` using the fused operation `op`. The accumulator is prepared by [`zero_for`](@ref) from the promoted result type -and the inputs (with the left operands adjointed for `add_dot`). +and the inputs. """ function fused_map_reduce(op::F, args::Vararg{Any,N}) where {F<:Function,N} _check_same_length(args...) T = promote_map_reduce(op, _concrete_eltype.(args)...) - accumulator = zero_for(T, _zero_inputs(op, args...)...) + accumulator = zero_for(T, map(_zero_entries, args)...) buffer = buffer_for(op, T, eltype.(args)...) for I in zip(eachindex.(args)...) accumulator = From 4dc86dff349a03724cbf63648c8801a530da2704 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 16 Oct 2026 09:18:45 +0200 Subject: [PATCH 5/6] Prepare contextual zeros from stored diagonal entries Read Diagonal.diag recursively to avoid constructing off-diagonal zeros and implicit zeros in sparse diagonal storage. Parent-dependent element types can then prepare an accumulator from the stored values. Add a regression ensuring the default type-based zero does not read diagonal entries. The reduction tests pass 65 assertions, and 47 focused JuMP rewrite, expression and constraint checks pass. --- src/implementations/LinearAlgebra.jl | 4 ++++ test/test_reduce.jl | 8 ++++++++ 2 files changed, 12 insertions(+) diff --git a/src/implementations/LinearAlgebra.jl b/src/implementations/LinearAlgebra.jl index c7d8691..f1f3852 100644 --- a/src/implementations/LinearAlgebra.jl +++ b/src/implementations/LinearAlgebra.jl @@ -383,6 +383,10 @@ function operate!( return buffered_operate!(buffer, add_mul, C, A, B) end +function _zero_entries(input::LinearAlgebra.Diagonal) + return _zero_entries(input.diag) +end + function _zero_for!(C::AbstractArray, inputs::Vararg{Any,N}) where {N} isempty(C) && return C # Prepare context once, then give each output entry its own zero. diff --git a/test/test_reduce.jl b/test/test_reduce.jl index 05f036d..b4c7df8 100644 --- a/test/test_reduce.jl +++ b/test/test_reduce.jl @@ -16,6 +16,12 @@ function Base.iterate(::UnreadableInput) return error("The type-based zero must not read inputs") end +struct UnreadableVector{T} <: AbstractVector{T} end +Base.size(::UnreadableVector) = (2,) +function Base.getindex(::UnreadableVector, ::Int) + return error("The type-based zero must not read diagonal entries") +end + @testset "Zeros from inputs" begin for T in (Int, Float64, BigInt) a = @inferred MA.zero_for(T, UnreadableInput()) @@ -29,6 +35,8 @@ end end @test (@inferred MA.operate(sum, T[1, 2, 3])) == 6 @test iszero(MA.operate(sum, T[])) + d = LinearAlgebra.Diagonal(UnreadableVector{T}()) + @test iszero(@inferred MA.zero_for(T, MA._zero_entries(d))) end @test MA.zero_for(Matrix{Int}, UnreadableInput()) isa MA.Zero end From 59e30b2f543fa20de5e31f7cec051b23a70b546e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 16 Oct 2026 09:27:35 +0200 Subject: [PATCH 6/6] Explain the unreadable input tests for zero preparation --- test/test_reduce.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_reduce.jl b/test/test_reduce.jl index b4c7df8..bac218d 100644 --- a/test/test_reduce.jl +++ b/test/test_reduce.jl @@ -11,6 +11,10 @@ import LinearAlgebra import MutableArithmetics as MA import SparseArrays +# StarAlgebras needs to read entries to determine the parent of its zero. +# Other types, such as JuMP expressions, can construct zero from the type alone. +# These fixtures ensure that the default zero_for path, including _zero_entries, +# does not read inputs, preserving that fast path even for structured arrays. struct UnreadableInput end function Base.iterate(::UnreadableInput) return error("The type-based zero must not read inputs")