Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LinearAlgebra: copyto! between banded matrix types #54041

Merged
merged 2 commits into from
Apr 23, 2024

Conversation

jishnub
Copy link
Contributor

@jishnub jishnub commented Apr 11, 2024

This specialized copyto! for combinations of banded structured matrix types so that the copy may be O(N) instead of the fallback O(N^2) implementation.

E.g.:

julia> T = Tridiagonal(zeros(999), zeros(1000), zeros(999));

julia> B = Bidiagonal(ones(1000), fill(2.0, 999), :U);

julia> @btime copyto!($T, $B);
  1.927 ms (0 allocations: 0 bytes) # master
  229.870 ns (0 allocations: 0 bytes) # PR

This also changes the copyto! implementation for mismatched matrix sizes, bringing it closer to the docstring. So, the following works on master:

julia> Ddest = Diagonal(zeros(4));

julia> Dsrc = Diagonal(ones(2));

julia> copyto!(Ddest, Dsrc)
4×4 Diagonal{Float64, Vector{Float64}}:
 1.0            
     1.0        
         0.0    
             0.0

but this won't work anymore with this PR. This was inconsistent anyway, as materializing the matrices produces a different result, which shouldn't be the case:

julia> copyto!(Matrix(Ddest), Dsrc)
4×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0

After this PR, the way to carry out the copy would be

julia> copyto!(Ddest, CartesianIndices(Dsrc), Dsrc, CartesianIndices(Dsrc))
4×4 Diagonal{Float64, Vector{Float64}}:
 1.0            
     1.0        
         0.0    
             0.0

This change fixes JuliaLang/LinearAlgebra.jl#932.

Also fixes JuliaLang/LinearAlgebra.jl#1061
After this,

julia> @btime copyto!(C, B) setup=(n = 1_000; B = Bidiagonal(randn(n), randn(n-1), :L); C = Bidiagonal(randn(n), randn(n-1), :L));
  158.405 ns (0 allocations: 0 bytes)

julia> @btime copyto!(C, B) setup=(n = 10_000; B = Bidiagonal(randn(n), randn(n-1), :L); C = Bidiagonal(randn(n), randn(n-1), :L));
  4.706 μs (0 allocations: 0 bytes)

julia> @btime copyto!(C, B) setup=(n = 100_000; B = Bidiagonal(randn(n), randn(n-1), :L); C = Bidiagonal(randn(n), randn(n-1), :L));
  120.880 μs (0 allocations: 0 bytes)

which is roughly linear scaling.

Taken along with #54027, the speed-ups would also apply to the adjoints of banded matrices.

@jishnub jishnub added linear algebra Linear algebra arrays [a, r, r, a, y, s] labels Apr 11, 2024
Copy link
Member

@dkarrasch dkarrasch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't checked every single line, but the design LGTM. Also, _copyto_banded! seems to cover all 16 possible cases.

@jishnub jishnub merged commit 2ca9e00 into master Apr 23, 2024
9 of 11 checks passed
@jishnub jishnub deleted the jishnub/copytostructured branch April 23, 2024 05:10
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
arrays [a, r, r, a, y, s] linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

special case copyto!(::Bidiagonal, ::Bidiagonal) Inconsistency in copyto! for structured matrices
2 participants