Skip to content

Commit

Permalink
Don't read destination indices when copying structured matrices (Juli…
Browse files Browse the repository at this point in the history 
…aLang#55322)

Fixes the following regression introduced in v1.11
```julia
julia> using LinearAlgebra

julia> D = Diagonal(rand(4));

julia> T = Tridiagonal(Vector{BigFloat}(undef, 3), Vector{BigFloat}(undef, 4), Vector{BigFloat}(undef, 3))
4×4 Tridiagonal{BigFloat, Vector{BigFloat}}:
 #undef  #undef     ⋅       ⋅ 
 #undef  #undef  #undef     ⋅ 
    ⋅    #undef  #undef  #undef
    ⋅       ⋅    #undef  #undef

julia> copyto!(T, D)
ERROR: UndefRefError: access to undefined reference
Stacktrace:
  [1] getindex
    @ ./essentials.jl:907 [inlined]
  [2] _broadcast_getindex
    @ ./broadcast.jl:644 [inlined]
  [3] _getindex
    @ ./broadcast.jl:675 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:650 [inlined]
  [5] getindex
    @ ./broadcast.jl:610 [inlined]
  [6] macro expansion
    @ ./broadcast.jl:973 [inlined]
  [7] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [8] copyto!
    @ ./broadcast.jl:972 [inlined]
  [9] copyto!
    @ ./broadcast.jl:925 [inlined]
 [10] materialize!
    @ ./broadcast.jl:883 [inlined]
 [11] materialize!
    @ ./broadcast.jl:880 [inlined]
 [12] _copyto_banded!(T::Tridiagonal{BigFloat, Vector{BigFloat}}, D::Diagonal{Float64, Vector{Float64}})
    @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/special.jl:323
 [13] copyto!(dest::Tridiagonal{BigFloat, Vector{BigFloat}}, src::Diagonal{Float64, Vector{Float64}})
    @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/special.jl:315
 [14] top-level scope
    @ REPL[4]:1
```
After this PR
```julia
julia> copyto!(T, D)
4×4 Tridiagonal{BigFloat, Vector{BigFloat}}:
 0.909968  0.0        ⋅         ⋅ 
 0.0       0.193341  0.0        ⋅ 
  ⋅        0.0       0.194794  0.0
  ⋅         ⋅        0.0       0.506905
```
The current implementation used an optimization that may not be
applicable for non-isbits types, and this PR ensures that we always read
from the source and write to the destination.
  • Loading branch information
@jishnub
jishnub authored Aug 6, 2024
1 parent 09e5c40 commit fb4e4e5
Show file tree
Hide file tree
Showing 5 changed files with 19 additions and 18 deletions.
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/bidiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -305,7 +305,7 @@ function _copyto_banded!(A::Bidiagonal, B::Bidiagonal)
if A.uplo == B.uplo
A.ev .= B.ev
elseif iszero(B.ev) # diagonal source
A.ev .= zero.(A.ev)
A.ev .= B.ev
else
zeroband = istriu(A) ? "lower" : "upper"
uplo = A.uplo
Expand Down
14 changes: 7 additions & 7 deletions stdlib/LinearAlgebra/src/special.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,20 +320,20 @@ function copyto!(dest::BandedMatrix, src::BandedMatrix)
end
function _copyto_banded!(T::Tridiagonal, D::Diagonal)
T.d .= D.diag
T.dl .= zero.(T.dl)
T.du .= zero.(T.du)
T.dl .= view(D, diagind(D, -1, IndexStyle(D)))
T.du .= view(D, diagind(D, 1, IndexStyle(D)))
return T
end
function _copyto_banded!(SymT::SymTridiagonal, D::Diagonal)
issymmetric(D) || throw(ArgumentError("cannot copy a non-symmetric Diagonal matrix to a SymTridiagonal"))
SymT.dv .= D.diag
_ev = _evview(SymT)
_ev .= zero.(_ev)
_ev .= view(D, diagind(D, 1, IndexStyle(D)))
return SymT
end
function _copyto_banded!(B::Bidiagonal, D::Diagonal)
B.dv .= D.diag
B.ev .= zero.(B.ev)
B.ev .= view(D, diagind(D, B.uplo == 'U' ? 1 : -1, IndexStyle(D)))
return B
end
function _copyto_banded!(D::Diagonal, B::Bidiagonal)
Expand Down Expand Up @@ -361,18 +361,18 @@ function _copyto_banded!(T::Tridiagonal, B::Bidiagonal)
T.d .= B.dv
if B.uplo == 'U'
T.du .= B.ev
T.dl .= zero.(T.dl)
T.dl .= view(B, diagind(B, -1, IndexStyle(B)))
else
T.dl .= B.ev
T.du .= zero.(T.du)
T.du .= view(B, diagind(B, 1, IndexStyle(B)))
end
return T
end
function _copyto_banded!(SymT::SymTridiagonal, B::Bidiagonal)
issymmetric(B) || throw(ArgumentError("cannot copy a non-symmetric Bidiagonal matrix to a SymTridiagonal"))
SymT.dv .= B.dv
_ev = _evview(SymT)
_ev .= zero.(_ev)
_ev .= B.ev
return SymT
end
function _copyto_banded!(B::Bidiagonal, T::Tridiagonal)
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/src/tridiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1040,7 +1040,7 @@ function _copyto_banded!(A::Tridiagonal, B::SymTridiagonal)
return A
end
function _copyto_banded!(A::SymTridiagonal, B::Tridiagonal)
issymmetric(B) || throw(ArgumentError("cannot copy a non-symmetric Tridiagonal matrix to a SymTridiagonal"))
issymmetric(B) || throw(ArgumentError("cannot copy an asymmetric Tridiagonal matrix to a SymTridiagonal"))
A.dv .= B.d
_evview(A) .= B.du
return A
Expand Down
17 changes: 9 additions & 8 deletions stdlib/LinearAlgebra/test/special.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -555,8 +555,8 @@ end
@testset "from Diagonal" begin
D = Diagonal(d)
@testset "to Bidiagonal" begin
BU = Bidiagonal(zero(d), oneunit.(du), :U)
BL = Bidiagonal(zero(d), oneunit.(dl), :L)
BU = Bidiagonal(similar(d, BigInt), similar(du, BigInt), :U)
BL = Bidiagonal(similar(d, BigInt), similar(dl, BigInt), :L)
for B in (BL, BU)
copyto!(B, D)
@test B == D
Expand All @@ -573,7 +573,7 @@ end
end
end
@testset "to Tridiagonal" begin
T = Tridiagonal(oneunit.(dl), zero(d), oneunit.(du))
T = Tridiagonal(similar(dl, BigInt), similar(d, BigInt), similar(du, BigInt))
copyto!(T, D)
@test T == D

Expand All @@ -586,8 +586,8 @@ end
end
end
@testset "to SymTridiagonal" begin
for du2 in (oneunit.(du), oneunit.(d))
S = SymTridiagonal(zero(d), du2)
for du2 in (similar(du, BigInt), similar(d, BigInt))
S = SymTridiagonal(similar(d), du2)
copyto!(S, D)
@test S == D
end
Expand Down Expand Up @@ -630,13 +630,14 @@ end
end
end
@testset "to Tridiagonal" begin
T = Tridiagonal(oneunit.(dl), zero(d), oneunit.(du))
T = Tridiagonal(similar(dl, BigInt), similar(d, BigInt), similar(du, BigInt))
for B in (BL, BU, BLones, BUones)
copyto!(T, B)
@test T == B
end

@testset "mismatched size" begin
T = Tridiagonal(oneunit.(dl), zero(d), oneunit.(du))
for uplo in (:L, :U)
T .= 0
copyto!(T, Bidiagonal([1], Int[], uplo))
Expand All @@ -647,8 +648,8 @@ end
end
end
@testset "to SymTridiagonal" begin
for du2 in (oneunit.(du), oneunit.(d))
S = SymTridiagonal(zero(d), du2)
for du2 in (similar(du, BigInt), similar(d, BigInt))
S = SymTridiagonal(similar(d, BigInt), du2)
for B in (BL, BU)
copyto!(S, B)
@test S == B
Expand Down
2 changes: 1 addition & 1 deletion stdlib/LinearAlgebra/test/tridiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -808,7 +808,7 @@ end
@test copyto!(zero(S), T) == T

T2 = Tridiagonal(ones(length(ev)), zero(dv), zero(ev))
@test_throws "cannot copy a non-symmetric Tridiagonal matrix to a SymTridiagonal" copyto!(zero(S), T2)
@test_throws "cannot copy an asymmetric Tridiagonal matrix to a SymTridiagonal" copyto!(zero(S), T2)

@testset "mismatched sizes" begin
dv2 = [4; @view dv[2:end]]
Expand Down

0 comments on commit fb4e4e5

Please sign in to comment.