7
$\begingroup$

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?

Improve this question
$\endgroup$
2

5 Answers 5

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
0
7
$\begingroup$

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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$ ? $\endgroup$
    – xslittlegrass
    Commented Mar 28, 2013 at 3:59
  • $\begingroup$ 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. $\endgroup$
    – Wolfgang Bangerth
    Commented Mar 28, 2013 at 12:27
1
$\begingroup$

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


call cc299blktriad(A,B,C,2,3,xb,x)

! answer
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

subroutine cc299blktriad(maind,lower,upper,id,md,xb,x)

 !|     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)

end subroutine cc299blktriad


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.

Improve this answer
$\endgroup$
0
$\begingroup$

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.

Improve this answer
$\endgroup$
-3
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – Christian Clason
    Commented Aug 12, 2014 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.