Unalias off-diagonals in Tridiagonal constructor (#51763)

Author: jishnub

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

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

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)

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)