PROGRAM FORWARD
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION TIMEA,TIMEB,TIMTOT,COORD(50),SP(50),Z2,ZBEST,EPS,
     1       A(50,50),ASUM(50),RNK(50,50),MAXB(50),MAXBND(50),BND(50,50)
     1       ,SLEFT(50),T(50),RSUM,RSUM1,RSUM2,RDUM,ZVAL,Z,ZBD,RDX1,XLL,
     1       RDX2,DELTA,K1A,K1B,K2A,K2B,Z1,ZS,ZT,XR3,XR2,XL,XS2,XS3,SL
      REAL S1
      INTEGER Q(50),D(50),X(50),UNSEL(50)
C
C #################################################################
C 9/24/02 FORWARD LEAST SQUARES UNIDIMENSIONAL SCALING ALGORITHM
C         - This program is similar to "Inward", except that branching
C           is strictly left-to-right.
C #################################################################
C
      EPS = 1.0D-9
      OPEN(1,FILE='AMAT.DAT')          ! Dissimilarity matrix
      OPEN(2,FILE='SEQ.OUT')           ! Output file
      READ(1,*) N                      ! Read number of objects
      WRITE(*,*) 'TYPE 1 FOR HALF MATRIX OR TYPE 2 FOR FULL MATRIX'
      READ(*,*) ITYPE
      IF(ITYPE.EQ.2) THEN
        READ(1,*) ((A(I,J),J=1,N),I=1,N)
      ELSE
        DO J = 2,N
          READ(1,*) (A(I,J),I=1,J-1)
        END DO
        DO J = 2,N
          DO I = 1,J-1
            A(J,I) = A(I,J)
          END DO
        END DO
      END IF
      IFLAG = 0
      DO J = 2,N
        DO I = 1,J-1
          IF(A(I,J).LT.EPS) IFLAG = 1
        END DO
      END DO
      IF(IFLAG.EQ.1) THEN
        WRITE(*,*) 'THERE IS A ZERO VALUE IN THE DISSIMILARITY MATRIX'
        WRITE(*,*) 'DO YOU WANT TO IMPOSE A MILD REGULARITY CONDITION'
        WRITE(*,*) ' 1 = YES AND 0 = NO'
        READ(*,*) IFLAG2
        IF(IFLAG2.EQ.1) THEN
          DO J = 2,N
            DO I = 1,J-1
              IF(A(I,J).LT.EPS) A(I,J) = EPS
              A(J,I) = A(I,J)
            END DO
          END DO
        END IF
      END IF
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEA=DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
      DO I = 1,N
        A(I,I) = 0.0D0
      END DO
      DO I = 1,N                       ! The value stored in ASUM(I) is the
        ASUM(I)=0.0D0                  ! sum of the dissimilarity measures
        DO J = 1,N                     ! for row I.
          ASUM(I)=ASUM(I)+A(I,J)
        END DO
      END DO
C
      ZBEST = 0.0D0
      DO 3500 JJJ = 1,20
        DO I = 1,N
          UNSEL(I) = I
          Q(I) = 0
        END DO
        NNSEL = N
 3501   CALL RANDOM(S1)
        ISEL = 1. + S1*FLOAT(NNSEL)
        IF(ISEL.GT.NNSEL) ISEL = NNSEL
        Q(NNSEL) = UNSEL(ISEL)
        DO J = ISEL,NNSEL-1
          UNSEL(J) = UNSEL(J+1)
        END DO
        NNSEL = NNSEL - 1
        IF(NNSEL.GT.0) GO TO 3501
        Z = 0.0D0
        DO I = 1,N
          SL = 0.0D0
          DO J = 1,I-1
            SL = SL + A(Q(I),Q(J))
          END DO
          Z = Z + (2*SL-ASUM(Q(I)))**2
        END DO
 3502   ITRIG = 0
        DO II = 1,N-1
          DO JJ = II+1,N
            R3 = Q(JJ)
            R2 = Q(II)
            DELTA=0.
            XR3 = 0.0D0
            DO LL = 1,JJ-1
              IF(LL.NE.II) XR3 = XR3 + A(R2,Q(LL))
            END DO
            XR3 = XR3 + A(R2,R3)
            DELTA = (2*XR3-ASUM(R2))**2
            XR3 = 0.0D0
            DO LL = 1,II-1
              XR3 = XR3 + A(R3,Q(LL))
            END DO
            DELTA = DELTA + (2*XR3-ASUM(R3))**2
            DO LL = II+1,JJ-1
              XR3 = 0.0D0
              DO L = 1,LL-1
                XR3 = XR3 + A(Q(LL),Q(L))
              END DO
              XR3 = XR3 + A(Q(LL),R3)-A(Q(LL),R2)
              DELTA = DELTA + (2*XR3-ASUM(Q(LL)))**2
            END DO
C
            XR3 = 0.0D0
            DO LL = 1,JJ-1
              XR3 = XR3 + A(R3,Q(LL))
            END DO
            DELTA = DELTA - (2*XR3-ASUM(R3))**2
            XR3 = 0.0D0
            DO LL = 1,II-1
              XR3 = XR3 + A(R2,Q(LL))
            END DO
            DELTA = DELTA - (2*XR3-ASUM(R2))**2
            DO LL = II+1,JJ-1
              XR3 = 0.0D0
              DO L = 1,LL-1
                XR3 = XR3 + A(Q(LL),Q(L))
              END DO
              DELTA = DELTA - (2*XR3-ASUM(Q(LL)))**2
            END DO
            IF(DELTA.GT.EPS) THEN
              Z = Z + DELTA
              Q(II) = R3
              Q(JJ) = R2
              ITRIG = 1
            END IF
          END DO
        END DO
        IF(ITRIG.EQ.1) GO TO 3502
        IF(Z.GT.ZBEST) ZBEST = Z
 3500 CONTINUE
      WRITE(*,3505) ZBEST
 3505 FORMAT(' HEURISTIC DEFAYS CRIT VALUE ',F20.4)
      Z = ZBEST-.01
      DO I = 1,N
        Q(I) = 0
      END DO
 31   FORMAT(F17.5)
C
C ####### FIND THE RANK ORDER FOR EACH OBJECT
C
      DO 60 I = 1,N
        ICT=0
        DO 61 J = 1,N
          IF(I.EQ.J) GO TO 61
          ICT=ICT+1
          RNK(I,ICT) = A(I,J)
 61     CONTINUE
        DO K = 1,ICT-1                      ! Each row of the matrix RNK
          DO L = K+1,ICT                    ! contains the rank-ordered
            IF(RNK(I,L).LT.RNK(I,K)) THEN   ! (ascending) matrix elements
              RDUM=RNK(I,K)                 ! for that row.
              RNK(I,K)=RNK(I,L)
              RNK(I,L)=RDUM
            END IF
          END DO
        END DO
 60   CONTINUE
C      WRITE(2,97) ((RNK(I,J),J=1,20),I=1,20)
C 97   FORMAT(20I4)
C
C ####### FIND AN UPPER BOUND FOR EACH OBJECT IN EACH POSITION
C
      DO I = 1,N
        DO J = 1,N                    ! The matrix element BND(I,J) contains
          RSUM1=0                     ! the very best contribution to the
          DO K = 1,J-1                ! index value that object I could make
            RSUM1=RSUM1+RNK(I,K)      ! if it were placed in position J.
          END DO
          RSUM2=ASUM(I)-RSUM1
          BND(I,J)=(RSUM2-RSUM1)**2
C
          RSUM1=0
          DO K = 1,J-1
            KK=N-K
            RSUM1=RSUM1+RNK(I,KK)
          END DO
          RSUM2=ASUM(I)-RSUM1
          IF((RSUM2-RSUM1)**2.GT.BND(I,J)) BND(I,J)=(RSUM2-RSUM1)**2
        END DO
      END DO
C
      DO J = 1,N
        MAXBND(J)=0
        DO I = 1,N
          IF(BND(I,J).GT.MAXBND(J)) MAXBND(J)=BND(I,J)
        END DO
      END DO
      DO J = 1,N-1
        DO I = J,N
          MAXB(J)=MAXB(J)+MAXBND(I)
        END DO
      END DO
C
C  ############### BRANCH-AND-BOUND ALGORITHM BEGINS HERE #############
C
C  Some quick notes:
C    1) M is the sequence position pointer.  Initially set to 1.
C    2) Q(M) is the object stored in position M. So if Q(2) = 6, this
C       means that object 6 is stored in position 2 of the sequence.
C    3) D(I) takes on a value of a if object I has been assigned a
C       position in the partial sequence, otherwise D(I) = 0.
C    4) ZVAL is the objective index value of the partial sequence.
C    5) TRIG = 1 if object 1 has been assigned, 0 if it has not.
C    6) SLEFT(I) = 0 if object I has been assigned in the partial sequence,
C       otherwise it equals the sum of dissimilarities between object I
C       and all assigned objects.
C    7) T(M) = the sum of the dissimilarities to the left of postion M.
C
      M=1                  ! M is initially set to 1
      Q(M)=1               ! Object 1 is placed in position 1
      D(1)=1               ! Object 1 is assigned
      TRIG=1
      DO I = 2,N
        SLEFT(I)=A(1,I)
      END DO
      T(M)=0
      ZVAL=ASUM(1)**2
      DO K = 2,N
        Q(K)=0
      END DO
C
  1   M = M + 1                             ! BRANCH STEP
C
  2   Q(M)=Q(M)+1                           ! NEXT OBJ. AT CURRENT TREE LEVEL
C
      IF(D(Q(M)).EQ.1) GO TO 2
      IF(M.EQ.1.AND.Q(M).GT.N) GO TO 9      ! TERMINATE
      IF(M.GT.1.AND.Q(M).GT.N) GO TO 7      ! GO TO RETRACTION
      R3=Q(M)
      IF(R3.EQ.2.AND.TRIG.EQ.0) GO TO 2     ! SYMMETRY FATHOM CHECK
  22  D(R3)=1
      T(M)=SLEFT(R3)
      IF(M.EQ.1) THEN                       ! IF M = 1 THEN BRANCH
        ZVAL=ZVAL+ASUM(R3)**2
        DO 220 I = 1,N
          IF(D(I).EQ.1) GO TO 220
          SLEFT(I)=SLEFT(I)+A(I,R3)
  220   CONTINUE
        GO TO 1
      END IF
      IF(M.EQ.N-1) THEN                    ! COMPLETE SEQUENCE - EVALUATE
        ZBD=0.0D0
        DO 185 I = 1,N
          IF(D(I).EQ.1) GO TO 185
          Q(N)=I
          GO TO 188
 185    CONTINUE
 188    T(N)=0.0D0
        R2=Q(N)
        DO I = N-1,1,-1
          R1=Q(I)
          T(N)=T(N)+A(R1,R2)
        END DO
        DO I = 1,N
          R1=Q(I)
          RSUM=ASUM(R1)-T(I)
          ZBD=ZBD+(RSUM-T(I))**2
        END DO
        IF(ZBD.GT.Z) THEN
          Z=ZBD
          write(*,31) z
          DO I = 1,N
            X(I)=Q(I)
          END DO
        END IF
        Q(N)=0
        D(R3)=0
        GO TO 2
C
C ############################################
C  THE BOUNDING RULES ARE ENACTED
C ############################################
C
      ELSE                                    ! ADJACENCY TEST
        R2=Q(M-1)
        RDX1=ASUM(R2)-2*T(M-1)-2*A(R2,R3)
        RDX2=ASUM(R3)-2*T(M)
        IF(RDX1.GE.RDX2) GO TO 98
        D(R3)=0
        GO TO 2
C
C        GO TO 198
 98     DO L = M-2,1,-1                       ! INTERCHANGE TEST
          R2=Q(L)
          DELTA=0.
          XR2 = 0.
          XR3 = 0.
          DO LL = L,M-1
            R1=Q(LL)
            XR3 = XR3 - A(R3,R1)
          END DO
          DO LL = L+1,M
            R1=Q(LL)
            XR2 = XR2 + A(R2,R1)
          END DO
          DELTA=DELTA+4*XR3*XR3+8*XR3*T(M)-4*XR3*ASUM(R3)
          DELTA=DELTA+4*XR2*XR2+8*XR2*T(L)-4*XR2*ASUM(R2)
          DO LL = L+1,M-1
            R1=Q(LL)
            XL = A(R1,R3) - A(R1,R2)
            DELTA=DELTA+4*XL*XL+8*XL*T(LL)-4*XL*ASUM(R1)
          END DO
          IF(DELTA.GT.0) THEN
            D(R3)=0
            GO TO 2
          END IF
        END DO
C
        GO TO 97
 198    DO L = M-2,1,-1                       ! INSERTION TEST
          DELTA=0.
          XR3 = 0.
          DO LL = M-1,L,-1
            R1=Q(LL)
            XR3 = XR3 - A(R3,R1)
          END DO
          DELTA=DELTA+4*XR3*XR3+8*XR3*T(M)-4*XR3*ASUM(R3)
          DO LL = M-1,L,-1
            R1=Q(LL)
            XL = A(R1,R3)
            DELTA=DELTA+4*XL*XL+8*XL*T(LL)-4*XL*ASUM(R1)
          END DO
          IF(DELTA.GT.0) THEN
            D(R3)=0
            GO TO 2
          END IF
        END DO
C
 97     Z1=RDX2**2
        ZS=ZVAL+Z1
C
        ZT=MAXB(M+1)                 ! SIMPLE BOUND TEST
        IF(ZS+ZT.LT.Z) THEN          ! Fails if < Z
          D(R3)=0
          GO TO 2
        END IF
C
  92    ZVAL=ZS                      ! If we make it to here, then all the
        DO 225 I = 1,N               ! interchange tests and bound tests
          IF(D(I).EQ.1) GO TO 225    ! have passed.  Thus, the search will
          SLEFT(I)=SLEFT(I)+A(I,R3)  ! move deeper into the tree and we
 225    CONTINUE                     ! branch by going to Step 1.
        SLEFT(R3)=0
        IF(R3.EQ.1) TRIG=1
        GO TO 1
      END IF
C
   7  D(Q(M))=0                   ! ACTUAL DEPTH RETRACTION
      Q(M)=0                      ! Object in position M set to ZERO
      M=M-1                       ! Back up one position in sequence (retract)
      R2=Q(M)                     ! Update values of SLEFT as needed
      RSUM2=ASUM(R2)-T(M)
      ZVAL=ZVAL-(RSUM2-T(M))**2
      DO 227 I = 1,N
        IF(D(I).EQ.1) GO TO 227
        SLEFT(I)=SLEFT(I)-A(I,R2)
 227  CONTINUE
      SLEFT(R2)=T(M)
      D(R2)=0
      IF(R2.EQ.1) TRIG=0
      GO TO 2
C
C  #############################################################
C  Branch-and-Bound is done when we get to statement 9 below.
C  The remainder of the program just computes the optimal
C  coordinates for each object given that the optimal ordering has
C  been found.
C ##############################################################
C
   9  CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)
      TIMEB=DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.
      TIMTOT=TIMEB-TIMEA
      WRITE(*,69) Z,TIMTOT
      WRITE(2,69) Z,TIMTOT
C
      DO I = 1,N-1
        K1A=0
        K2A=0
        R2=X(I)
        R3=X(I+1)
        DO J = 1,I-1
          R1=X(J)
          K1A=K1A+A(R1,R2)
          K2A=K2A+A(R1,R3)
        END DO
        K2A=K2A+A(R2,R3)
        K1B=ASUM(R2)-K1A
        K2B=ASUM(R3)-K2A
        SP(I)=(K1B-K1A)-(K2B-K2A)
        SP(I)=SP(I)/DFLOAT(N)
        COORD(I+1)=COORD(I)+SP(I)
      END DO
      DO I = 1,N
        WRITE(2,70) I,X(I),COORD(I)
      END DO
      Z2=0.0D0
      DO I = 1,N-1
        R1=X(I)
        DO J = I+1,N
          R2=X(J)
          Z2=Z2+(ABS(COORD(I)-COORD(J))-A(R1,R2))**2
        END DO
      END DO
      WRITE(*,71) Z2 
      WRITE(2,71) Z2
 69   FORMAT(' MAXIMUM DEFAYS CRITERION VALUE ',F17.5,' CPU TIME ',F8.2)
 70   FORMAT(2I3,F12.5)
 71   FORMAT(' MINIMUM LOSS FUNCTION VALUE ', F18.5)
C
      END