PROGRAM SIMULTANEOUS_EQ_BY_LDL_DECOMPOSITION
C==================================================================
C    SOLVE SIMULTANEOUS EQUATIONS [A]{X}={B} BY  
C              [A]=[L][D][L]TRANSPOSE DECOMPSITION
C              SYMMETRY [A] NON-POSITIVE DIFINITE
C     INPUT [A] ------------------------- SQUARE MATRIX
C     OUTPUT [A] = UPPER TRIANGLE
C     DIAGONAL OF [A] IS [D]
C     [A]t-[D] = [L]TRANSPOSE
C              MARCH 25, 2011 EIJI FUKUMORI
C==================================================================
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER ( MXN=100 )
      DIMENSION A(MXN,MXN), B(MXN)
C----------------------------------------------------------------
      CALL INPUT( MXN, NNODE, A, B )
C----------------------------------------------------------------
      OPEN ( 2, FILE='SOLUTION-BY-LDL.OUT', STATUS='UNKNOWN' )
C----------------------------------------------------------------
      WRITE (2,*) 'SOLUTION OF SIMULTANEOUS EQUATIONS BY LDL'
      WRITE (2,*) 'NUMBER OF DIMENSIONS=', NNODE
      WRITE (2,*) 'MATRIX OF SIMULTANEOUS EQUATIONS'
      WRITE (2,*) 'GIVEN MATRIX: FULL MATRIX [A]'
      DO I = 1 , NNODE
      WRITE (2,*) (A(I,J), J= 1 , NNODE)
      END DO
      WRITE (2,*)
      WRITE (2,*) 'RIGHT HAND SIDE EQUATIONS {B}'
      DO J = 1 , NNODE
      WRITE(2,*) B(J)
      END DO
C----------------------------------------------------------------
C------------- SOLUTION OF  [A]{X}={B} BY [L][D][L]T DECOMP. ----
      CALL DECOMPSQ ( MXN, NNODE, A    )
      CALL SYSDCPSQ ( MXN, NNODE, A, B )
C
      WRITE (2,*)
      WRITE (2,*) 'SOLUTION OF {X} BY [L][D][L]T DECOMPOSITION'
      WRITE (2,*) 'FULL MATRIX CASE'
      DO I = 1 , NNODE
      WRITE (2,*) B(I)
      END DO
C----------------------------------------------------------------
      CLOSE (2)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE SYSDCPSQ ( MXN, NNODE, A, B )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(MXN,MXN), B(MXN)
C
C======= STEP 1: SOLUTION OF [L]{X"}={B}     [L]=LOWER [L]
C      B(1) = B(1)
      DO I = 2 , NNODE
      SUM = 0.D0
      DO J = 1, I-1
      SUM = SUM + A(I,J)*B(J)
      END DO
      B(I) = B(I) - SUM
      END DO
C
C======= STEP 2: [D]{X"'}={X"}
      DO I = 1 , NNODE
      B(I) = B(I) / A(I,I)
      END DO
C
C======= STEP 3: [L]TRANSPOSE{X}={X"'}     [L]TRANSPOSE = UPPER [L]
C      B(NNODE) = B(NNODE)
      DO I = NNODE-1, 1, -1
      SUM = 0.D0
      DO J = I+1 , NNODE
      SUM = SUM + A(I,J)*B(J)
      END DO
      B(I) = B(I) - SUM
      END DO
      RETURN
      END
C
C
      SUBROUTINE DECOMPSQ ( MXN, NNODE, A )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(MXN,MXN)
C------- COMPUTATION OF UPPER L(I,J) --------
C-------I = 1
      DO J = 2, NNODE
         A(1,J) = A(1,J)/A(1,1)
      END DO
C------ I >= 2
      DO I = 2 , NNODE
         SUM = 0.D0
            DO M = 1, I-1
              SUM = SUM + A(M,I)**2*A(M,M)
            END DO
         A(I,I) = A(I,I) - SUM
      DO J = I+1, NNODE
         SUM = 0.D0
            DO M = 1, I-1
              SUM = SUM + A(M,I)*A(M,J)*A(M,M)
            END DO
            A(I,J) = (A(I,J)-SUM)/A(I,I)
      END DO
      END DO
C -----------  MAKE LOWER [L] -----------
      DO I = 1 , NNODE
      DO J = I+1, NNODE
      A(J,I) = A(I,J)
      END DO
      END DO
      RETURN
      END
C
C
      SUBROUTINE INPUT( MXN, NNODE , A , B )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION  A(MXN,MXN), B(MXN)
      OPEN ( 1, FILE='SYSFULL.DAT', STATUS='UNKNOWN' )
      READ (1,*) NNODE
      DO I = 1 , NNODE
      READ (1,*) (A(I,J),J=1,NNODE), B(I)
      END DO
      CLOSE (1)
      RETURN
      END