share/odepack/fortran/ddecbt.f

   1 *DECK DDECBT
   2       SUBROUTINE DDECBT (M, N, A, B, C, IP, IER)
   3       INTEGER M, N, IP(M,N), IER
   4       DOUBLE PRECISION A(M,M,N), B(M,M,N), C(M,M,N)
   5 C-----------------------------------------------------------------------
   6 C Block-tridiagonal matrix decomposition routine.
   7 C Written by A. C. Hindmarsh.
   8 C Latest revision:  November 10, 1983 (ACH)
   9 C Reference:  UCID-30150
  10 C             Solution of Block-Tridiagonal Systems of Linear
  11 C             Algebraic Equations
  12 C             A.C. Hindmarsh
  13 C             February 1977
  14 C The input matrix contains three blocks of elements in each block-row,
  15 C including blocks in the (1,3) and (N,N-2) block positions.
  16 C DDECBT uses block Gauss elimination and Subroutines DGEFA and DGESL
  17 C for solution of blocks.  Partial pivoting is done within
  18 C block-rows only.
  19 C
  20 C Note: this version uses LINPACK routines DGEFA/DGESL instead of
  21 C of dec/sol for solution of blocks, and it uses the BLAS routine DDOT
  22 C for dot product calculations.
  23 C
  24 C Input:
  25 C     M = order of each block.
  26 C     N = number of blocks in each direction of the matrix.
  27 C         N must be 4 or more.  The complete matrix has order M*N.
  28 C     A = M by M by N array containing diagonal blocks.
  29 C         A(i,j,k) contains the (i,j) element of the k-th block.
  30 C     B = M by M by N array containing the super-diagonal blocks
  31 C         (in B(*,*,k) for k = 1,...,N-1) and the block in the (N,N-2)
  32 C         block position (in B(*,*,N)).
  33 C     C = M by M by N array containing the subdiagonal blocks
  34 C         (in C(*,*,k) for k = 2,3,...,N) and the block in the
  35 C         (1,3) block position (in C(*,*,1)).
  36 C    IP = integer array of length M*N for working storage.
  37 C Output:
  38 C A,B,C = M by M by N arrays containing the block-LU decomposition
  39 C         of the input matrix.
  40 C    IP = M by N array of pivot information.  IP(*,k) contains
  41 C         information for the k-th digonal block.
  42 C   IER = 0  if no trouble occurred, or
  43 C       = -1 if the input value of M or N was illegal, or
  44 C       = k  if a singular matrix was found in the k-th diagonal block.
  45 C Use DSOLBT to solve the associated linear system.
  46 C
  47 C External routines required: DGEFA and DGESL (from LINPACK) and
  48 C DDOT (from the BLAS, or Basic Linear Algebra package).
  49 C-----------------------------------------------------------------------
  50       INTEGER NM1, NM2, KM1, I, J, K
  51       DOUBLE PRECISION DP, DDOT
  52       IF (M .LT. 1 .OR. N .LT. 4) GO TO 210
  53       NM1 = N - 1
  54       NM2 = N - 2
  55 C Process the first block-row. -----------------------------------------
  56       CALL DGEFA (A, M, M, IP, IER)
  57       K = 1
  58       IF (IER .NE. 0) GO TO 200
  59       DO 10 J = 1,M
  60         CALL DGESL (A, M, M, IP, B(1,J,1), 0)
  61         CALL DGESL (A, M, M, IP, C(1,J,1), 0)
  62  10     CONTINUE
  63 C Adjust B(*,*,2). -----------------------------------------------------
  64       DO 40 J = 1,M
  65         DO 30 I = 1,M
  66           DP = DDOT (M, C(I,1,2), M, C(1,J,1), 1)
  67           B(I,J,2) = B(I,J,2) - DP
  68  30       CONTINUE
  69  40     CONTINUE
  70 C Main loop.  Process block-rows 2 to N-1. -----------------------------
  71       DO 100 K = 2,NM1
  72         KM1 = K - 1
  73         DO 70 J = 1,M
  74           DO 60 I = 1,M
  75             DP = DDOT (M, C(I,1,K), M, B(1,J,KM1), 1)
  76             A(I,J,K) = A(I,J,K) - DP
  77  60         CONTINUE
  78  70       CONTINUE
  79         CALL DGEFA (A(1,1,K), M, M, IP(1,K), IER)
  80         IF (IER .NE. 0) GO TO 200
  81         DO 80 J = 1,M
  82  80       CALL DGESL (A(1,1,K), M, M, IP(1,K), B(1,J,K), 0)
  83  100    CONTINUE
  84 C Process last block-row and return. -----------------------------------
  85       DO 130 J = 1,M
  86         DO 120 I = 1,M
  87           DP = DDOT (M, B(I,1,N), M, B(1,J,NM2), 1)
  88           C(I,J,N) = C(I,J,N) - DP
  89  120      CONTINUE
  90  130    CONTINUE
  91       DO 160 J = 1,M
  92         DO 150 I = 1,M
  93           DP = DDOT (M, C(I,1,N), M, B(1,J,NM1), 1)
  94           A(I,J,N) = A(I,J,N) - DP
  95  150      CONTINUE
  96  160    CONTINUE
  97       CALL DGEFA (A(1,1,N), M, M, IP(1,N), IER)
  98       K = N
  99       IF (IER .NE. 0) GO TO 200
 100       RETURN
 101 C Error returns. -------------------------------------------------------
 102  200  IER = K
 103       RETURN
 104  210  IER = -1
 105       RETURN
 106 C----------------------- End of Subroutine DDECBT ----------------------
 107       END