# How to improve the Thomas algorithm for block tridiagonal matrices

**URL:** <https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894>\
**Category:** Numerics\
**Created:** [March 13, 2025, 4:00am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894 "2025-03-13T04:00:20Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [March 13, 2025, 4:00am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/1 "2025-03-13T04:00:20Z")

</div>

I did some searching but it seems that Julia doesn’t implement a specific algorithm for block tridiagonal matrices. So I implemented one myself. I would like to know:

1. Have I missed some existing algorithm implementations for block tridiagonal matrices?
2. is there any room for improvement of the solver function I wrote?

```julia
function block_thomas_tridiagonal(A::AbstractVector{T}, B::AbstractVector{T}, C::AbstractVector{T},
    d::AbstractVector) where {T}
    B = copy(B)
    n = size(B[1], 1) # rows of each block
    m = length(B) # number of blocks

    # Forward elimination
    for i = 2:m # (no speed up using inbounds+simd)
        multiplier = A[i-1] * inv(B[i-1])
        B[i] -= multiplier * C[i-1]
        d[n*(i-1)+1:n*i] .-= multiplier * d[n*(i-2)+1:n*(i-1)]
    end

    # Backward substitution
    x = similar(d)
    x[n*(m-1)+1:n*m] = B[m] \ d[n*(m-1)+1:n*m] # no speed up using .=
    for i = m-1:-1:1
        x[n*(i-1)+1:n*i] = B[i] \ (d[n*(i-1)+1:n*i] - C[i] * x[n*i+1:n*(i+1)])
    end

    return x
end

```

This function is specifically designed to handle the case where every block has the same size. It usually takes a Vector{Matrix{Float64}} as input and solves each block one by one.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [March 13, 2025, 8:37am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/2 "2025-03-13T08:37:46Z")

</div>

This parallel (thread+simd) C++ implementation may help you to implement a performant Julia implementation of Thomas Algorithm. Be aware that for large problem exceeding cache size, the algorihm is memory bound:

> **[GitHub - triscale-innov/Legolas](https://github.com/triscale-innov/Legolas)**
>
> Contribute to triscale-innov/Legolas development by creating an account on GitHub.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/1/e14be852310f8d4ee7f12b45fc002516c49d8ddd.jpeg)

Note also that this algorithm is well suited for GPU implementation.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 13, 2025, 11:18am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/3 "2025-03-13T11:18:29Z")

</div>

> [@karei](#):
>
> `multiplier = A[i-1] * inv(B[i-1])`

Don’t use explicit matrix inversion. You could use `A[i-1] / inv(B[i-1])`. If the blocks are all the same size, as you seem to be assuming here, you could preallocate buffers for the LU factorization, the multipliers, and the other matrix multiplications.

> [@karei](#):
>
> `# no speed up using .=`

That’s because `.=` is for elementwise/broadcasting operations, not matrix operations. If you want to do the latter in-place, use `ldiv!` and `mul!` etc.

Beware that this algorithm may be numerically unstable, IIRC, e.g. if any of the blocks are nearly singular.

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [March 13, 2025, 12:46pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/4 "2025-03-13T12:46:12Z")

</div>

> [@stevengj](#):
>
> Beware that this algorithm may be numerically unstable, IIRC, e.g. if any of the blocks are nearly singular.

It is unfortunately. It’s even unstable for a non-block symmetric tridiagonal. (And the block case is going be worse with the explicit inverse you pointed out.) You really need some form of pivoting for this to guarantee stability. I put an example for the factorization variant used in `ldlt` in a [post](https://discourse.julialang.org/t/solving-shifted-linear-system/111265/10) last year.

I just showed the factorization backward error, but the thing that seems more troubling to me is that the symmetric tridiagonal solver in Julia uses an unstable algorithm. With a relative residual for the solver, the example is:

```julia
using LinearAlgebra
A = SymTridiagonal([
  3e-10 1 0
  1 2 1 
  0 1 2
])

F = ldlt(A)
display(A - F.L * F.D * F.L')

b=[1.0, 2.0, 3.0]
x=A\b
@show norm(A*x-b) / opnorm(A) / norm(x)

```

which gives output

```julia
3×3 Matrix{Float64}:
 0.0 0.0 0.0
 0.0 4.76837e-7 0.0
 0.0 0.0 0.0
(norm(A * x - b) / opnorm(A)) / norm(x) = 8.692567662779919e-8
8.692567662779919e-8

```

LAPACK avoids this by using LU with partial pivoting for its tridiagonal solver. However, in addition to the permutation, you pick up an extra superdiagonal in U. So it isn’t a fully in place algorithm like the elegant algorithm implemented in `ldlt!`. The alternative is Bunch’s method that uses 2\times 2 diagonal blocks mentioned in my other post. You don’t need any permutations, but you do need to accept 2\times 2 blocks in D. Bunch’s algorithm has the further advantage of preserving symmetry, so by finding the eigenvalues in D you can reliably compute the inertia of a matrix. However, it also requires extra storage and you can’t represent the triangular factor as a `UnitLowerTriangular` a `SymTriadiagonal` as is done now. The structure of the factors doesn’t align as nicely with existing types in `LinearAlgebra`, which is annoying even if you are fine giving up on a fully in-place factorization.

---

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [March 16, 2025, 6:03am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/5 "2025-03-16T06:03:20Z")

</div>

> [@LaurentPlagne](#):
>
> [GitHub - triscale-innov/Legolas](https://github.com/triscale-innov/Legolas)

Thanks for the link. But I took a closer look in there and there is no implementation of Thomas’ algorithm for solving block tridiagonal matrices. There is only the normal Thomas algorithm. But it did help me improve the Thomas algorithm I originally wrote, thank you.

> [@stevengj](#):
>
> use `ldiv!` and `mul!` etc.

Thank you for the reply. I spent a little time on how to improve this algorithm based on your suggestion. Here is a comparison between the original and the improved version

original

```julia
function block_thomas_tridiagonal(L::AbstractVector{T}, D::AbstractVector{T}, U::AbstractVector{T},
    b::AbstractVector) where {T}
    D = deepcopy(D)
    n = size(D[1], 1) # each block has size n*n
    m = length(D) # num of blocks

    # Forward elimination
    for i = 1:m-1
        multiplier = L[i] * inv(D[i])
        D[i+1] -= multiplier * U[i]
        b[n*i+1:n*(i+1)] .-= multiplier * b[n*(i-1)+1:n*i]
    end

    # Backward substitution
    x = similar(b)
    x[n*(m-1)+1:n*m] = D[m] \ b[n*(m-1)+1:n*m] 
    for i = m-1:-1:1
        x[n*(i-1)+1:n*i] = D[i] \ (b[n*(i-1)+1:n*i] - U[i] * x[n*i+1:n*(i+1)])
    end

    return x
end

```

and improved

```julia
function block_thomas_tridiagonal!(x::AbstractVector, L::AbstractVector{T}, D::AbstractVector{T},
    U::AbstractVector{T}, b::AbstractVector) where {T}
    n = size(D[1], 1) # each block has size n*n
    m = length(D) # num of blocks
    (; lu_buf, D_buf) = get_block_thomas_buffer(n, m)

    # Forward elimination
    for i = 1:m-1
        lu_buf[i] = lu!(D[i])
        rdiv!(L[i], lu_buf[i])
        D[i+1] .-= mul!(D_buf[i], L[i], U[i])
        b[n*i+1:n*(i+1)] .-= L[i] * b[n*(i-1)+1:n*i]
    end

    # Backward substitution
    ldiv!(@view(x[n*(m-1)+1:n*m]), lu!(D[m]), b[n*(m-1)+1:n*m])
    for i = m-1:-1:1
        ldiv!(@view(x[n*(i-1)+1:n*i]), lu_buf[i], (b[n*(i-1)+1:n*i] - U[i] * x[n*i+1:n*(i+1)]))
    end

    return x
end

```

where the buffer is a vector of structures I wrote myself

```julia
mutable struct BlockThomas
    n::Int
    m::Int
    lu_buf::Vector{LU{Float64,Matrix{Float64},Vector{Int64}}}
    D_buf::Vector{Matrix{Float64}}
    function BlockThomas(n, m)
        lu_buf = Vector{LU{Float64,Matrix{Float64},Vector{Int64}}}(undef, m)
        D_buf = [mat(n, n) for _ = 1:m]
        return new(n, m, lu_buf, D_buf)
    end
end

const BlockThomasBuf = Vector{BlockThomas}(undef, nthreads())

function get_block_thomas_buffer(n::Int, m::Int)
    tid = threadid()
    return if !isassigned(BlockThomasBuf, tid)
        BlockThomasBuf[tid] = BlockThomas(n, m)
    elseif BlockThomasBuf[tid].n != n || BlockThomasBuf[tid].m != m
        BlockThomasBuf[tid] = BlockThomas(n, m)
    else
        BlockThomasBuf[tid]
    end
end

```

The average time for the benchmark changed from `460.868 μs Memory estimate: 692.58 KiB, allocs estimate: 1893` to `325.692 μs Memory estimate: 114.73 KiB, allocs estimate: 893.`, a 29% savings on this medium-sized example.

For the benefit of anyone who sees this post later, I’ve put both pieces of code in [my repository](https://github.com/abcdvvvv/Thomas-Tridiagonal-Algorithm).

The improved version does not explicitly use matrix inversion, but the two seem to be about the same in accuracy.

```julia
Random.seed!(1)
n = 10
m = 100
A = [rand(BigFloat, n, n) for _ = 1:m-1]
B = [rand(BigFloat, n, n) for _ = 1:m]
C = [rand(BigFloat, n, n) for _ = 1:m-1]
d = rand(BigFloat, m * n)
x = similar(d);
A_float64 = [Float64.(i) for i in A]
B_float64 = [Float64.(i) for i in B]
C_float64 = [Float64.(i) for i in C]
d_float64 = Float64.(d)
x_float64 = similar(d_float64)
original_result = block_thomas_tridiagonal(A, B, C, d)
result_float64 = block_thomas_tridiagonal(A_float64, B_float64, C_float64, d_float64)
error_norm = norm(original_result - result_float64, 1)
println("1-norm of the error: ", error_norm)

# original_result = block_thomas_tridiagonal!(x, A, B, C, d)
# result_float64 = block_thomas_tridiagonal!(x_float64, A_float64, B_float64, C_float64, d_float64)
# error_norm = norm(original_result - result_float64, 1)
# println("1-norm of the error: ", error_norm)

```

original: `1-norm of the error: 5.06330e-09`  
improved: `1-norm of the error: 1.1942802e-08`  
If you change the seed, sometimes the improved version has a smaller error.

But it does demonstrate that the Block Thomas algorithm has accumulated errors of the order of `1e-9` on this medium-sized example.

> [@stevengj](#):
>
> Beware that this algorithm may be numerically unstable, IIRC, e.g. if any of the blocks are nearly singular.

May I ask if there is a better algorithm for solving block tridiagonal matrix? I would like to trade-off accuracy and performance and utilize sparsity instead of solving the entire matrix.

---

<div class="post-metadata">

**Author:** ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)\
**Post date:** [March 16, 2025, 7:45am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/6 "2025-03-16T07:45:29Z")

</div>

If your block are sparse, you should just build the whole sparse matrix and use either KLU or Kyrlov. If your blocks are Dense, Thomas will be the best (even though not stable as said before).

---

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [March 16, 2025, 8:01am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/7 "2025-03-16T08:01:56Z")

</div>

That’s a really clear answer and solves my doubts, thank you!

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 16, 2025, 8:15am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/8 "2025-03-16T08:15:02Z")

</div>

Yes but if not stable why bother?

---

<div class="post-metadata">

**Author:** ![yolhan\_mannes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yolhan_mannes/32/220485_2.png) [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)\
**Post date:** [March 16, 2025, 8:19am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/9 "2025-03-16T08:19:11Z")

</div>

It depends but for discontinuous galerkin problem (block tridiag) you don’t need that much stability of the linear problem as an example. In this case, if you use monomial bases you get dense blocks and if you use lobatto-legendre you get diagonal matricies in block and I think it’s better to just make everything sparse in this case

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [March 16, 2025, 8:21am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/10 "2025-03-16T08:21:10Z")

</div>

Or maybe as a preconditioner

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 16, 2025, 12:02pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/11 "2025-03-16T12:02:24Z")

</div>

> [@karei](#):
>
> `b[n*i+1:n*(i+1)] .-= L[i] * b[n*(i-1)+1:n*i]`

Using `*` instead of `.*` means that the multiplication is not fused with the `.-=` — it allocates a separate array, in a separate loop. See the [“more dots” performance tip](https://docs.julialang.org/en/v1/manual/performance-tips/#More-dots:-Fuse-vectorized-operations). You might want to read [this article](https://julialang.org/blog/2017/01/moredots/) about Julia’s “dot” notation and what it does.

Also, `b[n*(i-1)+1:n*i]` allocates a new array, and also `b[n*i+1:n*(i+1)]` on the right-hand side, since you aren’t using `@views`. See the [“consider using views”](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-views) performance tip. Probably best to just put `@views` in front of `function` so that you use it everywhere in the function, since right now you are allocating _lots_ of copies for slices. (Then you can get rid of the `@view` calls.)

> [@karei](#):
>
> `lu_buf[i] = lu!(D[i])`

I thought you didn’t want to overwrite the input array?

> [@karei](#):
>
> `D[i+1] .-= mul!(D_buf[i], L[i], U[i])`

Why are you using a separate pre-allocated buffer `D_buf[i]` for every loop iteration, rather than just allocating a single buffer and re-using it?

> [@karei](#):
>
> `b[n*i+1:n*(i+1)] .-= L[i] * b[n*(i-1)+1:n*i]`

Not using `mul!` with a buffer? Similarly for the `U[i] *` operation later?

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [March 16, 2025, 12:14pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/12 "2025-03-16T12:14:05Z")

</div>

> [@karei](#):
>
> Thanks for the link. But I took a closer look in there and there is no implementation of Thomas’ algorithm for solving block tridiagonal matrices. There is only the normal Thomas algorithm. But it did help me improve the Thomas algorithm I originally wrote, thank you.

Actually there is but it is a bit tricky to understand. The “normal” Thomas Algorithm (TA) implementation is automatically transformed into a TA applied on fixed size pack of numbers where usual operations (+,-,\*,/…) are transformed (at compile time) into the corresponding SIMD operations. This generated SIMD version of TA is then applied in parallel (threads) to a collection of different data (pmap operation).  
More details here

- [Legolas/presentation.pdf at master · triscale-innov/Legolas · GitHub](https://github.com/triscale-innov/Legolas/blob/master/presentation.pdf) (slides)
- [https://www.researchgate.net/publication/317485219\_Portable\_vectorization\_and\_parallelization\_of\_C\_multi-dimensional\_array\_computations](https://www.researchgate.net/publication/317485219_Portable_vectorization_and_parallelization_of_C_multi-dimensional_array_computations) (paper)

The main idea here is to understand that the TA is intrinsically sequential but that it can be applied in parallel to different tridiagonal blocks of the matrix. It is pretty easy to do this using threads but the SIMD part (finer level of //ism) is more difficult and implies to adapt (interleave) the data layout of the inputs (matrix and vectors).

I understand that this C++ code is not trivial to reproduce/translate in Julia but I think that it may provide a rather solid upper bound for performance (answer to your initial question about potential performance improvement).

---

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [March 16, 2025, 1:14pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/13 "2025-03-16T13:14:09Z")

</div>

> [@stevengj](#):
>
> Using `*` instead of `.*` means that the multiplication is not fused with the `.-=` — it allocates a separate array, in a separate loop.

I know Fuse vectorized operations. But here `*` is the _vector product_ of the matrix `L` and the vector `b`, not the _scalar product_ (`.*`). Replacing vector product with scalar product will lead to wrong result.

> [@stevengj](#):
>
> Also, `b[n*(i-1)+1:n*i]` allocates a new array

I know this, but I also know that array views don’t participate in calculations as quickly as array slices (which create new arrays). So I’m hesitant to choose between slices or views here.  
But I just tested it and you are right, it’s faster to use the view here.

> [@stevengj](#):
>
> Probably best to just put `@views` in front of `function`

Thank you very much for this suggestion, it cut my allocation in half again (from `325.692 μs Memory estimate: 114.73 KiB, allocs estimate: 893` to `308.563 μs Memory estimate: 31.06 KiB, allocs estimate: 298`) Awesome!

> [@stevengj](#):
>
> rather than just allocating a single buffer and re-using it?

Yes, I simply didn’t think of that, thanks for the suggestion. I’ve revised it.

> [@stevengj](#):
>
> Not using `mul!` with a buffer?

Yes, I added it later and it’s an area that could be improved.

Now the code is highly optimized.

```julia
@views function block_thomas_tridiagonal!(x::AbstractVector, L::AbstractVector{T}, D::AbstractVector{T},
    U::AbstractVector{T}, b::AbstractVector) where {T}
    n = size(D[1], 1) # each block has size n*n
    m = length(D) # num of blocks
    (; lu_buf, D_buf, b_buf) = get_block_thomas_buffer(n, m)

    # Forward elimination
    for i = 1:m-1
        lu_buf[i] = lu!(D[i])
        rdiv!(L[i], lu_buf[i])
        D[i+1] .-= mul!(D_buf, L[i], U[i])
        b[n*i+1:n*(i+1)] .-= mul!(b_buf, L[i], b[n*(i-1)+1:n*i])
    end

    # Backward substitution
    ldiv!(x[n*(m-1)+1:n*m], lu!(D[m]), b[n*(m-1)+1:n*m])
    for i = m-1:-1:1
        ldiv!(x[n*(i-1)+1:n*i], lu_buf[i],
            b[n*(i-1)+1:n*i] - mul!(b_buf, U[i], x[n*i+1:n*(i+1)]))
    end

    return x
end

```

with a buffer system:

```julia
mutable struct BlockThomas
    n::Int
    m::Int
    lu_buf::Vector{LU{Float64,Matrix{Float64},Vector{Int64}}}
    D_buf::Matrix{Float64}
    b_buf::Vector{Float64}
    function BlockThomas(n, m)
        lu_buf = Vector{LU{Float64,Matrix{Float64},Vector{Int64}}}(undef, m)
        D_buf = Matrix{Float64}(undef, n, n)
        b_buf = Vector{Float64}(undef, n)
        return new(n, m, lu_buf, D_buf, b_buf)
    end
end

const BlockThomasBuf = Vector{BlockThomas}(undef, nthreads())

function get_block_thomas_buffer(n::Int, m::Int)
    tid = threadid()
    return if !isassigned(BlockThomasBuf, tid)
        BlockThomasBuf[tid] = BlockThomas(n, m)
    elseif BlockThomasBuf[tid].n != n || BlockThomasBuf[tid].m != m
        BlockThomasBuf[tid] = BlockThomas(n, m)
    else
        BlockThomasBuf[tid]
    end
end

```

Thanks for your help stevengj.

---

<div class="post-metadata">

**Author:** ![karei](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/karei/32/214809_2.png) [@karei](https://discourse.julialang.org/u/karei)\
**Post date:** [March 16, 2025, 1:20pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/14 "2025-03-16T13:20:31Z")

</div>

Thank you for your detailed response. It is a bit difficult for me to understand what you are talking about 😢 and I don’t feel that I can do the simd version at my level.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [March 16, 2025, 1:28pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/15 "2025-03-16T13:28:46Z")

</div>

No worry 😉 . Like I said it may only be useful if you really want to get an upper bound for the performance (I put a lot of work into this since it was an essential part of a large industrial software).

A Julia translation of Legolas++ has been on my mind for a few years now (discussed this with @Elrod some time ago), but I could not find the time yet …

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 16, 2025, 4:54pm UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/16 "2025-03-16T16:54:18Z")

</div>

> [@karei](#):
>
> But here `*` is the _vector product_ of the matrix `L` and the vector `b`,

In that case you could use `mul!` with a pre-allocated buffer.

> [@karei](#):
>
> I know this, but I also know that array views don’t participate in calculations as quickly as array slices (which create new arrays).

This shouldn’t be a problem for _contiguous_ views like the ones here.

---

<div class="post-metadata">

**Author:** ![randy854](https://avatars.discourse-cdn.com/v4/letter/r/71c47a/32.png) [@randy854](https://discourse.julialang.org/u/randy854)\
**Post date:** [March 17, 2025, 4:02am UTC](https://discourse.julialang.org/t/how-to-improve-the-thomas-algorithm-for-block-tridiagonal-matrices/126894/17 "2025-03-17T04:02:38Z")

</div>

> [@karei](#):
>
> I did some searching but it seems that Julia doesn’t implement a specific algorithm for block tridiagonal matrices. So I implemented one myself. I would like to know:
> 
> 1. Have I missed some existing algorithm implementations for block tridiagonal matrices?
> 2. is there any room for improvement of the solver function I wrote?
> 
> ```julia
> function block_thomas_tridiagonal(A::AbstractVector{T}, B::AbstractVector{T}, C::AbstractVector{T},
> d::AbstractVector) where {T}
> B = copy(B)
> n = size(B[1], 1) # rows of each block
> m = length(B) # number of blocks
> 
> # Forward elimination
> for i = 2:m # (no speed up using inbounds+simd)
> multiplier = A[i-1] * inv(B[i-1])
> B[i] -= multiplier * C[i-1]
> d[n*(i-1)+1:n*i] .-= multiplier * d[n*(i-2)+1:n*(i-1)]
> end
> 
> # Backward substitution
> x = similar(d)
> x[n*(m-1)+1:n*m] = B[m] \ d[n*(m-1)+1:n*m] # no speed up using .=
> for i = m-1:-1:1
> x[n*(i-1)+1:n*i] = B[i] \ (d[n*(i-1)+1:n*i] - C[i] * x[n*i+1:n*(i+1)])
> end
> 
> return x
> end
> 
> ```
> 
> This function is specifically designed to handle the case where every block has the same size. It usually takes a Vector{Matrix{Float64}} as input and solves each block one by one.

Your block tridiagonal solver is a good start; improve it by using LU factorization for efficiency, ensuring type stability, adding error handling, and considering SIMD/multithreading for further optimization, while maintaining clear documentation and thorough testing.