PROGRAM SIMLUTANEOUS 
C==================================================================
C   SOLUTION OF AX=B BY LDLT DECOMPOSITION 
C              [A]=[L][D][L]T
C     INPUT [A]
C     OUTPUT [A] = UPPER TRIANGLE
C     DIAGONAL OF [A] IS [D]
C     [A]-[D] = [L]TRANSPOSE
C  SUBROUTINE LLDECOMP3 HAND MADE
C
CXXXXX
CXXXXX
CXXXX
CXXX
CXX
CX
C              AUGUST 01, 2008 EIJI FUKUMORI
C==================================================================
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER ( MXN=200, MXW=50 )
      DIMENSION A(MXN,MXW),  B(MXN)
C-----------------------------------------------------------------
      CALL INPUT ( MXW, MXN, NNODE, IBAND, A , B )
C-----------------------------------------------------------------
      OPEN ( 21, FILE='SOLUTION-CRUSHVT.SIM', STATUS='UNKNOWN' )
      WRITE (21,*) 'SOLUTION OF SIMULTANEOUS EQUATIONS BY LDL'
      WRITE (21,*) 'NUMBER OF DIMENSIONS=', NNODE
      WRITE (21,*) 'MATRIX OF SIMULTANEOUS EQUATIONS'
      WRITE (21,*) 'GIVEN MATRIX: UPPER HALF OF [A]'
      DO I = 1,NNODE
      WRITE (21,*)  ( A(I,J), J = 1 ,IBAND  )
      END DO
      WRITE (21,*) 'RIGHT HAND SIDE OF EQUATION {B}'
      DO I = 1 , NNODE
      WRITE (21,*)  B(I)
      END DO
C------- COMPUTATION OF RLT(I,J)
C-----------------------------------------------------------------
      CALL LLDECOMPVT ( A, NNODE, IBAND, MXW, MXN )
C-----------------------------------------------------------------
      CALL SYSLDLVT ( MXN, MXW, NNODE, IBAND, A, B )
C-----------------------------------------------------------------
C
      WRITE (21,*) 'SOLUTION OF {X} BY [L][D][L]T DECOMPOSITION'
      WRITE (21,*) 'HALF BAND MATRIX CASE'
      DO I = 1 , NNODE
      WRITE (21,*)  B(I)
      END DO
C
      CLOSE (21)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE LLDECOMPVT ( A, NNODE, IBAND, MXW, MXN )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(MXN,MXW)
C-------I = 1
      DO J = 2, IBAND
         A(1,J) = A(1,J)/A(1,1)
      END DO
C------ I >= 2
      DO I = 2 , NNODE
         SUM = 0.D0
            DO M = 2, MIN0(I,IBAND)
              II = I - M + 1
              SUM = SUM + A(II,M)**2*A(II,1)
            END DO
         A(I,1) = A(I,1) - SUM
       DO J = 2, MIN0(IBAND,NNODE-I+1)
         SUM = 0.D0
            DO M = 2, MIN0( I,IBAND, IBAND-J+1 )
              II = I - M + 1
              L = J + M - 1
              SUM = SUM + A(II,1)*A(II,M)*A(II,L)
            END DO
            A(I,J) = (A(I,J)-SUM)/A(I,1)
       END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE INPUT ( MXW, MXN, NNODE, IBAND, A, B )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(MXN,MXW), B(MXN)
      OPEN ( 2, FILE='MTRIXLDLVT.DAT', STATUS='UNKNOWN' )
      READ (2,*) NNODE, IBAND
      DO I = 1,NNODE
      READ (2,*) ( A(I,J), J=1,IBAND )
      END DO
      READ (2,*) ( B(I), I= 1, NNODE )
      CLOSE (2)
      RETURN
      END
C
C
      SUBROUTINE SYSLDLVT ( MXN, MXW, NNODE, NBWDTH, GSMTX, V1 )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION GSMTX(MXN,MXW),  V1(MXN)
C------- SOLUTION OF [A]{X}={B} --------
C----- SEQUENCE 1:  [L]{X'}={B}
      DO I=2,NNODE
       DO J = 2, MIN0(I,NBWDTH)
        V1(I) = V1(I) - GSMTX(I-J+1,J)*V1(I-J+1)
       END DO
      END DO
C------ SEQUENCE 2 [D]{X"'}={X"}
      DO I = 1 , NNODE
       V1(I) = V1(I) / GSMTX(I,1)
      END DO
C------ SEQUENCE 3  [L]TRANSPOSE{X'}={X"'}   NOTE THAT L(I,I) = 1
      DO I = NNODE-1 , 1 , -1
       DO J = 2 , MIN0(NNODE-I+1,NBWDTH)
          V1(I) = V1(I) - GSMTX(I,J)*V1(I+J-1)
       END DO
      END DO
      RETURN
      END