2
$\begingroup$

In 1D problems, tridiagonal systems of equations are obtained when we use finite-difference or finite-volumes in a structured mesh. A wide solver is the TDMA algorithm here. In two-dimensional modelling, we achieve a block tridiagonal systems of equations. However, I found a lack of codes dealing with those block matrices. In book:

Pletcher R.H., Tannehill J.C., Anderson D.A. Computational Fluid Mechanics and Heat Transfer (series in computational and physical processes in mechanics and thermal sciences). 3rd ed. CRC Press Taylor & Francis Group, 2013.

there are the subroutines for making up a block-TDMA, code in Fortran. I transcript almost the entire appendix B:

Subroutine NBTRIP solves a block tridiagonal system of equations of the form $$\begin{pmatrix}B_{IL}&C_{IL}&&\\A_{I}&B_{I}&C_{I}&&\\&\ddots&\ddots&\ddots\\&&A_{IU}&B_{IU}\end{pmatrix}\begin{pmatrix}X_{IL}\\X_{I}\\\vdots\\X_{IU}\end{pmatrix} = \begin{pmatrix}D_{IL}\\D_{I}\\\vdots\\D_{IU}\end{pmatrix}$$ Subroutine PBTRIP solves a periodic block tridiagonal system of equations in the form $$\begin{pmatrix}B_{IL}&C_{IL}&&A_{IL}\\A_{I}&B_{I}&C_{I}&&\\&\ddots&\ddots&\ddots\\C_{IU}&&A_{IU}&B_{IU}\end{pmatrix}\begin{pmatrix}X_{IL}\\X_{I}\\\vdots\\X_{IU}\end{pmatrix} = \begin{pmatrix}D_{IL}\\D_{I}\\\vdots\\D_{IU}\end{pmatrix}$$ The block matrices $A$, $B$, and $C$ are $N × N$ matrices at every point $I$ with $N$ being an integer greater than 1. The right-hand side vector $D$ has length $N$ at each point $I$.Therefore, the total number of $I$ points at which the matrices are defined (denoted by $NI$) is given by $NI = IU - IL + 1$. The matrices $A$, $B$, and $C$ are dimensioned as $A(N,N,NI)$, $B(N,N,NI)$ and $C(N,N,NI)$ while the vector $D$ is dimensioned as $D(N,NI)$. The solution, $X$, is returned to the calling program by overwriting the $D$ vector with the $X$ vector.

FORTRAN 90 code:

!++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                                                      !
! SUBROUTINE TO SOLVE NON-PERIODIC BLOCK TRIDIAGONAL   !
! SYSTEM OF EQUATIONS WITHOUT PIVOTING STRATEGY        !
! WITH THE DIMENSIONS OF THE BLOCK MATRICES BEING      !
! N × N (N IS ANY NUMBER GREATER THAN 1).              !
!                                                      !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++!
SUBROUTINE NBTRIP (A, B, C, D, IL, IU, ORDER)
  INTEGER, INTENT(IN) :: IL, IU, ORDER
  REAL, INTENT(INOUT) :: A(1), B(1)
  REAL, INTENT(INOUT) :: C(1), D(1)
!...
!...A = SUB DIAGONAL MATRIX
!...B = DIAGONAL MATRIX
!...C = SUP DIAGONAL MATRIX
!...D = RIGHT HAND SIDE VECTOR
!...IL = LOWER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED
!...IU = UPPER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED
!... (SOLUTION IS SOUGHT FOR BTRI(A, B, C)*X = D
!...  FOR INDICES OF X BETWEEN IL AND IU (INCLUSIVE).
!...  SOLUTION WRITTEN IN D VECTOR (ORIGINAL CONTENTS
!...  ARE OVERWRITTEN)).
!...ORDER = ORDER OF A, B, C MATRICES AND LENGTH OF D VECTOR
!...        AT EACH POINT DENOTED BY INDEX I
!... (ORDER CAN BE ANY INTEGER GREATER THAN 1).
!...
!...THE MATRICES AND VECTORS ARE STORED IN SINGLE SUBSCRIPT FORM
!...
  INTEGER :: ORDSQ
  INTEGER :: I0MAT, I0MATJ, I0VEC
  INTEGER :: I1MAT, I1MATJ, I1VEC
  INTEGER :: I, J
!...
!...FORWARD ELIMINATION
!...
  ORDSQ = ORDER**2
  I = IL
  I0MAT = 1 + (I - 1)*ORDSQ
  I0VEC = 1 + (I - 1)*ORDER
  CALL LUDECO (B(I0MAT), ORDER)
  CALL LUSOLV (B(I0MAT), D(I0VEC), D(I0VEC), ORDER)
  DO J = 1, ORDER
    I0MATJ = I0MAT + (J - 1)*ORDER
    CALL LUSOLV (B(I0MAT), C(I0MATJ), C(I0MATJ), ORDER)
  END DO
  DO
    I = I + 1
    I0MAT = 1 + (I - 1)*ORDSQ
    I0VEC = 1 + (I - 1)*ORDER
    I1MAT = I0MAT - ORDSQ
    I1VEC = I0VEC - ORDER
    CALL MULPUT (A(I0MAT), D(I1VEC), D(I0VEC), ORDER)
    DO J = 1, ORDER
      I0MATJ = I0MAT + (J - 1)*ORDER
      I1MATJ = I1MAT + (J - 1)*ORDER
      CALL MULPUT (A(I0MAT), C(I1MATJ), B(I0MATJ), ORDER)
    END DO
    CALL LUDECO (B(I0MAT), ORDER)
    CALL LUSOLV (B(I0MAT), D(I0VEC), D(I0VEC), ORDER)
    IF(I == IU) EXIT
    DO J = 1, ORDER
      I0MATJ = I0MAT + (J - 1)*ORDER
      CALL LUSOLV (B(I0MAT), C(I0MATJ), C(I0MATJ), ORDER)
    END DO
  END DO
!...
!...BACK SUBSTITUTION
!...
  I = IU
  DO
    I = I - 1
    I0MAT = 1 + (I - 1)*ORDSQ
    I0VEC = 1 + (I - 1)*ORDER
    I1VEC = I0VEC + ORDER
    CALL MULPUT(C(I0MAT), D(I1VEC), D(I0VEC), ORDER)
    IF (I <= IL) EXIT
  END DO
!...
END SUBROUTINE NBTRIP
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                                                       !
! SUBROUTINE TO SOLVE PERIODIC BLOCK TRIDIAGONAL        !
! SYSTEM OF EQUATIONS WITHOUT PIVOTING STRATEGY.        !
! EACH BLOCK MATRIX MAY BE OF DIMENSION N WITH          !
! N ANY NUMBER GREATER THAN 1.                          !
!                                                       !
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++!
SUBROUTINE PBTRIP (A, B, C, D, IL, IU, ORDER)
  INTEGER, INTENT(IN) :: IL, IU, ORDER
  REAL, INTENT(INOUT) :: A(1), B(1)
  REAL, INTENT(INOUT) :: C(1), D(1)
!...
!...A = SUB DIAGONAL MATRIX
!...B = DIAGONAL MATRIX
!...C = SUP DIAGONAL MATRIX
!...D = RIGHT HAND SIDE VECTOR
!...IL = LOWER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED
!...IU = UPPER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED
!... (SOLUTION IS SOUGHT FOR BTRI(A, B, C)*X = D
!...  FOR INDICES OF X BETWEEEN IL AND IU (INCLUSIVE).
!...  SOLUTION WRITTEN IN D VECTOR (ORIGINAL CONTENTS
!...  ARE OVERWRITTEN)).
!...ORDER = ORDER OF A, B, C MATRICES AND LENGTH OF D VECTOR
!...        AT EACH POINT DENOTED BY INDEX I
!... (ORDER CAN BE ANY INTEGER GREATER THAN 1)
!...
!...
!...THE MATRICES AND VECTORS ARE STORED IN SINGLE SUBSCRIPT FORM
!...
  INTEGER :: ORDSQ
  INTEGER :: IS, IE, IUMAT, IUVEC, IEMAT, IEVEC
  INTEGER :: I0MAT, I0VEC, I1MAT, I1VEC
  INTEGER :: I0MATJ, I0VECJ, I1MATJ, I1VECJ, IUMATJ, IEMATJ
  INTEGER :: I, J, IBAC
  REAL :: AD(1), CD(1)
!...
  IS = IL + 1
  IE = IU - 1
  ORDSQ = ORDER**2
  IUMAT = 1 + (IU - 1)*ORDSQ
  IUVEC = 1 + (IU - 1)*ORDER
  IEMAT = 1 + (IE - 1)*ORDSQ
  IEVEC = 1 + (IE - 1)*ORDER
!...
!...FORWARD ELIMINATION
!...
  I = IL
  I0MAT = 1 + (I - 1)*ORDSQ
  I0VEC = 1 + (I - 1)*ORDER
  CALL LUDECO (B(I0MAT), ORDER)
  CALL LUSOLV(B(I0MAT), D(I0VEC), D(I0VEC), ORDER)
  DO J = 1, ORDER
    I0MATJ = I0MAT + (J - 1)*ORDER
    CALL LUSOLV (B(I0MAT), C(I0MATJ), C(I0MATJ), ORDER)
    CALL LUSOLV (B(I0MAT), A(I0MATJ), A(I0MATJ), ORDER)
  END DO
!...
  DO I = IS, IE
    I0MAT = 1  +(I - 1)*ORDSQ
    I0VEC = 1 + (I - 1)*ORDER
    I1MAT = I0MAT - ORDSQ
    I1VEC = I0VEC - ORDER
    DO J = 1, ORDSQ
      I0MATJ = J - 1 + I0MAT
      IUMATJ = J - 1 + IUMAT
      AD (J) = A(I0MATJ)
      CD (J) = C(IUMATJ)
      A(I0MATJ) = 0.0
      C(IUMATJ) = 0.0
    END DO
    CALL MULPUT (AD, D(I1VEC), D(I0VEC), ORDER)
    DO J = 1, ORDER
      I0MATJ = I0MAT + (J - 1)*ORDER
      I1MATJ = I1MAT + (J - 1)*ORDER
      CALL MULPUT (AD, C(I1MATJ), B(I0MATJ), ORDER)
      CALL MULPUT (AD, A(I1MATJ), A(I0MATJ), ORDER)
    END DO
    CALL LUDECO (B(I0MAT), ORDER)
    CALL LUSOLV (B(I0MAT), D(I0VEC), D(I0VEC), ORDER)
    DO J = 1, ORDER
      I0MATJ = I0MAT + (J - 1)*ORDER
      CALL LUSOLV (B(I0MAT), C(I0MATJ), C(I0MATJ), ORDER)
      CALL LUSOLV (B(I0MAT), A(I0MATJ), A(I0MATJ), ORDER)
    END DO
    CALL MULPUT (CD, D(I1VEC), D(IUVEC), ORDER)
    DO J = 1, ORDER
      IUMATJ = IUMAT + (J - 1)*ORDER
      I1MATJ = I1MAT + (J - 1)*ORDER
      CALL MULPUT (CD, A(I1MATJ), B(IUMATJ), ORDER)
      CALL MULPUT (CD, C(I1MATJ), C(IUMATJ), ORDER)
    END DO
  END DO
!...
  DO J = 1, ORDSQ
    IUMATJ = J - 1 + IUMAT
    AD(J) = A(IUMATJ) + C(IUMATJ)
  END DO
  CALL MULPUT (AD, D(IEVEC), D(IUVEC), ORDER)
  DO J = 1, ORDER
    IUMATJ = IUMAT + (J - 1)*ORDER
    IEMATJ = IEMAT + (J - 1)*ORDER
    CALL MULPUT (AD, C(IEMATJ), B(IUMATJ), ORDER)
    CALL MULPUT (AD, A(IEMATJ), B(IUMATJ), ORDER)
  END DO
  CALL LUDECO (B(IUMAT), ORDER)
  CALL LUSOLV (B(IUMAT), D(IUVEC), D(IUVEC), ORDER)
!...
!...BACK SUBSTITUTION
!...
  DO IBAC = IL, IE
    I = IE - IBAC + IL
    I0MAT = 1 + (I - 1)*ORDSQ
    I0VEC = 1 + (I - 1)*ORDER
    I1VEC = I0VEC + ORDER
    CALL MULPUT (A(I0MAT), D(IUVEC), D(I0VEC), ORDER)
    CALL MULPUT (C(I0MAT), D(I1VEC), D(I0VEC), ORDER)
  END DO
!...
END SUBROUTINE PBTRIP
!+++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                                                     !
! SUBROUTINE TO CALCULATE L-U DECOMPOSITION           !
! OF A GIVEN MATRIX A AND STORE RESULT IN A           !
! (NO PIVOTING STRATEGY IS EMPLOYED)                  !
!                                                     !
!+++++++++++++++++++++++++++++++++++++++++++++++++++++!
SUBROUTINE LUDECO (A, ORDER)
  INTEGER, INTENT(IN) :: ORDER
  INTEGER :: JR, JC, JM, JRJC, JRJCM1, JRJCP1
  REAL, INTENT(INOUT) :: A(ORDER, 1)
  REAL :: SUM
!...
  DO JC = 2, ORDER
    A(1, JC) = A(1, JC)/A(1,1)
  END DO
  JRJC = 1
  DO
    JRJC = JRJC + 1
    JRJCM1 = JRJC - 1
    JRJCP1 = JRJC + 1
    DO JR = JRJC, ORDER
      SUM = A(JR, JRJC)
      DO JM = 1, JRJCM1
        SUM = SUM - A(JR, JM)*A(JM, JRJC)
      END DO
      A(JR, JRJC) = SUM
    END DO
    IF (JRJC == ORDER) EXIT
    DO JC = JRJCP1, ORDER
      SUM = A(JRJC, JC)
      DO JM = 1, JRJCM1
        SUM = SUM - A(JRJC, JM)*A(JM, JC)
      END DO
      A(JRJC, JC) = SUM/A(JRJC, JRJC)
    END DO
  END DO
!...
END SUBROUTINE LUDECO
!++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                                                      !
! SUBROUTINE TO SOLVE LINEAR ALGEBRAIC SYSTEM OF       !
! EQUATIONS A*C=B AND STORE RESULTS IN VECTOR C.       !
! MATRIX A IS INPUT IN L-U DECOMPOSITION FORM.         !
! (NO PIVOTING STRATEGY HAS BEEN EMPLOYED TO           !
! COMPUTE THE L-U DECOMPOSITION OF THE MATRIX A).      !
!                                                      !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++!
SUBROUTINE LUSOLV (A, B, C, ORDER)
  INTEGER, INTENT(IN) :: ORDER
  INTEGER :: JR, JM, JRM1, JRP1, JRJR, JMJM
  REAL, INTENT(IN) :: A(ORDER, 1), B(1)
  REAL, INTENT(INOUT) :: C(1)
  REAL :: SUM
!...
!...FIRST L(INV)*B
!...
  C(1) = C(1)/A(1,1)
  DO JR = 2, ORDER
    JRM1 = JR - 1
    SUM = B(JR)
    DO JM = 1, JRM1
      SUM = SUM - A (JR, JM)*C(JM)
    END DO
    C(JR) = SUM/A(JR, JR)
  END DO
!...
!...NEXT U(INV) OF L(INV)*B
!...
  DO JRJR = 2, ORDER
    JR = ORDER - JRJR + 1
    JRP1 = JR + 1
    SUM = C(JR)
    DO JMJM = JRP1, ORDER
      JM = ORDER - JMJM + JRP1
      SUM = SUM - A(JR, JM)*C(JM)
    END DO
    C(JR) = SUM
  END DO
!...
END SUBROUTINE LUSOLV
!++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                                                      !
! SUBROUTINE TO MULTIPLY A VECTOR B BY A MATRIX A      !
! SUBTRACT RESULT FROM ANOTHER VECTOR C AND STORE      !
! RESULT IN C. THUS VECTOR C IS OVERWRITTEN.           !
!                                                      !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++!
SUBROUTINE MULPUT (A, B, C, ORDER)
  INTEGER, INTENT(IN) :: ORDER
  INTEGER :: JR, JC, IA
  REAL, INTENT(IN) :: A(1), B(1)
  REAL, INTENT(INOUT) :: C(1)
  REAL :: SUM
!...
  DO JR = 1, ORDER
    SUM = 0.0
    DO JC = 1, ORDER
      IA = JR + (JC - 1)*ORDER
      SUM = SUM + A(IA)*B(JC)
    END DO
    C(JR) = C(JR) - SUM
  END DO
!...
END SUBROUTINE MULPUT

However, someone advertised me that there is something wrong in LUDECO and LUSOLV subroutines that weren't fixed them. Have anyone used and tested this code? Is there any other block-TDMA in Fortran?

Thanks in advance.

$\endgroup$
2
$\begingroup$

Something between a comment and an answer – a couple of points and links to internal Computational Science resources that should be helpful and relevant.

In this question, the 3D finite-difference is discussed and it is pointed out that 2D and 3D discretizations are not actually block-tridiagonal. More details there.

This question discusses how to solve block-tridiagonal algorithms using the Thomas algorithm with links and even some Fortran code. And here, the derivation of the block algorithm is shown for a more complicated case.

So, I would reassess the structure of the matrix you have, feasibility of it being actually solved by a block-tridiagonal algorithm, and potentially resort to very efficient sparse solvers.

Unfortunately, I don't know anything about the subroutines you've listed above.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ The matrix for discrete 2D Poisson equation by numbering in the correct sequence is block-tridiagonal, for instance. $\endgroup$ – V.J. Dec 27 '19 at 13:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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