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.
