From 00bb4f2b8591168470be76d8b544ef6812076356 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreas@noack.dk>
Date: Wed, 17 Oct 2018 10:41:58 +0200
Subject: [PATCH] Use a stable inverse for Tridiagonal and SymTridiagonal
 (#29667)

* Use a stable inverse for Tridiagonal and SymTridiagonal

Fixes #29630

* Fix ldlt in the complex case
---
 stdlib/LinearAlgebra/src/ldlt.jl     |  4 +--
 stdlib/LinearAlgebra/src/tridiag.jl  | 39 +++-------------------------
 stdlib/LinearAlgebra/test/tridiag.jl | 20 ++++++++++++++
 3 files changed, 25 insertions(+), 38 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/ldlt.jl b/stdlib/LinearAlgebra/src/ldlt.jl
index 61c282bf524e3..da7fbd3fb723f 100644
--- a/stdlib/LinearAlgebra/src/ldlt.jl
+++ b/stdlib/LinearAlgebra/src/ldlt.jl
@@ -48,13 +48,13 @@ julia> S
   ⋅        0.545455  3.90909
 ```
 """
-function ldlt!(S::SymTridiagonal{T,V}) where {T<:Real,V}
+function ldlt!(S::SymTridiagonal{T,V}) where {T,V}
     n = size(S,1)
     d = S.dv
     e = S.ev
     @inbounds @simd for i = 1:n-1
         e[i] /= d[i]
-        d[i+1] -= abs2(e[i])*d[i]
+        d[i+1] -= e[i]^2*d[i]
     end
     return LDLt{T,SymTridiagonal{T,V}}(S)
 end
diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl
index 483f1a9b2d76d..40a0af469e1e2 100644
--- a/stdlib/LinearAlgebra/src/tridiag.jl
+++ b/stdlib/LinearAlgebra/src/tridiag.jl
@@ -323,44 +323,13 @@ function Base.replace_in_print_matrix(A::SymTridiagonal, i::Integer, j::Integer,
     i==j-1||i==j||i==j+1 ? s : Base.replace_with_centered_mark(s)
 end
 
-#Implements the inverse using the recurrence relation between principal minors
+# Implements the determinant using principal minors
 # a, b, c are assumed to be the subdiagonal, diagonal, and superdiagonal of
 # a tridiagonal matrix.
 #Reference:
 #    R. Usmani, "Inversion of a tridiagonal Jacobi matrix",
 #    Linear Algebra and its Applications 212-213 (1994), pp.413-414
 #    doi:10.1016/0024-3795(94)90414-6
-function inv_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
-    @assert !has_offset_axes(a, b, c)
-    n = length(b)
-    θ = zeros(T, n+1) #principal minors of A
-    θ[1] = 1
-    n>=1 && (θ[2] = b[1])
-    for i=2:n
-        θ[i+1] = b[i]*θ[i]-a[i-1]*c[i-1]*θ[i-1]
-    end
-    φ = zeros(T, n+1)
-    φ[n+1] = 1
-    n>=1 && (φ[n] = b[n])
-    for i=n-1:-1:1
-        φ[i] = b[i]*φ[i+1]-a[i]*c[i]*φ[i+2]
-    end
-    α = Matrix{T}(undef, n, n)
-    for i=1:n, j=1:n
-        sign = (i+j)%2==0 ? (+) : (-)
-        if i<j
-            α[i,j]=(sign)(prod(c[i:j-1]))*θ[i]*φ[j+1]/θ[n+1]
-        elseif i==j
-            α[i,i]=                       θ[i]*φ[i+1]/θ[n+1]
-        else #i>j
-            α[i,j]=(sign)(prod(a[j:i-1]))*θ[j]*φ[i+1]/θ[n+1]
-        end
-    end
-    α
-end
-
-#Implements the determinant using principal minors
-#Inputs and reference are as above for inv_usmani()
 function det_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
     @assert !has_offset_axes(a, b, c)
     n = length(b)
@@ -369,13 +338,12 @@ function det_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}}
         return θa
     end
     θb = b[1]
-    for i=2:n
-        θb, θa = b[i]*θb-a[i-1]*c[i-1]*θa, θb
+    for i in 2:n
+        θb, θa = b[i]*θb - a[i-1]*c[i-1]*θa, θb
     end
     return θb
 end
 
-inv(A::SymTridiagonal) = inv_usmani(A.ev, A.dv, A.ev)
 det(A::SymTridiagonal) = det_usmani(A.ev, A.dv, A.ev)
 
 function getindex(A::SymTridiagonal{T}, i::Integer, j::Integer) where T
@@ -656,7 +624,6 @@ end
 ==(A::Tridiagonal, B::SymTridiagonal) = (A.dl==A.du==B.ev) && (A.d==B.dv)
 ==(A::SymTridiagonal, B::Tridiagonal) = (B.dl==B.du==A.ev) && (B.d==A.dv)
 
-inv(A::Tridiagonal) = inv_usmani(A.dl, A.d, A.du)
 det(A::Tridiagonal) = det_usmani(A.dl, A.d, A.du)
 
 AbstractMatrix{T}(M::Tridiagonal) where {T} = Tridiagonal{T}(M)
diff --git a/stdlib/LinearAlgebra/test/tridiag.jl b/stdlib/LinearAlgebra/test/tridiag.jl
index bef772968778d..b4691c81b2262 100644
--- a/stdlib/LinearAlgebra/test/tridiag.jl
+++ b/stdlib/LinearAlgebra/test/tridiag.jl
@@ -417,4 +417,24 @@ end
           "LinearAlgebra.LU{Float64,LinearAlgebra.Tridiagonal{Float64,SparseArrays.SparseVector")
 end
 
+@testset "Issue 29630" begin
+    function central_difference_discretization(N; dfunc = x -> 12x^2 - 2N^2,
+                                               dufunc = x -> N^2 + 4N*x,
+                                               dlfunc = x -> N^2 - 4N*x,
+                                               bfunc = x -> 114ℯ^-x * (1 + 3x),
+                                               b0 = 0, bf = 57/ℯ,
+                                               x0 = 0, xf = 1)
+        h = 1/N
+        d, du, dl, b = map(dfunc, (x0+h):h:(xf-h)), map(dufunc, (x0+h):h:(xf-2h)),
+                       map(dlfunc, (x0+2h):h:(xf-h)), map(bfunc, (x0+h):h:(xf-h))
+        b[1] -= dlfunc(x0)*b0     # subtract the boundary term
+        b[end] -= dufunc(xf)*bf   # subtract the boundary term
+        Tridiagonal(dl, d, du), b
+    end
+
+    A90, b90 = central_difference_discretization(90)
+
+    @test A90\b90 ≈ inv(A90)*b90
+end
+
 end # module TestTridiagonal