CFD Online Discussion Forums - Subroutine for Block tridiagonal System

Here is the routine out of the Tannehill book "Computational Fluid Mechanics and Heat Transfer"

This will be an exact copy, and it looks like fortran77

c..subroutine to solve non-period block tridiaganol c..system of equations without pivoting strategy c..with the dimension of the block matrices being c..n x n (n is any number greater than 1)

SUBROUTINE NBTRIP(A,B,C,D,IL,IU,ORDER) INTEGER ORDER,ORDSQ DIMENSION A(1),B(1),D(1)

C..A = SUB DIAG. MATRIX C..B = DIAG. MATRIX C..C = SUP DIAG. MATRIX C..D = RIGHT HAND SIDE OF VECTOR C..IL = LOWER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED C..IU = UPPER VALUE OF ' ' ' ' C.. (SOLUTION IS SOUGHT FOR BTRI(A,B,C)*X=D C.. FOR INDICES OF X BETWEEN IL AND IU (INCLUSIVE). C.. SOLUTION WRITTEN IN D VECTOR (ORIGINAL CONTENTS C.. ARE OVERWRITTEN)). C..ORDER = ORDER OF A,B,C MATRICES AND LENGTH OF D VECTOR C.. AT EACH POINT DENOTED BY INDEX I C.. C..THE MATRICES AND VECTORS ARE STORED IN SINGLE C..SUBSCRIPT FORM

ORDSQ = ORDER**2

C.. C.. FORWARD ELIMINATION C..

I=IL

IOMAT = 1+(I-1)*ORDSQ

IOVEC = 1+(I-2)*ORDER

CALL LUDECO(B(IOMAT),ORDER)

CALL LUSOLV(B(IOMAT),D(IOVEC),D(IOVEC),ORDER)

DO 100 J=1,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),C(IOMAT),C(IOMATJ),ORDER) 100 CONTINUE 200 CONTINUE

I=I+1

IOMAT = 1+(I-1)*ORDSQ

IOVEC = 1+(I-1)*ORDER

I1MAT = IOMAT-ORDSQ

I1VEC = IOVEC-ORDER

CALL MULPUT(A(IOMAT),D(I1VEC),D(IOVEC),ORDER)

DO 300 J=1,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

I1MATJ = I1MAT+(J-1)*ORDER

CALL MULPUT(A(IOMAT),C(I1MATJ),B(IOMATJ),ORDER) 300 CONTINUE

CALL LUDECO(B(IOMAT),ORDER)

CALL LUSOLV(B(IOMAT),D(IOVEC),D(IOVEC),ORDER)

IF(I.EQ.IU) GOTO 500

DO 400 J=1,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),C(IOMATJ),C(IOMATJ),ORDER) 400 CONTINUE

GOTO 200 500 CONTINUE

C.. C.. BACK SUBSTITUTION C..

I=IU 600 CONTINUE

I=I-1

IOMAT=1+(I-1)*ORDSQ

IOVEC=1+(I-1)*ORDER

I1VEC=IOVEC+ORDER

CALL MULPUT(C(IOMAT),D(I1VEC),D(IOVEC),ORDER)

IF(I.GT.IL) GOTO 600 C..

RETURN

END

C.. SUBROUTINE TO CALCULATE L-U DECOMPOSITION C.. OF A GIVEN MATRIX A AND STORE RESULT IN A

SUBROUTINE LUDECO(A,ORDER)

INTEGER ORDER

DIMENSION A(ORDER,1)

DO 8 JC=2,ORDER 8 A(1,JC) = A(1,JC)/A(1,1)

JRJC = 1 10 CONTINUE

JRJC = JRJC+1

JRJCM1 = JRJC-1

JRJCP1 = JRJC+1

DO 14 JR=JRJC,ORDER

SUM = A(JR,JRJC)

DO 12 JM=1,JRJCM1 12 SUM = SUM-A(JR,JM)*A(JM,JRJC) 14 A(JR,JRJC) = SUM

IF (JRJC.EQ.ORDER) RETURN

DO 18 JC=JRJCP1,ORDER

SUM = A(JRJC,JC)

DO 16 JM=1,JRJCM1 16 SUM = SUM-A(JRJC,JM)*A(JM,JC) 18 A(JRJC,JC) = SUM/A(JRJC,JRJC)

GOTO 10

END

C.. SUBROUTINE TO MULTIPLY A VECTOR B BY A MATRIX A C.. SUBTRACT RESULT FROM ANOTHER VECTOR C AND STORE C.. RESULT IN C. THUS VECTOR C IS OVERWRITTEN

SUBROUTINE MULPUT(A,B,C,ORDER)

INTEGER ORDER

DIMENSION A(1),B(1),C(1)

DO 200 JR=1,ORDER

SUM = 0.0

DO 100 JC=1,ORDER

IA = JR+(JC-1)*ORDER 100 SUM = SUM+A(IA)*B(JC) 200 C(JR) = C(JR)-SUM

RETURN

END

C.. SUBROUTINE TO SOLVE LINEAR ALGEBRAIC SYSTEM OF C.. EQUATION A*C=B AND STORE RESULTS IN VECTOR C C.. MATRIX A INPUT IN L-U DECOMPOSITION FORM

SUBROUTINE LUSOLV(A,B,C,ORDER)

INTEGER ORDER

DIMENSION A(ORDER,1),B(1),C(1)

C.. FIRST L(INV)*B

C(1) = C(1)/A(1,1)

DO 14 JR=2,ORDER

JRM1 = JR-1

SUM = B(JR)

DO 12 JM=1,JRM1 12 SUM = SUM-A(JR,JM)*C(JM) 14 C(JR) = SUM/A(JR,JR)

C.. NEXT U(INV) OF L(INV)*B

DO 18 JRJR=1,ORDER

JR = ORDER,JRJR+1

JRP1 = JR+1

SUM = C(JR)

DO 16 JMJM = JRP1,ORDER

JM = ORDER-JMJM+JRP1 16 SUM = SUM-A(JR,JM)*C(JM) 18 C(JR) = SUM

RETURN

END

I guess if I get really bored tomorrow, I'll post the periodic subroutine, which uses the same sub..subroutines.