Fix size-1 StructuredMatrix's broadcast.#54190
Conversation
1. size-1 StructuredMatrix should behave like scalar during broadcast. Thus their `fzero` should return the only element. (fix #54087) 2. But for simple broadcast with only one StructuredMatrix, we can keep stability as the structure is "preserved" even for size-1 case. Thus `count_structedmatrix` is added to check that. 3. `count_structedmatrix` is fused to keep `Bidiagonal` stability. (replace JuliaLang#54067)
|
The impact of the type-instability is small on performance julia> D = Diagonal(ones(10000));
julia> @btime $D .+ $D;
7.359 μs (5 allocations: 78.22 KiB)
julia> @btime $D + $D;
7.330 μs (3 allocations: 78.19 KiB)but I wonder if this might lead to cascading type-instabilities downstream if one tries to store the result in a struct/array. The single-element structured matrix is a rather niche use case as such. I wonder if this should just be deprecated, with a suggestion to use a |
I agree the size-1 case has no realistic usage, but throwing for inputs which are valid for general fallback looks too broken.
I guess there's no need to worry too much for this. julia> bad(T) = T .* 2
bad (generic function with 1 method)
julia> good(T) = T .* 2.
good (generic function with 1 method)
julia> T = Tridiagonal(1:1,1:2,1:1)
2×2 Tridiagonal{Int64, UnitRange{Int64}}:
1 1
1 2
julia> @code_warntype good(T)
MethodInstance for good(::Tridiagonal{Int64, UnitRange{Int64}})
from good(T) @ Main REPL[2]:1
Arguments
#self#::Core.Const(good)
T::Tridiagonal{Int64, UnitRange{Int64}}
Body::Tridiagonal{Float64, Vector{Float64}}
1 ─ %1 = Main.:*::Core.Const(*)
│ %2 = Base.broadcasted(%1, T, 2.0)::Core.PartialStruct(Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Tridiagonal}, Nothing, typeof(*), Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Float64}}, Any[Core.Const(LinearAlgebra.StructuredMatrixStyle{Tridiagonal}()), Core.Const(*), Core.PartialStruct(Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Float64}, Any[Tridiagonal{Int64, UnitRange{Int64}}, Core.Const(2.0)]), Nothing])
│ %3 = Base.materialize(%2)::Tridiagonal{Float64, Vector{Float64}}
└── return %3
julia> @code_warntype bad(T)
MethodInstance for bad(::Tridiagonal{Int64, UnitRange{Int64}})
from bad(T) @ Main REPL[1]:1
Arguments
#self#::Core.Const(bad)
T::Tridiagonal{Int64, UnitRange{Int64}}
Body::Union{Tridiagonal{Int64, Vector{Int64}}, Matrix{Int64}}
1 ─ %1 = Main.:*::Core.Const(*)
│ %2 = Base.broadcasted(%1, T, 2)::Core.PartialStruct(Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{Tridiagonal}, Nothing, typeof(*), Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Int64}}, Any[Core.Const(LinearAlgebra.StructuredMatrixStyle{Tridiagonal}()), Core.Const(*), Core.PartialStruct(Tuple{Tridiagonal{Int64, UnitRange{Int64}}, Int64}, Any[Tridiagonal{Int64, UnitRange{Int64}}, Core.Const(2)]), Nothing])
│ %3 = Base.materialize(%2)::Union{Tridiagonal{Int64, Vector{Int64}}, Matrix{Int64}}
└── return %3Users can always use typeassert to restore stability where possible. |
|
I'm not really in favor of making this type-unstable to address a rarely encountered case, but I wonder what @mbauman might have to say on this |
|
We have moved the LinearAlgebra stdlib to an external repo: https://github.com/JuliaLang/LinearAlgebra.jl @N5N3 If you think that this PR is still relevant, please open a new PR on the LinearAlgebra.jl repo. |
size-1 StructuredMatrix should behave like scalar during broadcast. Thus their
fzeroshould return the only element.But for simple broadcast with only one StructuredMatrix, we can keep stability as the structure is "preserved" even for size-1 case. Thus
count_structedmatrixis added to check that.count_structedmatrixis fused to keepBidiagonalstability.fix JuliaLang/LinearAlgebra.jl#1063
replace #54067