Use a stable inverse for Tridiagonal and SymTridiagonal (JuliaLang#29667

)

* Use a stable inverse for Tridiagonal and SymTridiagonal

Fixes JuliaLang#29630

* Fix ldlt in the complex case

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

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)

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