PROGRAM OPTIMZR4
C=======================================================================
C     EIJI FUKUMORI DECEMBER 29, 1993, 2008 JUNE 13
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      INCLUDE 'PARAM.DAT'
CCCC      PARAMETER ( ND=8, MXE=800000, MXN=800000 )
      PARAMETER ( ND2=ND*ND )
C======================================================================
C
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), ZCOORD(MXN)
      DIMENSION IREL(MXN,ND2),INOD(MXN),NEW(MXN),JNT(MXN),
     *   NJNT(MXN),TEMP(MXN), IP(ND,ND), IBNDW(MXE)
      CHARACTER INPFILE*12, OUTFILE*12, EXFILE*3
      LOGICAL YES
C=======================================================================
      WRITE (*,*) '======== WELCOME TO BANDWIDTH CRUNCHER ========'
C=======================================================================
C--------------------------- FILE MANAGEMENT ---------------------------
      INPFILE = 'NONE'
      IF ( ND .EQ. 8  ) THEN
          INPFILE = 'EIGN3D8.DAT'
          OUTFILE = 'EIGN3D8.DAT'
      END IF
      IF ( ND .EQ. 27 ) THEN
          INPFILE = 'EIGN3D27.DAT'
          OUTFILE = 'EIGN3D27.DAT'
      END IF
      IF ( INPFILE .EQ. 'NONE' ) STOP 'NO INPUTFILE CHECK ND'
C=======================================================================
      CALL PERMX ( ND, IP )
C=======================================================================
      CALL INPUT  ( ND,MXE,MXN,NE,NNODE,DELTA,MXNEIGEN,
     * NODEX, XCOORD,YCOORD, ZCOORD, INPFILE, WSPD )
C=======================================================================
      CALL BNDWDTH (ND, NODEX, MXE, NE, IBNDW ) 
C=======================================================================
      CALL OPTIMX(ND,ND2,NODEX,IREL,INOD,NEW,JNT,NJNT,XCOORD,YCOORD,
     * ZCOORD,TEMP,MXE,MXN,NE,NNODE,IMP,IP,IBNDW )
      IF ( IMP .GT. 0 ) THEN
      EXFILE = 'NEW'
      IW = 2
      OPEN ( IW, FILE=OUTFILE, STATUS='UNKNOWN' )
      WRITE (IW,*) DELTA, MXNEIGEN, WSPD
      WRITE (IW,*) NE
      DO I = 1 , NE
      WRITE (IW,*) I, (NODEX(I,J),J=1,ND)
      END DO
      WRITE (IW,*) NNODE
      DO J = 1 , NNODE
      WRITE (IW,*) J, XCOORD(J), YCOORD(J), ZCOORD(J)
      END DO
      CLOSE (IW)
      WRITE (*,*) 'OPTIMIZED. WELL DONE'
      ELSE
      WRITE (*,*) 'SOORY NO IMPROVEMENT'
      END IF
C
      IW = 7
      OPEN ( IW, FILE='BNDWDTH.FEM', STATUS='UNKNOWN' )
      CALL BNDWDTH (ND, NODEX, MXE, NE, IBNDW ) 
      WRITE (IW,*) 'ELEMENT BANDWIDTH'
      DO I = 1 , NE
      WRITE (IW,*) I , IBNDW(I)+1
      END DO
      CLOSE (IW)
      STOP 'NORMAL TERMINATION'
      END
C
C
      SUBROUTINE INPUT ( ND,MXE,MXN,NE,NNODE,DELTA,MXNEIGEN,
     * NODEX, XCOORD,YCOORD,ZCOORD, INPFILE, WSPD )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN)
      LOGICAL YES
      CHARACTER INPFILE*12
      IR = 1
      INQUIRE ( FILE=INPFILE, EXIST=YES )
      IF ( YES ) OPEN ( IR, FILE=INPFILE, STATUS='OLD' )
      IF ( .NOT. YES ) STOP' NO INPUT FILE'
C========> PARAMETER
C--- APPROPRIATE VALUES: DELTA=0.1,  MXNEIGEN=40
      READ (1,*) DELTA, MXNEIGEN, WSPD
C========> ELEMENT
      READ (1,*) NE
C------MEMORY CHECK
      IF ( NE .GT. MXE ) THEN
      WRITE (*,*) 'NE=',NE
      WRITE (*,*) 'NE > MXE'
      STOP
      END IF
      DO I = 1 , NE
      READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
      END DO
C========> COORDINATES
      READ (1,*) NNODE
C------MEMORY CHECK
      IF ( NNODE .GT. MXN ) THEN
      WRITE (*,*) 'NNODE=',NNODE
      WRITE (*,*) 'NNODE > MXN'
      STOP
      END IF
      DO I = 1 , NNODE
      READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE), ZCOORD(NODE)
      END DO
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE OPTIMX ( ND,ND2,NODEX,IREL,INOD, NEW, JNT, NJNT,XCOORD,
     *  YCOORD,ZCOORD,TEMP,MXE,MXN,NE,NNODE,IMP,IP,IBNDW )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION NODEX(MXE,ND), IREL(MXN,ND2), INOD(MXN), 
     & NEW(MXN),JNT(MXN),NJNT(MXN),XCOORD(MXN),YCOORD(MXN),ZCOORD(MXN),
     *  IP(ND,ND), TEMP(MXN),JTEMP(ND2), IBNDW(MXE)
C-------------------------------------------------------------------------
      IW = 4
      OPEN ( IW, FILE='OPTIMIZATION-STATUS.FEM', STATUS='UNKNOWN' )
C-------------------------------------------------------------------------
      WRITE (*,*) 'OPTIMIZATION STARTED'
      CALL SETUP( ND, ND2,NODEX,IREL,INOD,MXE,MXN,NE,NNODE,IBWOR,IP,
     *            IBNDW )
C-------------------------------------------------------------------------
      IMP = 0
      CALL OPTIM(ND2,IREL,INOD,JNT,NJNT,NEW,MXN,NNODE,IBWOR,MAXDIF,IMP) 
      IBW = MAXDIF
      DO I = 1, NNODE 
      DO J = 1, INOD(I)
      JTEMP(J) = IREL(I,INOD(I)-J+1)
      END DO
      DO K = 1, INOD(I) 
      IREL(I,K) = JTEMP(K) 
      END DO
      END DO
      WRITE (*,*) 'SECOND TRIAL'
      CALL OPTIM(ND2,IREL,INOD,JNT,NJNT,NEW,MXN,NNODE,IBW,MAXDIF,IMP) 
      IF ( IMP .LE. 0 ) THEN
      IBOLD = IBWOR + 1 
      WRITE(IW,*) '===== NODE NUMBERING CANNOT BE IMPROVED.====='
      WRITE(IW,*) '===== HALF-BANDWIDTH REMAINS AT ', IBOLD,' ====='
      RETURN
      END IF
      IBOLD = IBWOR + 1 
      MAXBWH = MAXDIF + 1 
      WRITE(IW,*) 'ORIGINAL HALF-BANDWIDTH =', IBOLD
      WRITE(IW,*) 'REDUCED  HALF-BANDWIDTH =', MAXBWH 
      DO I = 1, NE 
      DO J = 1, ND 
      NODEX(I,J) = NEW( NODEX(I,J) ) 
      END DO
      END DO
      DO I = 1, NNODE
      TEMP(I) = XCOORD(I) 
      END DO
      DO I = 1, NNODE
      XCOORD(NEW(I)) = TEMP(I)
      TEMP(I) = YCOORD(I) 
      END DO
      DO I = 1, NNODE
      YCOORD(NEW(I)) = TEMP(I)
      TEMP(I) = ZCOORD(I) 
      END DO
      DO I = 1, NNODE
      ZCOORD(NEW(I)) = TEMP(I) 
      END DO
C
      WRITE(IW,*) 'NEW OPTIMIZED NODE NUMBERING:'
      WRITE(IW,*) '(OLD NODE NO.--> NEW NODE NO.)'
      WRITE(IW,'(5(1H(,I7,4H -->,I7,4H )  ))') (I,NEW(I),I=1,NNODE) 
      RETURN
      END
C 
C 
      SUBROUTINE OPTIM ( ND2,IREL,INOD,JNT,NJNT,NEW,MXN,NNODE,IDIFF,
     &  MAXDIF, IMP )
      DIMENSION IREL(MXN,ND2),INOD(MXN),JNT(MXN),NJNT(MXN), NEW(MXN) 
      MAXDIF = IDIFF
      DO 60 IX = 1, NNODE
       IF ( IX/1000*1000 .EQ. IX ) WRITE(*,*)'OPTIMIZATION AT NODE =',IX
       DO J = 1, NNODE 
       JNT(J) = 0
       NJNT(J) = 0 
       END DO
       MAX = 0 
       I = 1 
       NJNT(1) = IX
       JNT(IX) = 1 
       K = 1
       DO
         DO JJ = 1, INOD(NJNT(I))
           KX = IREL(NJNT(I),JJ) 
           IF ( JNT(KX) .LE. 0 ) THEN
             K = K + 1 
             NJNT(K) = KX
             JNT(KX) = K 
             NDIFF = IABS ( I - K )
             IF ( NDIFF .GE. MAXDIF ) GO TO 60 
              MAX = MAX0 ( NDIFF, MAX ) 
           END IF
         END DO
         IF ( K .EQ. NNODE ) EXIT 
         I = I + 1 
       END DO
       MAXDIF = MAX
       DO J = 1, NNODE 
       NEW(J) = JNT(J) 
       END DO
       IMP = IMP + 1 
   60 CONTINUE
      RETURN
      END
C 
C 
      SUBROUTINE SETUP (ND,ND2,NODEX,IREL,INOD,MXE,MXN,NE,NNODE,IBWOR,
     *                  IP,  IBNDW ) 
      DIMENSION NODEX(MXE,ND),IREL(MXN,ND2),INOD(MXN),IP(ND,ND),
     *          IBNDW(MXE)
C
      IBWOR = 0
      DO NODE = 1, NNODE
       KOUNT = 0
       DO IEL = 1, NE
         DO K = 1, ND
           IF ( NODE .EQ. NODEX(IEL,IP(K,1)) ) THEN
           DO J = 2, ND
           CALL SQUENCE ( ND2,MXN,NODE,NODEX(IEL,IP(K,J)),KOUNT,IREL )
           END DO
           END IF
         END DO
         IBWOR = MAX0 ( IBNDW(IEL), IBWOR ) 
       END DO
       INOD(NODE) = KOUNT
       IF ( NODE/1000*1000 .EQ. NODE ) WRITE(*,*)'SETUP AT NODE =',NODE
      END DO
      RETURN
      END
C
C
      SUBROUTINE SQUENCE ( ND2, MXN, NODE, IPNODE, KOUNT, IREL )
      DIMENSION IREL(MXN,ND2)
      IF ( KOUNT .GT. 0 ) THEN
      DO L = 1, KOUNT
      IF ( IREL(NODE,L) .EQ. IPNODE ) RETURN
      END DO
      END IF
      KOUNT = KOUNT + 1 
      IREL(NODE,KOUNT) = IPNODE
      RETURN
      END
C
C
      SUBROUTINE PERMX ( ND, IP )
      DIMENSION IP(ND,ND)
      DO I = 1 , ND
      K = I
      IP(I,1) = I
      DO J = 2 , ND
      K = K + 1
      IF ( K .GT. ND ) K = K - ND
      IP(I,J) = K
      END DO
      END DO
      RETURN
      END
C 
C 
      SUBROUTINE BNDWDTH (ND, NODEX, MXE, NE, IBNDW ) 
      DIMENSION NODEX(MXE,ND), IBNDW(MXE)
C
      DO IEL = 1, NE
      IBNDW(IEL) = 0
      DO  J = 1 , ND-1
      DO  K = J+1 , ND
      IBNDW(IEL) = MAX0(IBNDW(IEL),IABS(NODEX(IEL,J)-NODEX(IEL,K)))
      END DO
      END DO
      END DO
      RETURN
      END