diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 4afb855117c14..25e32a8ab0828 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -520,9 +520,6 @@ end # Initialize variables info = 0 - if dl === du - throw(ArgumentError("off-diagonals must not alias")) - end fill!(du2, 0) @inbounds begin diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl index d71519cfa5365..235fa5636877a 100644 --- a/stdlib/LinearAlgebra/src/tridiag.jl +++ b/stdlib/LinearAlgebra/src/tridiag.jl @@ -472,13 +472,13 @@ struct Tridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T} "lengths of subdiagonal, diagonal and superdiagonal: ", "($(length(dl)), $(length(d)), $(length(du)))"))) end - new{T,V}(dl, d, du) + new{T,V}(dl, d, Base.unalias(dl, du)) end # constructor used in lu! function Tridiagonal{T,V}(dl, d, du, du2) where {T,V<:AbstractVector{T}} require_one_based_indexing(dl, d, du, du2) # length checks? - new{T,V}(dl, d, du, du2) + new{T,V}(dl, d, Base.unalias(dl, du), du2) end end @@ -491,6 +491,10 @@ solvers, but may be converted into a regular matrix with [`convert(Array, _)`](@ref) (or `Array(_)` for short). The lengths of `dl` and `du` must be one less than the length of `d`. +!!! note + The subdiagonal `dl` and the superdiagonal `du` must not be aliased to each other. + If aliasing is detected, the constructor will use a copy of `du` as its argument. + # Examples ```jldoctest julia> dl = [1, 2, 3]; @@ -912,9 +916,6 @@ function ldiv!(A::Tridiagonal, B::AbstractVecOrMat) dl = A.dl d = A.d du = A.du - if dl === du - throw(ArgumentError("off-diagonals of `A` must not alias")) - end @inbounds begin for i in 1:n-1 diff --git a/stdlib/LinearAlgebra/test/lu.jl b/stdlib/LinearAlgebra/test/lu.jl index 2facbc55c8aff..d96510549caa5 100644 --- a/stdlib/LinearAlgebra/test/lu.jl +++ b/stdlib/LinearAlgebra/test/lu.jl @@ -97,7 +97,9 @@ dimg = randn(n)/2 dlu = convert.(eltya, [1, 1]) dia = convert.(eltya, [-2, -2, -2]) tri = Tridiagonal(dlu, dia, dlu) - @test_throws ArgumentError lu!(tri) + L = lu(tri) + @test lu!(tri) == L + @test UpperTriangular(tri) == L.U end end @testset for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl index 68a792b55d396..2d846683e38c3 100644 --- a/stdlib/LinearAlgebra/test/tridiag.jl +++ b/stdlib/LinearAlgebra/test/tridiag.jl @@ -82,7 +82,7 @@ end @test TT == Matrix(TT) @test TT.dl === y @test TT.d === x - @test TT.du === y + @test TT.du == y @test typeof(TT)(TT) === TT end ST = SymTridiagonal{elty}([1,2,3,4], [1,2,3]) @@ -521,6 +521,14 @@ end @test Tridiagonal(4:5, 1:3, 1:2) == [1 1 0; 4 2 2; 0 5 3] end +@testset "Prevent off-diagonal aliasing in Tridiagonal" begin + e = ones(4) + f = e[1:end-1] + T = Tridiagonal(f, 2e, f) + T ./= 10 + @test all(==(0.1), f) +end + @testset "Issue #26994 (and the empty case)" begin T = SymTridiagonal([1.0],[3.0]) x = ones(1)