PROGRAM SIMULTANEOUS
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]-[D] = [L]TRANSPOSE
C              AUGUST 01, 2008 EIJI FUKUMORI
C==================================================================
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER ( MXN=100,MXW=20, NEMX=MXN*(3*MXN-1)/2 )
      DIMENSION A(MXN,MXN), IA(NEMX), B(MXN), ALDL(MXW,MXN),
     * VRTCL(MXN,MXN), BACK(MXW,MXN)
C
      NNODE = 8
      IBAND = 5
      CALL MTXMAKE ( MXN,MXW,NEMX,NNODE,IBAND,IA,A,B,ALDL,
     * BACK,VRTCL )
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE MTXMAKE ( MXN,MXW,NEMX,NNODE,IBAND,IA,A,B,ALDL,
     * BACK,VRTCL )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(MXN,MXN), IA(NEMX), B(MXN), ALDL(MXW,MXN)
      DIMENSION VRTCL(MXN,MXN),  BACK(MXW,MXN)
      NELEMNT = NNODE*(3*NNODE-1)/2 + NNODE
      OPEN ( 1, FILE='PI.DAT', STATUS='OLD' )
      READ (1,*)
      READ (1,*)
      READ(1,'(4(5I2,1X),5I2)') ( IA(I), I= 1, NELEMNT )
C
C-------- MATRIX [A]
      M = 0
      DO I = 1 , NNODE
      MAXJ = IBAND
      IF (NNODE-I+1 .LT. IBAND ) MAXJ = NNODE-I+1
      DO J = I , I+MAXJ-1
      M = M + 1
      A(I,J) = IA(M)
      A(J,I) = A(I,J)
      END DO
      END DO
C
C------- VECTOR {B}
      DO I = 1 , NNODE
      M = M + 1
      B(I) = IA(M)
      END DO
C==================================================================
C------- CREATING FULL MATRIX DATA FILE
      OPEN ( 2, FILE='SYSFULL.DAT', STATUS='UNKNOWN' )
      WRITE (2,*) NNODE
      DO I = 1 , NNODE
      WRITE (2,*) (A(I,J),J=1,NNODE), B(I)
      END DO
      CLOSE (2)
C==================================================================
C--------- CREATING CRUSH MATRIX (VERTICAL)---------
      DO I = 1 , NNODE
      DO J = 1 , IBAND
      VRTCL(I,J) = 0.D0
      END DO
      END DO
      DO I = 1 , NNODE
      DO J = I , NNODE
      JJ = J - I + 1
      VRTCL(I,JJ) = A(I,J)
      END DO
      END DO
      OPEN ( 2, FILE='MTRIXLDLVT.DAT', STATUS='UNKNOWN' )
      WRITE (2,*) NNODE, IBAND
      DO I = 1 , NNODE
      WRITE (2,'(20G14.5)') ( VRTCL(I,J), J=1,IBAND )
      END DO
      WRITE (2,'(20G14.5)') ( B(I), I=1,NNODE )
      CLOSE (2)
      RETURN
      END