C SUBROUTINE BLKTRI(A,B,C,F,IF,I1,I2,KF,K1,K2,BETA,GAMMA,Y,X) implicit none C **** C **** This procedure solves (I2-I1+1) tridiagonal block matrix C **** systems in which all blocks are 2 x 2 matrices. C **** C **** Each system may be written: C **** C **** A(K) * X(K-1) + B(K) * X(K) + Z(K) * X(K+1) = F(K) C **** C **** where: C **** C **** K = K1,K2,1 C **** C **** A(K), B(K), C(K) are given (2 x 2) matrices. C **** C **** The F(k) are given two componente vectors. C **** C **** The system is to be solved for the two component C **** vectors, X(K). C **** C **** A(K1) = C(K2) = 0. C **** C **** BETA(K), GAMMA(K), (K = K1,K2,1), are work space for C **** (2 x 2) matrices. C **** C **** Y(K), (K = K1,K2,1), is work space for two component C **** vectors. C **** C **** Algorithm: (See Isaacson and Keller p55) C **** C **** Forward sweep from K = K1 to K = K2: C **** C **** BETA(K1) = B(K1)**(-1) C **** C **** Y(K1) = BETA(K1)*F(K1) C **** C **** GAMMA(K) = BETA(K)*C(K), K = K1,(K2-1),1 C **** C **** BETA(K) = (B(K) - A(K)*GAMMA(K-1))**(-1), C **** K = K1+1,K2,1 C **** C **** Y(K) = BETA(K)*(F(K) - A(K)*Y(K-1)), C **** K = K1+1,K2,1 C **** C **** Backward sweep, K = K2,K1,-1 C **** C **** X(K2) = Y(K2) C **** C **** X(K) = Y(K) - GAMMA(K)*X(K+1), K = K2-1,K1,-1 C **** C **** Dimension statements: C **** C **** Block matrices are dimensioned thus: C **** C **** MATRIX(2,2,IF,KF) C **** C **** Two component vectors are similarly treated: C **** C **** VECTOR(2,IF,KF) C **** C **** Our block matrix scheme spans the range, (K = K1,K2,1), C **** where (1 .LE. K1 .LT. K2 .LE. KF) C **** C **** Similarly, we are solving (I1-I2+1) systems C **** simultaneously as the index, I, spans the range, C **** (I = I1,I2,1), where (1 .LE. I1 .LT. I2 .LE. IF) C **** ! ! Args: real,intent(in) :: a,b,c,f real,intent(out) :: x,y,beta,gamma integer,intent(in) :: if,i1,i2,kf,k1,k2 ! ! Local: integer :: i,k ! DIMENSION A(2,2,IF,KF), B(2,2,IF,KF), C(2,2,IF,KF), F(2,IF,KF), 1 BETA(2,2,IF,KF), GAMMA(2,2,IF,KF), Y(2,IF,KF), X(2,IF,KF) C **** C **** Lower boundary at K = K1 C **** DO I = I1,I2 C **** C **** Y(1,I,K1) = determinant(B(K1)) C **** Y(1,I,K1) = B(1,1,I,K1)*B(2,2,I,K1) - B(1,2,I,K1)*B(2,1,I,K1) C **** C **** BETA(K1) = B(K1)**(-1) C **** BETA(1,1,I,K1) = B(2,2,I,K1)/Y(1,I,K1) BETA(1,2,I,K1) = -B(1,2,I,K1)/Y(1,I,K1) BETA(2,1,I,K1) = -B(2,1,I,K1)/Y(1,I,K1) BETA(2,2,I,K1) = B(1,1,I,K1)/Y(1,I,K1) C **** C **** Y(K1) = BETA(K1)*F(K1) C **** Y(1,I,K1) = BETA(1,1,I,K1)*F(1,I,K1) + BETA(1,2,I,K1)*F(2,I,K1) Y(2,I,K1) = BETA(2,1,I,K1)*F(1,I,K1) + BETA(2,2,I,K1)*F(2,I,K1) ENDDO C **** C **** Now deal with levels (K1+1),K2,1 C **** DO K = K1+1,K2 DO I = I1,I2 C **** C **** GAMMA(K-1) = BETA(K-1)*C(K-1) C **** GAMMA(1,1,I,K-1) = BETA(1,1,I,K-1)*C(1,1,I,K-1) + 1 BETA(1,2,I,K-1)*C(2,1,I,K-1) GAMMA(1,2,I,K-1) = BETA(1,1,I,K-1)*C(1,2,I,K-1) + 1 BETA(1,2,I,K-1)*C(2,2,I,K-1) GAMMA(2,1,I,K-1) = BETA(2,1,I,K-1)*C(1,1,I,K-1) + 1 BETA(2,2,I,K-1)*C(2,1,I,K-1) GAMMA(2,2,I,K-1) = BETA(2,1,I,K-1)*C(1,2,I,K-1) + 1 BETA(2,2,I,K-1)*C(2,2,I,K-1) C **** C **** GAMMA(K) = B(K) - A(K)*GAMMA(K-1) C **** GAMMA(1,1,I,K) = B(1,1,I,K) - A(1,1,I,K)*GAMMA(1,1,I,K-1) - 1 A(1,2,I,K)*GAMMA(2,1,I,K-1) GAMMA(1,2,I,K) = B(1,2,I,K) - A(1,1,I,K)*GAMMA(1,2,I,K-1) - 1 A(1,2,I,K)*GAMMA(2,2,I,K-1) GAMMA(2,1,I,K) = B(2,1,I,K) - A(2,1,I,K)*GAMMA(1,1,I,K-1) - 1 A(2,2,I,K)*GAMMA(2,1,I,K-1) GAMMA(2,2,I,K) = B(2,2,I,K) - A(2,1,I,K)*GAMMA(1,2,I,K-1) - 1 A(2,2,I,K)*GAMMA(2,2,I,K-1) C **** C **** Y(1,I,K) = determinant(GAMMA(K)) C **** Y(1,I,K) = GAMMA(1,1,I,K)*GAMMA(2,2,I,K) - 1 GAMMA(1,2,I,K)*GAMMA(2,1,I,K) C **** C **** BETA(K) = GAMMA(K)**(-1) C **** BETA(1,1,I,K) = GAMMA(2,2,I,K)/Y(1,I,K) BETA(1,2,I,K) = -GAMMA(1,2,I,K)/Y(1,I,K) BETA(2,1,I,K) = -GAMMA(2,1,I,K)/Y(1,I,K) BETA(2,2,I,K) = GAMMA(1,1,I,K)/Y(1,I,K) C **** C **** X(K) = F(K) - A(K)*Y(K-1) C **** X(1,I,K) = F(1,I,K) - A(1,1,I,K)*Y(1,I,K-1) - 1 A(1,2,I,K)*Y(2,I,K-1) X(2,I,K) = F(2,I,K) - A(2,1,I,K)*Y(1,I,K-1) - 1 A(2,2,I,K)*Y(2,I,K-1) C **** C **** Y(K) = BETA(K)*X(K) C **** Y(1,I,K) = BETA(1,1,I,K)*X(1,I,K) + BETA(1,2,I,K)*X(2,I,K) Y(2,I,K) = BETA(2,1,I,K)*X(1,I,K) + BETA(2,2,I,K)*X(2,I,K) ENDDO ENDDO C **** C **** Backward sweep to determine final solution, X(K) for C **** K = K2,K1,-1 C **** C **** X(K2) = Y(K2) C **** DO I = I1,I2 X(1,I,K2) = Y(1,I,K2) X(2,I,K2) = Y(2,I,K2) ENDDO C **** C **** X(K) = Y(K) - GAMMA(K)*X(K+1) C **** DO K = K2-1,K1,-1 DO I = I1,I2 X(1,I,K) = Y(1,I,K) - GAMMA(1,1,I,K)*X(1,I,K+1) - 1 GAMMA(1,2,I,K)*X(2,I,K+1) X(2,I,K) = Y(2,I,K) - GAMMA(2,1,I,K)*X(1,I,K+1) - 1 GAMMA(2,2,I,K)*X(2,I,K+1) ENDDO ENDDO RETURN END