LinearAlgebra - inv() producing NaNs when called on a Tridiagonal

[MichaelDS]

I’m creating a 90x90 tridiagonal matrix using the function below, which I store in the variable A90. When I call inv(A90), I get a matrix of all NaNs, but if I call inv(Matrix(A90)), the inverse is successfully computed. Is this behavior expected?

using LinearAlgebra
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)

typeof(A90)  # Tridiagonal{Float64,Array{Float64,1}}

inv(A90)  # Produces a matrix of all NaNs
inv(Matrix(A90))  # Successfully computes the inverse

A80, b80 = central_difference_discretization(80)
inv(A80)  # Works fine on Tridiagonals smaller than 82x82

@tpapp suggested that it might be as issue related to conditioning and the Tridiagonal inversion method, which I suspect as well because my code generates an ill-conditioned matrix by design.

[RalphAS]

This issue arises because the special-case method for inverting tridiagonals (inv_usmani) has no scaling protection. Using the factorization instead

lu(A90) \ Matrix(1.0I,size(A90))

gives a good result via a stable banded scheme (faster than working with a full matrix).

Furthermore, inv_usmani is currently type-unstable, so it takes much longer.

[andreasnoack]

@RalphAS Thanks for identifying the problem. I think we should just remove inv_usmani and instead use the fallback which solves against the identity.