# How to solve block tridiagonal matrix using Thomas algorithm

Thomas algorithm can be used to solve a tridiagonal matrix:

$$\begin{bmatrix} {b_ 1} & {c_ 1} & { } & { } & { 0 } \\ {a_ 2} & {b_ 2} & {c_ 2} & { } & { } \\ { } & {a_ 3} & {b_ 3} & \ddots & { } \\ { } & { } & \ddots & \ddots & {c_{n-1}}\\ { 0 } & { } & { } & {a_n} & {b_n}\\ \end{bmatrix} \cdot \begin{bmatrix} {x_ 1 } \\ {x_ 2 } \\ {x_ 3 } \\ \vdots \\ {x_n } \\ \end{bmatrix} = \begin{bmatrix} {d_ 1 } \\ {d_ 2 } \\ {d_ 3 } \\ \vdots \\ {d_n } \\ \end{bmatrix}$$

The solution can be get by

$$c'_i= \begin{cases} \begin{array}{lcl} \cfrac{c_i}{b_i} &&; i = 1 \\ \cfrac{c_i}{b_i - c'_{i - 1} a_i} &&; i = 2, 3, \ \dots, n-1 \\ \end{array} \end{cases}$$

$$d'_i= \begin{cases} \begin{array}{lcl} \cfrac{d_i}{b_i} &&; i = 1 \\ \cfrac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i} &&; i = 2, 3, \dots, n. \\ \end{array} \end{cases}$$

$$\begin{array}{lcl} x_n&=&d'_n\\ x_i&=&d'_i-c'_ix_{i+1} \end{array}$$

How to generalize this algorithm to block tridiagonal matrix, ie. both a,b,c are small square matrix itself. Do I just replace "times" by "dot" and "divide" by "dot inverse" in the above formula, and everything is OK?

I think it will be too much time consuming, if I try to give you answer and to correct superscripts,subscripts,capital and small letters of your query. Instead of that you can refer a detailed description of Thomas Algorithm from a book "Computational Fluid Dynamics" by John Anderson. You can directly refer page number "534" for it (Publisher: McGraw-Hill Science/Engineering/Math; 1st edition (February 1, 1995)). Its very easy to understand and implement it from there. All the best.

Yes. Except it's not "dot" but instead the product of matrices, resulting in another matrix.

• Thanks for the answer, but what is the correct order of these matrix operations? For example, in the above formula we have $c'_i=\cfrac{c_i}{b_i - c'_{i - 1} a_i}$, should it be translated to $c_i\cdot(b_i - c'_{i - 1}\cdot a_i)^{-1}$ or $c_i\cdot(b_i - a_i\cdot c'_{i - 1})^{-1}$ or $(b_i - c'_{i - 1}\cdot a_i)^{-1}\cdot c_i$ ? – xslittlegrass Mar 28 '13 at 3:59
• It depends. If you follow the derivation of the algorithm, you typically multiply each row of the matrix with a number so that, for example, you get a one on the diagonal (or something else you need to make the next step work). In each step, you will have to figure out what that means in the block matrix case. I suppose in most cases to multiply a block row you will need to multiply from the left by an inverse matrix. Just what it is you'll need to figure out by following how the algorithm works. – Wolfgang Bangerth Mar 28 '13 at 12:27

Maybe those subroutines will help you. They are not as efficient as I would like and one of the LAPACK calls simply destroys the solution (the commented one) but as it is below, it works!

program testando

implicit none

integer(4) :: i,j
real(8), dimension(2,2,3) :: A
real(8), dimension(2,2,3) :: B,C
real(8), dimension(2,3) :: xb
real(8), dimension(2,3) :: x

do i = 1,3
A(1,1,i) = 4.0d0
A(1,2,i) = 3.0d0
A(2,1,i) = 3.0d0
A(2,2,i) = 2.0d0
end do

do i = 1,3

B(1,1,i) = 8.0d0
B(1,2,i) = 6.0d0
B(2,1,i) = 6.0d0
B(2,2,i) = 4.0d0

C(1,1,i) = 16.0d0
C(1,2,i) = 9.0d0
C(2,1,i) = 9.0d0
C(2,2,i) = 6.0d0
end do

do i = 1, 2
do j = 1, 3
xb(i,j) = 3.0d0
end do
end do

write(*,*) "x(1) = 5.8571429"
write(*,*) "x(2) = -8.9220779"
write(*,*) "x(3) = 0.5714286"
write(*,*) "x(4) = -0.3116883"
write(*,*) "x(5) = 1.8571429"
write(*,*) "x(6) = -2.3766234"

write(*,*) "Calculated..."

do j = 1, 3
do i = 1, 2
write(*,*) "x(",i,") = ",x(i,j)
end do
end do
end program testando

!|     B(1)    C(1)             | | xb(1)  |
!| A(2)  B(2)    C(2)           | |        |
!|   A(3)  B(3)                 | |        |
!|     .     .                  |*|        | = B[1:mb*3]
!|       .     .                | |        |
!|               .       C(mb-1)| |        |
!|         A(mb)  B(mb)         | |xb(n*id)|

! id = inner matrices dimension.
! md = number matrices.
! maind = main diagonal of matrices  format: maind(id,id,md)
! lower = lower diagonal of matrices format: lower(id,id,2:md)
! upper = upper diagonal of matrices format: maind(id,id,md-1)
! xb    = B vector in Ax=B           format: xb(md*id)
! x     = x vector in Ax=B           format: x(md*id)
!

implicit none

! +++ Inputs +++
!
! Scalar input variables.
integer(kind=4) :: id,md

! Main diagonal of matrices.
real(kind=8), dimension(id,id,md) :: maind

! Lower diagonal of matrices.
real(kind=8), dimension(id,id,md) :: lower

! Upper diagonal of matrices.
real(kind=8), dimension(id,id,md) :: upper

! Vector of equalties B in Ax=B.
real(kind=8), dimension(id,md)    :: xb

! Vector of answers x in Ax=B.
real(kind=8), dimension(id,md)    :: x

! ++ Inside variables ++
!
! Scalar variables
integer(kind=4) :: i,ii,jj

! Array of gamma coefficients.
real(kind=8), dimension(id,id,md) :: gamm

! Array of beta coefficients.
real(kind=8), dimension(id,md) :: beta

! Auxiliary arrays.
real(kind=8), allocatable, dimension(:,:) :: aux_copy
real(kind=8), allocatable, dimension(:,:) :: aux_mult
real(kind=8), allocatable, dimension(:,:) :: aux_summ
real(kind=8), allocatable, dimension(:)   :: aux_dumm

allocate(aux_mult(id,id))
allocate(aux_summ(id,id))
allocate(aux_copy(id,id))
allocate(aux_dumm(id))

!--------------------------------------------------------------------------!
!                  Step 1: BLOCK TRIANGULARIZATION                         !
!--------------------------------------------------------------------------!

! Lets first get our first gamma.

call dlacpy('A',id,id,maind(:,:,1),id,aux_copy,id)

call inv(aux_copy,aux_copy,id)

call dgemm('N','N',id,id,id,1.0d0,aux_copy,id, &
upper(:,:,1),id,1.0d0,gamm(:,:,1),id)

! Now that we have our first gamma, let's get the rest of them.

do i = 2, md-1

call dgemm('N','N',id,id,id,1.0d0,lower(:,:,i),id, &
gamm(:,:,i-1),id,1.0d0,aux_mult,id)

do jj = 1, id
do ii = 1, id
aux_summ(ii,jj) = maind(ii,jj,i) - aux_mult(ii,jj)
end do
end do

call inv(aux_summ,aux_summ,id)

call dgemm('N','N',id,id,id,1.0d0,aux_summ,id, &
upper(:,:,i),id,1.0d0,gamm(:,:,i),id)

end do

! Now that we have our gammas, let's get the betas, starting from the first
! ones. Note that now the calls done by the Lapack library will get a bit
! more complicated so let's use matmul...

call dlacpy('A',id,id,maind(:,:,1),id,aux_copy,id)

call inv(aux_copy,aux_copy,id)

beta(:,1) = matmul(aux_copy,xb(:,1))

! We now have our first beta, lets get the rest.

do i = 2, md

! This LAPACK call does not work .... I dont know why.
!call dgemm('N','N',id,id,id,1.0d0,lower(:,:,i),id, &
!    gamm(:,:,i-1),id,1.0d0,aux_mult(:,:),id)

aux_mult(:,:) = matmul(lower(:,:,i),gamm(:,:,i-1))

do jj = 1, id
do ii = 1, id
aux_summ(ii,jj) = maind(ii,jj,i) - aux_mult(ii,jj)
end do
end do

call inv(aux_summ,aux_summ,id)

aux_dumm(:) = xb(:,i) - matmul(lower(:,:,i),beta(:,i-1))

beta(:,i) = matmul(aux_summ(:,:),aux_dumm(:))

end do

!--------------------------------------------------------------------------!
!                  Step 2: BACKWARD SWEEP                                  !
!--------------------------------------------------------------------------!

! How cool is that, lets start build our solution vector... iupiiii!

x(:,md) = beta(:,md)

do i = md-1,1,-1

aux_dumm(:) = matmul(gamm(:,:,i),x(:,i+1))

do ii = 1, id
x(ii,i) = beta(ii,i) - aux_dumm(ii)
end do

end do

deallocate(aux_mult)
deallocate(aux_summ)
deallocate(aux_copy)
deallocate(aux_dumm)

subroutine inv(A,A_inv,m)

Implicit none
integer :: m
real(8), dimension(m,m)::A, A_inv
real(8),dimension(m)::WORK
integer,dimension(m)::IPIV
integer info

A_inv = A

call DGETRF(M,M,A_inv,M,IPIV,info)

if (info /=  0) then
write(*,*)"DGETRF: Failed during matrix factorization"
stop
end if

call DGETRI(M,A_inv,M,IPIV,WORK,M,info)

if (info /=  0) then
write(*,*)"DGETRI: Failed during matrix inversion."
stop
end if

end subroutine inv


You have to link LAPACK to run it. With gfortran and a proper LAPACK installation in OpenSUSE the line is:

gfortran -g -llapack -lblas <program> -o block


Any tips on how to improve the performance of this subroutine will be welcome as well as why the LAPACK call is not working.

Yes as @wolfgang_bangerth says in his answer, you simply replace division by multiplication of matrix inverse and times by matrix multiplication. Things get interesting when the blocks become comparable or larger than available memory in which case the block operations can be parallelized by ScaLapack or other distributed linear algebra LU decomposition. For serial solvers, individual LU decompositions of the blocks can be written to disk to during forward phase of the solve to create an out-of-core solver that get's around memory limitations - at least for a some scales. Even more interesting, the tri-diagonal solver can be done parallelized itself with a block-cyclic technique while the block operations can or can not be parallelized and this can result in some truly fast solvers - see this paper for details.

No, you cannot use Thomas algorithm for block matrixes in general.

Go to wikipedia http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm, you will see that this formula is also there, but this form is derived from form with ~, ~ form may be used for block tridiagonal problems in some special cases.

Thomas algorithm is basically Gauss elimination method, problem with block matrixes is that product of two matrixes is not commutative (in general), so Gauss elimination can't work, for block problem you can only multiply from the left!

If Ai block are commutative with Bi and Ci and Bi~ blocks (this is the greatest problem since it would be hard to check commutativity because you derive Bi blocks in every iteration) then you can use Thomas ~ algorithm explained on wikipedia (I think you can't use form above because that form is derived from ~ form by multiplying every equation in matrix with 1/Bi, when you deal with matrix blocks that approach is not safe, it is better to use ~, it's not hard to derive formulas for ~ form).

Cases when this method can be used for tridiagonal block matrix problem are quite rare, because of commutativity with Bi~ blocks is needed. Luckily for you many differential equations have Ai i Ci blocks in a form const*eye(n)=const * I (in many cases Ai=Ci) and eye matrix is commutative with everything :D

• I don't see why you need commutativity (although you do need to worry about the difference between left- and right-inverses). The Thomas algorithm is just a special case of the LU decomposition, for which block versions exist. – Christian Clason Aug 12 '14 at 17:02