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

Inconsistency in copyto! for structured matrices #932

Closed
jishnub opened this issue Jul 12, 2022 · 1 comment · Fixed by JuliaLang/julia#54041
Closed

Inconsistency in copyto! for structured matrices #932

jishnub opened this issue Jul 12, 2022 · 1 comment · Fixed by JuliaLang/julia#54041
Labels
arrays [a, r, r, a, y, s]

Comments

@jishnub
Copy link
Collaborator

jishnub commented Jul 12, 2022

MWE

julia> D1 = Diagonal(ones(4))
4×4 Diagonal{Float64, Vector{Float64}}:
 1.0            
     1.0        
         1.0    
             1.0

julia> D2 = Diagonal(2ones(3))
3×3 Diagonal{Float64, Vector{Float64}}:
 2.0        
     2.0    
         2.0

julia> copyto!(Matrix(D1), Matrix(D2))
4×4 Matrix{Float64}:
 2.0  2.0  2.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

julia> copyto!(D1, D2)
4×4 Diagonal{Float64, Vector{Float64}}:
 2.0            
     2.0        
         2.0    
             1.0

The docstring of copyto! states

  copyto!(dest::AbstractArray, src) -> dest

  Copy all elements from collection src to array dest, whose length must be greater than or equal to the length n of src. The first n elements of dest are overwritten, the other
  elements are left untouched.

This linear copying doesn't seem to be followed by structured matrices.

Version info:

julia> versioninfo()
Julia Version 1.7.3
Commit 742b9abb4d (2022-05-06 12:58 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, tigerlake)
Environment:
  JULIA_EDITOR = vim

but the same behavior is there on

Julia Version 1.9.0-DEV.888
Commit eb72c2aaf2 (2022-07-05 12:07 UTC)
@N5N3
Copy link
Member

N5N3 commented Jul 12, 2022

It seems that we should perform an axes check before trying the structural copyto!

List to be fixed:

copyto!(A::T, B::T) where T<:Union{UnitUpperTriangular, UpperTriangular}  # miss unaliasing
copyto!(A::T, B::T) where T<:Union{LowerTriangular, UnitLowerTriangular}  # miss unaliasing
copyto!(D1::Diagonal, D2::Diagonal)
copyto!(dest::Tridiagonal, src::Tridiagonal)
copyto!(dest::SymTridiagonal, src::SymTridiagonal)
copyto!(dest::Symmetric, src::Symmetric)   # miss unaliasing in `transpose!` branch
copyto!(dest::Hermitian, src::Hermitian)   # miss unaliasing in `adjoint!` branch

@N5N3 N5N3 added the arrays [a, r, r, a, y, s] label Jul 12, 2022
N5N3 referenced this issue in N5N3/julia Jul 25, 2022
N5N3 referenced this issue in N5N3/julia Aug 16, 2022
N5N3 referenced this issue in N5N3/julia Sep 13, 2022
N5N3 referenced this issue in N5N3/julia Oct 29, 2022
N5N3 referenced this issue in N5N3/julia Nov 30, 2022
N5N3 referenced this issue in N5N3/julia Nov 30, 2022
N5N3 referenced this issue in N5N3/julia Dec 8, 2022
N5N3 referenced this issue in N5N3/julia Dec 10, 2022
N5N3 referenced this issue in N5N3/julia Jan 6, 2023
N5N3 referenced this issue in N5N3/julia Jan 12, 2023
jishnub referenced this issue in JuliaLang/julia Apr 23, 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
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
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
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
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 https://github.com/JuliaLang/julia/issues/46005.

Also fixes https://github.com/JuliaLang/julia/issues/53997
After this,
```julia
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.
KristofferC referenced this issue Nov 14, 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
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
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
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
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 https://github.com/JuliaLang/julia/issues/46005.

Also fixes https://github.com/JuliaLang/julia/issues/53997
After this,
```julia
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 JuliaLang/julia#54027, the
speed-ups would also apply to the adjoints of banded matrices.
@KristofferC KristofferC transferred this issue from JuliaLang/julia Nov 26, 2024
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]
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants