      PROGRAM STAGE4A
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(1500,1500),U(1500),junk,
     1                 CBAR(1500),LAMB,UU(1500),
     1                 MEDIAN
      INTEGER Y(1500),XSET(1500),DENOM,V,XFILL(1500,1500),YSEL,G(1500)
      OPEN(1, FILE='AMAT')
      OPEN(2, FILE='LMULT')
      OPEN(3, FILE='HEUR')
      OPEN(5, FILE='BBFLAG')
      OPEN(6, FILE='EXEMPLARS')
      OPEN(7, FILE='PERMUTE')
      READ(1,*) M,K,MEDIAN
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEC = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
      READ(5,*) IFLAG
      IF(IFLAG.EQ.0) THEN
        WRITE(*,*) ' NO BRANCH AND BOUND REQUIRED '
        GO TO 889
      END IF
      READ(3,*) zbest         ! ZBEST = UPPER BOUND FROM TEITZ & BART
      DO J = 1,M              ! UU(J) = INITIAL LAGR. MULTIPLIERS
C        READ(2,*) UU(J)
        UU(J) = 0
        G(J) = 0
      END DO
      EPS = 1.0D-07           ! EPSILON FOR PRECISION
      READ(1,*) ((C(I,J),J=1,M),I=1,M)  ! A = VARIABLE MEAS. MATRIX
      READ(7,*) INUM
      DO J = 1,INUM
        READ(7,*) IJNK
        G(IJNK) = 1
      END DO
C
      DO 1100 IJ = 1,5          ! RUN 5 RESTARTS OF LAGR. RELAXATION
        ZDUB = -9.9D10          ! ZDUB = BEST FOUND LB FOR RESTART
        DO I = 1,M
          U(I) = UU(I)          ! REINSTALL BEST FOUND LARG. MULTIPLIERS
        END DO
        YSEL = 0
        DO 230 J = 1,M
          IF(G(J).EQ.0) GO TO 230
          CBAR(J) = MEDIAN      ! CBAR(J) WILL COLLECT SUM OF C(I,J)-U(J)
          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
          IF(CBAR(J).LT.0) THEN
            YSEL = YSEL + 1
            Y(YSEL)=J
          END IF
  230   CONTINUE
        DO I = 1,M                       
          XSET(I) = 0
          DO JJ = 1,YSEL
            J = Y(JJ)
            IF(C(I,J)-U(I).LT.EPS) THEN  
              XSET(I) = XSET(I) + 1
              XFILL(I,XSET(I)) = J
            END IF
          END DO
        END DO
        LAMB = 2.0D0  ! INITIAL SETTING OF 2 FOR LAMBDA PARAMETER
        MAXNI = 200   ! MAX # OF ITERATIONS WITH NO IMPROVEMENT
        NOIMP = 0     ! CURRENT # OF ITERATIONS WITH NO IMPROVEMENT
C
C  BEGIN THE INITIALIZATION PROCESS
C
        DO IJKL = 1,2000        ! RUN 2000 TRIALS PER RESTART
          ZDU = MEDIAN*DFLOAT(YSEL)  ! ZDU = THE CURRENT LB
          DENOM = 0             ! DENOM = SQRD DEVIATION FROM CONST.
          DO I = 1,M            !         SATISFACTION
            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(DENOM.EQ.0.OR.ZDU.GT.ZBEST-.00001) THEN
            IFLAG=0
            GO TO 80
          END IF
          IF(ZDU.GT.ZDUB+EPS) THEN    ! THIS IF BLOCK CHECKS FOR
            ZDUB = ZDU                ! IMPROVEMENT IN THE BOUND
            NOIMP = 0                 ! INCREMENTS NOIMP ACCORDINGLY
          ELSE
            NOIMP = NOIMP + 1
            IF(NOIMP .EQ.MAXNI) THEN  ! LAMBDA IN HALVED WHEN NOIMP
              NOIMP = 0               ! REACHES THE MAXIMUM
              LAMB = LAMB/2.0D0
            END IF
          END IF
          T = LAMB*(ZBEST-ZDU)/DFLOAT(DENOM)  ! COMPUTE T PARAMETER
          DO I = 1,M
            U(I) = U(I) + T*(1-XSET(I)) ! UPDATE LAGR. MULTIPLIERS
          END DO
          YSEL = 0
          DO 240 J = 1,M
            IF(G(J).EQ.0) GO TO 240
            CBAR(J) = MEDIAN          ! CBAR(J) WILL COLLECT SUM OF C(I,J)-U(J)
            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
            IF(CBAR(J).LT.0) THEN
              YSEL = YSEL + 1
              Y(YSEL)=J
            END IF
 240      CONTINUE
          DO I = 1,M                       
            XSET(I) = 0
            DO JJ = 1,YSEL
              J = Y(JJ)
              IF(C(I,J)-U(I).LT.EPS) THEN    ! C(I,J) - U(I) .LT. 0)
                XSET(I) = XSET(I) + 1
                XFILL(I,XSET(I)) = J
              END IF
            END DO
          END DO
        END DO
        DO I = 1,M
          UU(I) = U(I)
        END DO
 1100 CONTINUE
   80 rewind 2
      do j = 1,m
        write(2,96) u(j)
      end do
      CALL GETTIM (IHR, IMIN, ISEC, I100)  
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEB = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
      WRITE(*,90) IJ,ZDU, TIMEB-TIMEC
      DO I = 1,YSEL
        WRITE(6,*) Y(I)
      END DO
   90 FORMAT(' BEST SUBSET RELAXATION ',I4,F20.5,F10.2)
   96 FORMAT(F20.5)
  889 CONTINUE
      END
