share/odepack/fortran/dsolbt.f

   1 *DECK DSOLBT
   2       SUBROUTINE DSOLBT (M, N, A, B, C, Y, IP)
   3       INTEGER M, N, IP(M,N)
   4       DOUBLE PRECISION A(M,M,N), B(M,M,N), C(M,M,N), Y(M,N)
   5 C-----------------------------------------------------------------------
   6 C Solution of block-tridiagonal linear system.
   7 C Coefficient matrix must have been previously processed by DDECBT.
   8 C M, N, A,B,C, and IP  must not have been changed since call to DDECBT.
   9 C Written by A. C. Hindmarsh.
  10 C Input:
  11 C     M = order of each block.
  12 C     N = number of blocks in each direction of matrix.
  13 C A,B,C = M by M by N arrays containing block LU decomposition
  14 C         of coefficient matrix from DDECBT.
  15 C    IP = M by N integer array of pivot information from DDECBT.
  16 C     Y = array of length M*N containg the right-hand side vector
  17 C         (treated as an M by N array here).
  18 C Output:
  19 C     Y = solution vector, of length M*N.
  20 C
  21 C External routines required: DGESL (LINPACK) and DDOT (BLAS).
  22 C-----------------------------------------------------------------------
  23 C
  24       INTEGER NM1, NM2, I, K, KB, KM1, KP1
  25       DOUBLE PRECISION DP, DDOT
  26       NM1 = N - 1
  27       NM2 = N - 2
  28 C Forward solution sweep. ----------------------------------------------
  29       CALL DGESL (A, M, M, IP, Y, 0)
  30       DO 30 K = 2,NM1
  31         KM1 = K - 1
  32         DO 20 I = 1,M
  33           DP = DDOT (M, C(I,1,K), M, Y(1,KM1), 1)
  34           Y(I,K) = Y(I,K) - DP
  35  20       CONTINUE
  36         CALL DGESL (A(1,1,K), M, M, IP(1,K), Y(1,K), 0)
  37  30     CONTINUE
  38       DO 50 I = 1,M
  39         DP = DDOT (M, C(I,1,N), M, Y(1,NM1), 1)
  40      1     + DDOT (M, B(I,1,N), M, Y(1,NM2), 1)
  41         Y(I,N) = Y(I,N) - DP
  42  50     CONTINUE
  43       CALL DGESL (A(1,1,N), M, M, IP(1,N), Y(1,N), 0)
  44 C Backward solution sweep. ---------------------------------------------
  45       DO 80 KB = 1,NM1
  46         K = N - KB
  47         KP1 = K + 1
  48         DO 70 I = 1,M
  49           DP = DDOT (M, B(I,1,K), M, Y(1,KP1), 1)
  50           Y(I,K) = Y(I,K) - DP
  51  70       CONTINUE
  52  80     CONTINUE
  53       DO 100 I = 1,M
  54         DP = DDOT (M, C(I,1,1), M, Y(1,3), 1)
  55         Y(I,1) = Y(I,1) - DP
  56  100    CONTINUE
  57       RETURN
  58 C----------------------- End of Subroutine DSOLBT ----------------------
  59       END