      SUBROUTINE TRDBLK(M,NBLOKS,A,B,C,IPIVOT,IFLAG,Y)
C
C***********************************************************************
C
C  THIS PROGRAM SOLVES THE LINEAR SYSTEM  T*X = Y, WHERE T IS BLOCK
C  TRIDIAGONAL OF THE FORM
C
C                    A(1)   B(1)
C
C                    C(1)   A(2)   B(2)
C
C                         .      .      .
C
C                             .      .       .
C
C                                 .      .      B(NBLOKS-1)
C
C                                 C(NBLOKS-1)   A(NBLOKS)
C
C  WHERE
C                   A(K), K = 1,NBLOKS, ARE  M  BY  M,
C                   B(K), K = 1,NBLOKS-1, ARE  M  BY  M,
C                   C(K), K = 1,NBLOKS-1, ARE  M  BY  M,
C
C
C  THE METHOD USED IS A VARIANT OF GAUSS ELIMINATION.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C               *****  PARAMETERS  *****
C
C       *** ON ENTRY ...
C
C               M  -  INTEGER
C                      THE ORDER OF EACH BLOCK
C
C          NBLOKS  -  INTEGER
C                      NUMBER OF DIAGONAL BLOCKS. (MUST BE GREATER
C                      THAN OR EQUAL TO 3)
C
C               A  -  DOUBLE PRECISION(M,M,NBLOKS)
C                      A(,,K) CONTAINS THE K-TH M BY M DIAGONAL BLOCK
C                      OF THE MATRIX T
C
C               B  -  DOUBLE PRECISION(M,M,NBLOKS-1)
C                      B(,,K) CONTAINS THE K-TH SUPERDIAGONAL BLOCK
C                      OF THE MATRIX T
C
C               C  -  DOUBLE PRECISION(M,M,NBLOKS-1)
C                      C(,K) CONTAINS THE K-TH SUBDIAGONAL BLOCK
C                      OF THE MATRIX T
C
C          IPIVOT  -  INTEGER(M,NBLOKS)
C                      WORK SPACE
C
C               Y  -  DOUBLE PRECISION(M,NBLOKS)
C                      THE RIGHT HAND SIDE VECTOR
C
C
C       ***  ON RETURN  ...
C
C         A, B, C  -  ARRAYS CONTAINING THE DESIRED BLOCK  LU
C                      DECOMPOSITION OF THE MATRIX T (IF IFLAG = 0)
C
C          IPIVOT  -  INTEGER(M,NBLOKS)
C                      IPIVOT(,K) RECORDS THE PIVOTING INDICES
C                      DETERMINED IN THE DECOMPOSITION OF K-TH
C                      DIAGONAL BLOCK.
C
C               Y  -  DOUBLE PRECISION(M,NBLOKS)
C                      THE SOLUTION VECTOR (IF IFLAG = 0)
C
C           IFLAG  -  INTEGER
C                      = -1, IF INPUT PARAMETERS ARE INVALID
C                      =  K, IF THE K-TH DIAGONAL BLOCK IS SINGULAR
C                      =  0, OTHERWISE
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C             *****  AUXILIARY PROGRAMS  *****
C
C         TRDFAC(M,NBLOKS,A,B,C,IPIVOT,IFLAG)
C              - DECOMPOSES THE MATRIX  T  USING BLOCK GAUSSIAN
C                ELIMINATION AND IS USED FOR THIS PURPOSE IN
C                T R D B L K.
C                THE ARGUMENTS ARE AS IN T R D B L K.
C
C         TRDSOL(M,NBLOKS,A,B,C,Y,IPIVOT)
C              - SOLVES THE SYSTEM T*X = Y ONCE T  IS DECOMPOSED.
C                THE ARGUMENTS ARE AS IN T R D B L K.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  THE SUBROUTINE  T R D B L K  AUTOMATICALLY SOLVES THE INPUT
C  SYSTEM WHEN IFLAG=0.  T R D B L K  IS CALLED ONLY ONCE FOR A
C  GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT HAND
C  SIDES CAN BE OBTAINED BY ONE CALL TO  T R D B L K  AND P-1
C  CALLS TO TRDSOL ONLY. SINCE THE ARRAYS A, B, C AND IPIVOT CONTAIN
C  THE DECOMPOSITION OF THE GIVEN COEFFICIENT MATRIX AND PIVOTING
C  INFORMATION ON RETURN FROM T R D B L K, THEY MUST NOT BE
C  ALTERED BETWEEN SUCCESSIVE CALLS TO T R D B L K WITH THE SAME
C  LEFT HAND SIDES. FOR THE SAME REASON, IF THE USER WISHES TO SAVE
C  THE COEFFICIENT MATRIX, THE ARRAYS A, B, C  MUST BE COPIED BEFORE A
C  CALL TO  T R D B L K .
C
C***********************************************************************
C
      DOUBLE PRECISION A,B,C,Y
      INTEGER M,NBLOKS,IPIVOT,IFLAG
      DIMENSION A(M,M,NBLOKS),B(M,M,NBLOKS-1),C(M,M,NBLOKS-1),
     *  IPIVOT(M,NBLOKS),Y(M,NBLOKS)
      CALL TRDFAC(M,NBLOKS,A,B,C,IPIVOT,IFLAG)
      IF(IFLAG.NE.0)  RETURN
      CALL TRDSOL(M,NBLOKS,A,B,C,Y,IPIVOT)
      RETURN
      END

      SUBROUTINE TRDFAC(M,NBLOKS,A,B,C,IPIVOT,IFLAG)
C
C***********************************************************************
C
C  T R D F A C  DECOMPOSES THE BLOCK TRIDIAGONAL MATRIX  T  USING A
C  BLOCK VERSION OF GAUSSIAN ELIMINATION.  THE MATRIX  T  IS STORED
C  IN THE ARRAYS  A, B, C.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C            *****  PARAMETERS  *****
C
C      *** ON ENTRY ...
C
C              M  -  INTEGER
C                     THE ORDER OF EACH BLOCK
C
C         NBLOKS  -  INTEGER
C                     NUMBER OF DIAGONAL BLOCKS. (MUST BE GREATER
C                     THAN OR EQUAL TO 3)
C
C              A  -  DOUBLE PRECISION(M,M,NBLOKS)
C                     A(,,K) CONTAINS THE K-TH M BY M DIAGONAL BLOCK
C                     OF THE MATRIX T
C
C              B  -  DOUBLE PRECISION(M,M,NBLOKS-1)
C                     B(,,K) CONTAINS THE K-TH SUPERDIAGONAL BLOCK
C                     OF THE MATRIX T
C
C              C  -  DOUBLE PRECISION(M,M,NBLOKS-1)
C                     C(,K) CONTAINS THE K-TH SUBDIAGONAL BLOCK
C                     OF THE MATRIX T
C
C         IPIVOT  -  INTEGER(M,NBLOKS)
C                     WORK SPACE
C
C
C      ***  ON RETURN  ...
C
C        A, B, C  -  ARRAYS CONTAINING THE DESIRED  BLOCK  LU
C                     DECOMPOSITION OF THE MATRIX T ( IF IFLAG = 0)
C
C         IPIVOT  -  INTEGER(M,NBLOKS)
C                     IPIVOT(,K) RECORDS THE PIVOTING INDICES DETERMINED
C                     IN THE DECOMPOSITION OF K-TH DIAGONAL BLOCK.
C
C         IFLAG  -  INTEGER
C                     = -1, IF INPUT PARAMETERS ARE INVALID
C                     =  K, IF THE K-TH DIAGONAL BLOCK IS SINGULAR
C                     =  0, OTHERWISE
C
C***********************************************************************
C
      DOUBLE PRECISION A,B,C,DP
      INTEGER M,NBLOKS,IPIVOT,IFLAG,I,J,K,L,KMIN1
      DIMENSION A(M,M,NBLOKS),B(M,M,NBLOKS-1),C(M,M,NBLOKS-1),
     *  IPIVOT(M,NBLOKS)
      IFLAG = 0
C
C***********************************************************************
C
C         ***  CHECK VALIDITY OF INPUT PARAMETERS ....
C
C              IF PARAMETERS ARE INVALID THEN TERMINATE.
C
C***********************************************************************
C
      IF (M .LT. 1 .OR. NBLOKS .LT. 3)  THEN
         IFLAG = -1
         RETURN
      ENDIF
C
C***********************************************************************
C
C         ***  PROCESS THE FIRST BLOCK ROW ....
C
C***********************************************************************
C
      CALL COLDEC(M,M,A,IPIVOT,IFLAG)
C
C     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C       CHECK FOR SINGULARITY OF FIRST DIAGONAL BLOCK:
C
C                 IF SINGULAR TERMINATE ....
C
C     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IF (IFLAG .NE. 0) THEN
         IFLAG = 1
         RETURN
      ENDIF
      DO 10 J = 1,M
         CALL COLSOL(M,M,A,B(1,J,1),IPIVOT)
10    CONTINUE
C
C***********************************************************************
C
C          ***  MAIN LOOP : PROCESS BLOCK ROWS 2 TO NBLOKS ....
C
C***********************************************************************
C
      DO 60 K = 2,NBLOKS
         KMIN1 = K-1
         DO 40 J = 1,M
            DO 30 L = 1,M
               DP = B(L,J,KMIN1)
               DO 20 I = 1,M
                  A(I,J,K) = A(I,J,K) - C(I,L,KMIN1)*DP
20             CONTINUE
30          CONTINUE
40       CONTINUE
         CALL COLDEC(M,M,A(1,1,K),IPIVOT(1,K),IFLAG)
C
C        CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C          CHECK FOR SINGULARITY OF K-TH DIAGONAL BLOCK:
C
C                    IF SINGULAR TERMINATE ....
C
C        CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
         IF (IFLAG .NE. 0) THEN
            IFLAG = K
            RETURN
         ENDIF
         IF(K .LT. NBLOKS) THEN
            DO 50 J = 1,M
               CALL COLSOL(M,M,A(1,1,K),B(1,J,K),IPIVOT(1,K))
50          CONTINUE
         ENDIF
60    CONTINUE
      RETURN
      END
      SUBROUTINE TRDSOL(M,NBLOKS,A,B,C,Y,IPIVOT)
C
C***********************************************************************
C
C  T R D S O L  SOLVES THE LINEAR SYSTEM
C
C                     T*X = Y
C
C  USING THE DECOMPOSITION ALREADY GENERATED IN T R D F A C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C            *****  PARAMETERS  *****
C
C      *** ON ENTRY ...
C
C              M  -  INTEGER
C                     THE ORDER OF EACH BLOCK
C
C         NBLOKS  -  INTEGER
C                     NUMBER OF DIAGONAL BLOCKS. (MUST BE GREATER
C                     THAN OR EQUAL TO 3)
C
C              A  -  DOUBLE PRECISION(M,M,NBLOKS)
C                     OUTPUT FROM TRDFAC
C
C              B  -  DOUBLE PRECISION(M,M,NBLOKS-1)
C                     OUTPUT FROM TRDFAC
C
C              C  -  DOUBLE PRECISION(M,M,NBLOKS-1)
C                     OUTPUT FROM TRDFAC
C
C         IPIVOT  -  INTEGER(M,NBLOKS)
C                     THE PIVOT INFORMATION FROM TRDFAC
C
C              Y  -  DOUBLE PRECISION(M,NBLOKS)
C                     THE RIGHT HAND SIDE VECTOR
C
C
C      ***  ON RETURN  ...
C
C              Y  -  DOUBLE PRECISION(M,NBLOKS)
C                     THE SOLUTION VECTOR
C
C***********************************************************************
C
      DOUBLE PRECISION A,B,C,Y,DP
      INTEGER M,NBLOKS,IPIVOT,I,J,K,KB,KPLS1,KMIN1
      DIMENSION A(M,M,NBLOKS),B(M,M,NBLOKS-1),C(M,M,NBLOKS-1),
     *  Y(M,NBLOKS),IPIVOT(M,NBLOKS)
C
C***********************************************************************
C
C         ***  FORWARD SOLUTION ....
C
C***********************************************************************
C
      CALL COLSOL(M,M,A,Y,IPIVOT)
      DO 30 K = 2,NBLOKS
         KMIN1 = K-1
         DO 20 J = 1,M
            DP = Y(J,KMIN1)
            DO 10 I = 1,M
               Y(I,K) = Y(I,K) - C(I,J,KMIN1)*DP
10          CONTINUE
20       CONTINUE
         CALL COLSOL(M,M,A(1,1,K),Y(1,K),IPIVOT(1,K))
30    CONTINUE
C
C***********************************************************************
C
C         ***  BACKWARD SOLUTION....
C
C***********************************************************************
C
      DO 60 KB = 1,NBLOKS-1
      K = NBLOKS-KB
      KPLS1 = K+1
         DO 50 J = 1,M
            DP = Y(J,KPLS1)
            DO 40 I = 1,M
               Y(I,K) = Y(I,K) - B(I,J,K)*DP
40          CONTINUE
50       CONTINUE
60    CONTINUE
      RETURN
      END
      SUBROUTINE COLDEC(LDA,M,A,IP,IFLAG)
C
C***********************************************************************
C
C     C O L D E C FACTORS THE MATRIX  A  USING COLUMN PIVOTING AND
C     COLUMN ELIMINATION.
C
C***********************************************************************
C
      DOUBLE PRECISION A,ZERO,ONE,ROWMAX,HOLD,AMULT,COLPIV
      INTEGER LDA,M,IP,IFLAG,I,J,MCOL,IPIVOT,IPLUS1,JPIVOT,K
      DIMENSION A(LDA,M),IP(M)
      DATA ZERO,ONE/0.D0,1.D0/
      IF (M .EQ. 1) RETURN
      DO 10 I = 1,M
         IP(I) = I
10    CONTINUE
      DO 50 I = 1,M-1
         ROWMAX = ZERO
         DO 20 J = I,M
            HOLD = DABS(A(I,IP(J)))
            IF (HOLD .GT. ROWMAX) THEN
               ROWMAX = HOLD
               MCOL = J
            ENDIF
20       CONTINUE
         IF (ONE + ROWMAX .EQ. ONE) THEN
            IFLAG = -1
            RETURN
         ENDIF
         IPIVOT = IP(MCOL)
         IF (MCOL .NE. I) THEN
            IP(MCOL) = IP(I)
            IP(I) = IPIVOT
         ENDIF
         IPLUS1 = I+1
         COLPIV =  ONE/A(I,IPIVOT)
         DO 40 J = IPLUS1,M
            JPIVOT = IP(J)
            AMULT = A(I,JPIVOT)*COLPIV
            A(I,JPIVOT) = AMULT
            DO 30 K = IPLUS1,M
               A(K,JPIVOT) = A(K,JPIVOT) - AMULT*A(K,IPIVOT)
30          CONTINUE
40       CONTINUE
50    CONTINUE
      IF (ONE + DABS(A(M,IP(M))) .EQ. ONE) IFLAG = -1
      RETURN
      END
      SUBROUTINE COLSOL(LDA,M,A,B,IP)
C
C***********************************************************************
C
C     C O L S O L SOLVES THE SYSTEM OF EQUATIONS  A*X = B USING
C     THE FACTORISATION FROM C O L D E C.
C
C***********************************************************************
C
      DOUBLE PRECISION A,B,HOLD,BJ,BNJ
      INTEGER LDA,M,IP,K,J,I,NJ,L
      DIMENSION A(LDA,M),B(M),IP(M)
      IF (M .EQ. 1) THEN
         B(1) = B(1)/A(1,1)
         RETURN
      ENDIF
      DO 20 J = 1,M-1
         B(J) = B(J)/A(J,IP(J))
         BJ = -B(J)
         DO 10 I = J+1,M
            B(I) = B(I) + A(I,IP(J))*BJ
 10      CONTINUE
 20    CONTINUE
      B(M) = B(M)/A(M,IP(M))
      DO 40 J = 1,M-1
         NJ = M-J+1
         BNJ = -B(NJ)
         DO 30 L = 1,NJ-1
            B(L) = B(L) + A(L,IP(NJ))*BNJ
 30   CONTINUE
 40   CONTINUE
      DO 50 I = 1,M
         IF (IP(I) .NE. I) THEN
            HOLD = B(IP(I))
            B(IP(I)) = B(I)
            B(I) = HOLD
            K = IP(I)
            IP(K) = K
         ENDIF
50    CONTINUE
      RETURN
      END

