*DECK DSOLBT SUBROUTINE DSOLBT (M, N, A, B, C, Y, IP) INTEGER M, N, IP(M,N) DOUBLE PRECISION A(M,M,N), B(M,M,N), C(M,M,N), Y(M,N) C----------------------------------------------------------------------- C Solution of block-tridiagonal linear system. C Coefficient matrix must have been previously processed by DDECBT. C M, N, A,B,C, and IP must not have been changed since call to DDECBT. C Written by A. C. Hindmarsh. C Input: C M = order of each block. C N = number of blocks in each direction of matrix. C A,B,C = M by M by N arrays containing block LU decomposition C of coefficient matrix from DDECBT. C IP = M by N integer array of pivot information from DDECBT. C Y = array of length M*N containg the right-hand side vector C (treated as an M by N array here). C Output: C Y = solution vector, of length M*N. C C External routines required: DGESL (LINPACK) and DDOT (BLAS). C----------------------------------------------------------------------- C INTEGER NM1, NM2, I, K, KB, KM1, KP1 DOUBLE PRECISION DP, DDOT NM1 = N - 1 NM2 = N - 2 C Forward solution sweep. ---------------------------------------------- CALL DGESL (A, M, M, IP, Y, 0) DO 30 K = 2,NM1 KM1 = K - 1 DO 20 I = 1,M DP = DDOT (M, C(I,1,K), M, Y(1,KM1), 1) Y(I,K) = Y(I,K) - DP 20 CONTINUE CALL DGESL (A(1,1,K), M, M, IP(1,K), Y(1,K), 0) 30 CONTINUE DO 50 I = 1,M DP = DDOT (M, C(I,1,N), M, Y(1,NM1), 1) 1 + DDOT (M, B(I,1,N), M, Y(1,NM2), 1) Y(I,N) = Y(I,N) - DP 50 CONTINUE CALL DGESL (A(1,1,N), M, M, IP(1,N), Y(1,N), 0) C Backward solution sweep. --------------------------------------------- DO 80 KB = 1,NM1 K = N - KB KP1 = K + 1 DO 70 I = 1,M DP = DDOT (M, B(I,1,K), M, Y(1,KP1), 1) Y(I,K) = Y(I,K) - DP 70 CONTINUE 80 CONTINUE DO 100 I = 1,M DP = DDOT (M, C(I,1,1), M, Y(1,3), 1) Y(I,1) = Y(I,1) - DP 100 CONTINUE RETURN C----------------------- End of Subroutine DSOLBT ---------------------- END