Dear All,

I am trying to solve linear systems `Ax = rhs`

where `A`

has a tridiagonal shape. When the type of `A`

is tridiagonal, one can use the thomas algo to solve it. I find it faster than defining `A`

as a `Tridiagonal`

and calling `\`

. Indeed, the actual implementation of `\`

uses a LU decomposition followed by a solve. Maybe the choice of not relying on `\`

for `Tridiagonal`

is numerical stability of the algorithm.

I want to do something slighly different. I want to solve the same problem for a banded matrix with a diagonal band and two other bands. When the bandwith is 1, one recovers the case above of the Thomas algo.

I have the MWE below, but it does not work. As you can see, I first compute the LU decomposition of the matrix and then solve the linear systems. I have two issues I cannot solve

- Firstly despite my analytical formula, my LU decomposition seems wrong.
- Secondly, I am able to invert L but I am unable to invert U which I am ashamed of.

If anyone has an idea, Iād be delighted.

Thank you a lot.

```
using LinearAlgebra, SparseArrays
function thomas2bands!(a::vec, b::vec, c::vec, rhs::vec) where {T, vec <: Vector{T}}
@assert length(a) == length(c)
n = length(b)
p = length(a)
# we first derive the LU decomposition
# we have L and U written 8 lines below, L = [l, 1] and U = [v, c]
v = copy(b)
l = a ./ v[1:p]
for ii = n-p+1:n
v[ii] = b[ii] - l[ii-n+p] * c[ii-n+p]
end
# up to here, this has been check, the LU decomposition seems OK
@show n-p
L = spdiagm(-n+p => l, 0 => fill(1.0, n))
U = spdiagm(0 => v, n-p => c)
A = spdiagm(-n+p => a, 0 => b, n-p => c)
println("-> A - LU = ", norm(L * U - A, Inf64))
# we want to solve LUx = rhs
# we first solve Ly = rhs
y = similar(rhs)
for ii=1:n-p
y[ii] = rhs[ii]
end
for ii=n-p+1:n
y[ii] = rhs[ii] - l[ii-n+p] * y[ii-n+p]
end
println("-> L \\ rhs = ",norm(y - (L \ rhs), Inf64))
# we solve Ux = y
x = similar(rhs)
for ii = n-p+1:n
x[ii] = y[ii] / v[ii]
end
for ii=1:n-p
x[ii] = y[ii] / v[ii] - x[ii+n-p] * c[ii] / v[ii]
end
println("-> y - U * x = ",norm((y - (U * x)), Inf64))
println("-> U \\ y = ",norm(x - (U \ y),Inf64))
return x
end
sol = @time thomas2bands!(rand(10^5-10^3), 2 .+ rand(10^5), rand(10^5-10^3), rand(10^5) )
```