From fb4e4e5bddf5896bf6c8596026ecd275e599b9af Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 6 Aug 2024 19:04:06 +0000 Subject: [PATCH] Don't read destination indices when copying structured matrices (#55322) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- stdlib/LinearAlgebra/src/bidiag.jl | 2 +- stdlib/LinearAlgebra/src/special.jl | 14 +++++++------- stdlib/LinearAlgebra/src/tridiag.jl | 2 +- stdlib/LinearAlgebra/test/special.jl | 17 +++++++++-------- stdlib/LinearAlgebra/test/tridiag.jl | 2 +- 5 files changed, 19 insertions(+), 18 deletions(-) diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl index 0caaec02b5056..24958422015ab 100644 --- a/stdlib/LinearAlgebra/src/bidiag.jl +++ b/stdlib/LinearAlgebra/src/bidiag.jl @@ -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 diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index 9633594574055..28a1b4b1b2eab 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -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) @@ -361,10 +361,10 @@ 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 @@ -372,7 +372,7 @@ 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) diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl index 2ff688f4b4ed1..e217425402df9 100644 --- a/stdlib/LinearAlgebra/src/tridiag.jl +++ b/stdlib/LinearAlgebra/src/tridiag.jl @@ -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 diff --git a/stdlib/LinearAlgebra/test/special.jl b/stdlib/LinearAlgebra/test/special.jl index 2f870373c9586..9bb84ba0e9d03 100644 --- a/stdlib/LinearAlgebra/test/special.jl +++ b/stdlib/LinearAlgebra/test/special.jl @@ -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 @@ -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 @@ -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 @@ -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)) @@ -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 diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl index 5dc1d01e850d8..41a28631b27a0 100644 --- a/stdlib/LinearAlgebra/test/tridiag.jl +++ b/stdlib/LinearAlgebra/test/tridiag.jl @@ -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]]