      PROGRAM STAGE4B
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(950,950),CBAR(950),U(950)
      INTEGER S(950),SBEST(950),JST(950),JSL(950),Y(950),YSEL,G(950)
      OPEN(1, FILE= 'AMAT')
      OPEN(2, FILE = 'LMULT')
      OPEN(4, FILE = 'PERMUTE')
      OPEN(5, FILE = 'BBFLAG')
      READ(5,*) IFLAG
      IF(IFLAG.EQ.0) THEN
        WRITE(*,*) ' NO BRANCH AND BOUND REQUIRED '
        GO TO 889
      END IF
      READ(1,*) M,KCEN,PREF
      DO I = 1,M
        READ(2,*) U(I)
        G(I) = 0
      END DO
      EPS = 1.0D-07
      READ(1,*) ((C(I,J),J=1,M),I=1,M)
      READ(4,*) INUM
      DO I = 1,INUM
        READ(4,*) IJNK
        G(IJNK) = 1
      END DO

C ########################################################################
C OCTOBER 3, 2005: GALVAO BRANCH AND BOUND ALGORITHM WITH
C                  LAGRANGIAN RELAXATION
C ########################################################################
C
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEA = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
C
C ###################################################################
C BEGIN THE BRANCH AND BOUND ALGORITHM
C ###################################################################
C
C   ######################
C   STEP 1: INITIALIZE
C   ######################
C
      L=0
      DO K = 1,M
        S(K)=0
        JST(K) = 2  ! UNSELECTED
        JSL(K) = 0  ! ALL VARIABLES UNDECIDED
      END DO
C
      DO J = 1,M
        CBAR(J) = PREF
        DO I = 1,M
          IF(C(I,J)-U(I).LT.EPS) THEN
            CBAR(J) = CBAR(J)+C(I,J)-U(I)
          END IF
        END DO
      END DO
      YSEL = 0
      CBARMX = -0.000001   ! 99999999.0D0
      DO 213 J = 1,M
        IF(G(J).EQ.0) GO TO 213
        IF(CBAR(J).LT.CBARMX) THEN
          YSEL = YSEL + 1
          Y(YSEL) = J
          JSL(J) = 1
        END IF
 213  CONTINUE
      Z = DFLOAT(YSEL)*PREF
      DO I = 1,M
        ZIND = 9.9D10
        DO J = 1,YSEL
          JC = Y(J)
          IF(C(I,JC).LT.ZIND) ZIND = C(I,JC)
        END DO
        Z = Z+ZIND
      END DO
      Z = Z + .0001
      ZBEST = Z
      JBEST = YSEL
      WRITE(*,*) YSEL,ZBEST
      DO I = 1,YSEL
        SBEST(I)=Y(I)
      END DO
C
C
      GO TO 300
C
C   ###############################
C   STEP 2: SOLVE LAGRANGIAN DUAL
C   ###############################
C
 200  IFLAG=0
      CALL LAGRANG(M,C,IFLAG,S,L,ZBEST,U,PREF,CBAR,JSL,JST,YSEL,Y)
C      WRITE(*,211) IFLAG,JST(S(L)),L,YSEL,(s(i),i=1,l)
C 211  format(14i4)
      IF(IFLAG.EQ.1.OR.L.EQ.M) GO TO 400    ! FATHOM PARTIAL SOLUTION
      Z = DFLOAT(YSEL)*PREF
      DO I = 1,M
        ZIND = 9.9D10
        DO J = 1,YSEL
          JC = Y(J)
          IF(C(I,JC).LT.ZIND) ZIND = C(I,JC)
        END DO
        Z = Z+ZIND
      END DO
      IF(Z.LT.ZBEST) THEN
        ZBEST = Z
        JBEST = YSEL
        WRITE(*,*) YSEL,ZBEST
        DO I = 1,YSEL
          SBEST(I)=Y(I)
        END DO
      END IF
C
C   #####################################
C   STEP 3: CHOOSE A BRANCHING VARIABLE
C   #####################################
C
 300  CMAX = 999999999.
      DO 301 J = 1,M
        IF(JST(J).LE.1) GO TO 301  ! CANDIDATE FOR SELECTION
        IF(JSL(J).EQ.1) GO TO 301  ! NOT IN CURRENT SOLUTION
        IF(CBAR(J).LT.CMAX) THEN
          CMAX = CBAR(J)
          JSEL = J
        END IF
 301  CONTINUE
C
      L = L+1                  
      JST(JSEL) = 1                ! JSEL IS THE BRANCHING VARIABLE
      S(L) = JSEL
      GO TO 200
C
C   ###############################
C   STEP 4: TERMINATION CHECK
C   ###############################
C
 400  IF(L.EQ.0) GO TO 879
C
C   ###############################
C   STEP 5: BRANCH SIDEWAYS
C   ###############################
C
 500  IF(JST(S(L)).EQ.1) THEN
        JST(S(L)) = 0
c        write(*,*) ' here ', jst(s(l)
        GO TO 200
      ELSE
c        write(*,*) ' here ', jst(s(l))
        GO TO 600
      END IF
C
C   ###############################
C   STEP 6: RETRACTION
C   ###############################
C
 600  JST(S(L)) = 2
      S(L) = 0
      L = L - 1
      GO TO 400
 879  WRITE(*,*) ZBEST
C
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEB = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
      WRITE(*,834) ZBEST,TIMEB-TIMEA
      WRITE(*,835) (SBEST(K),K=1,JBEST)
  834 FORMAT(' GLOBAL OPTIMUM ',F20.5,F15.2)
  835 FORMAT(12I5)
  889 CONTINUE
      END
C
      SUBROUTINE LAGRANG(M,C,IG,S,L,ZBEST,U,PREF,CBAR,JSL,JST,YSEL,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(950,950),U(950),
     1        CBAR(950),LAMB
      INTEGER S(950),AVAIL(950),XSET(950),Y(950),XFILL(950,950),
     1        YSEL,JST(950),JSL(950)
C
      EPS = 1.0D-07
      ZDUB = -9.9D10
      IFLAG = 0
      IG = IFLAG
      DO I = 1,M
C        U(I) = 0  !UU(I) ! 0.0D0
        JSL(I) = 0
      END DO
      DO J = 1,M
        CBAR(J) = PREF
        DO I = 1,M
          IF(C(I,J)-U(I).LT.EPS) THEN
            CBAR(J) = CBAR(J)+C(I,J)-U(I)
          END IF
        END DO
      END DO
      YSEL = 0
      DO II = 1,L
        IF(JST(S(II)).EQ.1) THEN
          YSEL = YSEL + 1
          Y(YSEL)=S(II)
          JSL(S(II)) = 1
        END IF
      END DO
      CBARMX =  -0.00001   ! 99999999.0D0
      DO 446 J = 1,M
        IF(JST(J).LE.1) GO TO 446
        IF(CBAR(J).LT.CBARMX) THEN
          YSEL = YSEL + 1
          Y(YSEL) = J
          JSL(J) = 1
        END IF
 446  CONTINUE
      DO I = 1,M
        XSET(I) = 0
        DO JJ = 1,YSEL
          J = Y(JJ)
          IF(C(I,J)-U(I).LE.EPS) THEN
            XSET(I) = XSET(I) + 1
            XFILL(I,XSET(I)) = J
          END IF
        END DO
      END DO
      LAMB = 2.0D0
      MAXNI = 50    ! MAX # OF ITERATIONS WITH NO IMPROVEMENT
      NOIMP = 0     ! 
C
C  BEGIN THE INITIALIZATION PROCESS
C
      DO IJKL = 1,500
        ZDU = PREF*DFLOAT(YSEL)
        DENOM = .00001  !0.
        DO I = 1,M
          ZDU = ZDU + U(I)
          DO II = 1,XSET(I)
            J = XFILL(I,II)
            ZDU = ZDU + (C(I,J)-U(I))
          END DO
          DENOM = DENOM + (XSET(I)-1)**2
        END DO
        IF(ZDU.GT.ZBEST-EPS) THEN
          IFLAG=1           ! PRUNE IF BOUND EXCEEDS THE
          GO TO 889         ! BEST KNOWN
        END IF
        IF(ZDU.GT.ZDUB+EPS) THEN
          ZDUB = ZDU
          NOIMP = 0
        ELSE
          NOIMP = NOIMP + 1
          IF(NOIMP .EQ.MAXNI) THEN
            NOIMP = 0
            LAMB = LAMB/2.0D0
          END IF
        END IF
        T = LAMB*(ZBEST-ZDU)/DENOM
        DO I = 1,M
          U(I) = U(I) + T*(1-XSET(I))
        END DO
        DO J = 1,M
          CBAR(J) = PREF
          DO I = 1,M
            IF(C(I,J)-U(I).LT.EPS) THEN
              CBAR(J) = CBAR(J)+C(I,J)-U(I)
            END IF
          END DO
        END DO
        DO I = 1,M
          JSL(I) = 0
        END DO
        YSEL = 0
        DO II = 1,L
          IF(JST(S(II)).EQ.1) THEN
            YSEL = YSEL + 1
            Y(YSEL)=S(II)
            JSL(S(II)) = 1
          END IF
        END DO
        CBARMX = -0.00001   ! 99999999.0D0
        DO 447 J = 1,M
          IF(JST(J).LE.1) GO TO 447
          IF(CBAR(J).LT.CBARMX) THEN
            YSEL = YSEL + 1
            Y(YSEL) = J
            JSL(J) = 1
          END IF
 447    CONTINUE
        DO I = 1,M
          XSET(I) = 0
          DO JJ = 1,YSEL
            J = Y(JJ)
            IF(C(I,J)-U(I).LE.EPS) THEN
              XSET(I) = XSET(I) + 1
              XFILL(I,XSET(I)) = J
            END IF
          END DO
        END DO
      END DO
  889 CONTINUE
c      WRITE(*,765)  ZDU,ijkl
      IG = IFLAG
      WRITE(*,99) ZDU,ZBEST,L,JST(L),IFLAG  !(JST(S(I)),I=1,10)
   99 FORMAT(2F15.3,3I4)
  765 FORMAT(f20.5,i9)
      RETURN
      END
