Preserve structured matrices for zero-absorbing operations#1549
Preserve structured matrices for zero-absorbing operations#1549max-vassili3v wants to merge 8 commits intoJuliaLang:masterfrom
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## master #1549 +/- ##
==========================================
+ Coverage 94.03% 94.07% +0.03%
==========================================
Files 35 35
Lines 15986 16001 +15
==========================================
+ Hits 15033 15053 +20
+ Misses 953 948 -5 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
Is the untested function definition actually needed? |
src/structuredbroadcast.jl
Outdated
| # broadcasts with SymTridiagonal and arrays with more than one entry will generally break symmetry. | ||
| # This is useful for knowing whether to materialize a Tridiagonal for zero-preserving functions. | ||
| function issymmetrybreaking(bc::Broadcasted{StructuredMatrixStyle{SymTridiagonal}}) | ||
| any(x -> !(x isa SymTridiagonal || all(y -> length(y) == 1, axes(x))), bc.args) |
There was a problem hiding this comment.
Is there much of an advantage in preserving the structure for 1x1 matrices here, given that this comes at the price of type stability?
There was a problem hiding this comment.
I made it so that 1x1 SymTridiagonals always become Tridiagonals, and have changed the unit tests accordingly
|
I noticed that the following seems to have broken on this PR but is working on master: julia> using LinearAlgebra
julia> D = Diagonal(1:4)
4×4 Diagonal{Int64, UnitRange{Int64}}:
1 ⋅ ⋅ ⋅
⋅ 2 ⋅ ⋅
⋅ ⋅ 3 ⋅
⋅ ⋅ ⋅ 4
julia> (1:4) .* D .* (1:4)'
ERROR: MethodError: no method matching diagview(::UnitRange{Int64}, ::Int64)
The function `diagview` exists, but no method is defined for this combination of argument types.
Closest candidates are:
diagview(::Adjoint, ::Integer)
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/adjtrans.jl:579
diagview(::Transpose, ::Integer)
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/adjtrans.jl:578
diagview(::AbstractMatrix, ::Integer)
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/dense.jl:358
...
Stacktrace:
[1] diagview(A::Adjoint{Int64, UnitRange{Int64}}, k::Int64)
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/adjtrans.jl:579
[2] _preprocess_broadcasted(::Type{Diagonal}, d::Adjoint{Int64, UnitRange{Int64}})
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:286
[3] preprocess_broadcasted(::Type{Diagonal}, A::Adjoint{Int64, UnitRange{Int64}})
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:276
[4] (::LinearAlgebra.var"#preprocess_broadcasted##0#preprocess_broadcasted##1"{Diagonal})(x::Adjoint{Int64, UnitRange{Int64}})
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:278
[5] map
@ ./tuple.jl:358 [inlined]
[6] preprocess_broadcasted(::Type{…}, bc::Base.Broadcast.Broadcasted{…})
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:278
[7] copy(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{…}, Tuple{…}, typeof(*), Tuple{…}})
@ LinearAlgebra ~/Dropbox/JuliaPackages/LinearAlgebra.jl/src/structuredbroadcast.jl:291
[8] materialize(bc::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{…}, Nothing, typeof(*), Tuple{…}})
@ Base.Broadcast ./broadcast.jl:893
[9] top-level scope
@ REPL[2]:1
Some type information was truncated. Use `show(err)` to see complete types.This isn't the fault of this PR, so I'll try to address this separately. |
|
Ah sorry didn't mean to close this |
Should fix JuliaLang#1549 (comment). After this, the following works: ```julia julia> a = (1:4)' 1×4 adjoint(::UnitRange{Int64}) with eltype Int64: 1 2 3 4 julia> diagview(a) 1-element view(::UnitRange{Int64}, 1:2:1) with eltype Int64: 1 ```
This PR fixes #1543 (comment).
It overloads
fzerosuch that broadcasts involving structured matrices and zero-absorbing operations (*,dot, etc.) will returntrueforfzeropreserving. Logic was also added for left-absorbing operations (/,rem, etc.) where we attempt to preserve structure if the structured matrix is on the left and as with the scalar case, fall back to a dense matrix if division by zero occurs. There is also a check in the case ofSymTridiagonal, where broadcasts with zero-absorbing operations can avoid materializing a dense matrix, but might need to become aTridiagonal.Unit tests have been added to cover the left-absorbing functionality, and some existing ones have been changed. Specifically,
typeof(broadcast(*, s, fV, fA, X))should now be whateverStructuredMatrixX is, except forUnit***erTriangularwhich becomes***erTriangularandSymTridiagonalwhich becomesTridiagonalforN > 1.