PROGRAM INWARD
      IMPLICIT INTEGER(A-Z)
      DOUBLE PRECISION TIMEA,TIMEB,TIMTOT,COORD(50),SP(50),Z2,Z3,Z1,Z9,
     1                 RSET,A(50,50),ASUM(50),T(50),Z,ZBD,ZS,ZT,Z6,
     1        RNK(50,50),BND(50,50),MAXBND(50),MAXBND2(50),MAXB(50),
     1        SLEFT(50),SRIGHT(50),IDX1,IDX2,K1A,K2A,K1B,SL,XR3,EPS,
     1        K2B,ISUM,ZVAL,DELTA,ISUM1,ISUM2,IBEST,IBEST2,JKL,ZBEST
      REAL S1
      INTEGER Q(50),D(50),UNSEL(50),X(50),XX(50)
C
C #################################################################
C 9/24/02 This program fits least-squares unidimensional scaling in
C         the L1 norm using an "ends-in" branching approach similar toC         cuts.  An improved symmetry check was also incorporated.
C         Defays(1978) approach.
C    - The new features of this program are:
C      1) Improved bounding procedures
C      2) The "Interchange Test"
C #################################################################
C
      EPS = 1.0D-09
      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.1.0D-12) 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
C        write(*,*) 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
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
          DO L = K+1,ICT
            IF(RNK(I,L).LT.RNK(I,K)) THEN
              IDUM=RNK(I,K)
              RNK(I,K)=RNK(I,L)
              RNK(I,L)=IDUM
            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
          ISUM1=0
          DO K = 1,J-1
            ISUM1=ISUM1+RNK(I,K)
          END DO
          ISUM2=ASUM(I)-ISUM1
          BND(I,J)=(ISUM2-ISUM1)**2
C
          ISUM1=0
          DO K = 1,J-1
            KK=N-K
            ISUM1=ISUM1+RNK(I,KK)
          END DO
          ISUM2=ASUM(I)-ISUM1
          IF((ISUM2-ISUM1)**2.GT.BND(I,J)) BND(I,J)=(ISUM2-ISUM1)**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
C        WRITE(2,*) J,MAXBND(J)
      END DO
      ISET=FLOAT(N)/2.+.6
      ICT=-1
      DO I = 1,ISET
        ICT=ICT+2
        MAXBND2(ICT)=MAXBND(I)
      END DO
      ICT=0
      DO I = N,ISET+1,-1
        ICT=ICT+2
        MAXBND2(ICT)=MAXBND(I)
      END DO
C
      DO J = 1,N-1
        DO I = J,N
          MAXB(J)=MAXB(J)+MAXBND2(I)
        END DO
      END DO
C
C
C  ############### BRANCH-AND-BOUND ALGORITHM BEGINS HERE #############
C
      M=1
      Q(M)=1
      D(1)=1
      DO I = 2,N
        SLEFT(I)=A(1,I)
        SRIGHT(I)=0
      END DO
      T(M)=0
      SIGN=1
      ZVAL=ASUM(1)**2
      DO K = 2,N
        Q(K)=0
      END DO
C
  1   M = M + 1                             ! BRANCH STEP
      SIGN=1-SIGN
      GO TO 2
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
      IF(M.EQ.2.AND.Q(2).LT.Q(1)) GO TO 2   ! SYMMETRY FATHOM CHECK
 22   R3=Q(M)
      T(M)=SLEFT(R3)*SIGN+SRIGHT(R3)*(1-SIGN)
      D(R3)=1
      IF(M.LE.2) THEN                      ! IF M <=2 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)*SIGN
          SRIGHT(I)=SRIGHT(I)+A(I,R3)*(1-SIGN)
  220   CONTINUE
        GO TO 1
      END IF
      IF(M.EQ.N-1) THEN                    ! COMPLETE SEQUENCE - EVALUATE
        ZBD=0
        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
        R2=Q(N)
        DO I = N-2,1,-2
          R1=Q(I)
          T(N)=T(N)+A(R1,R2)
        END DO
        DO I = 1,N
          R1=Q(I)
          ISUM=ASUM(R1)-T(I)
          ZBD=ZBD+(ISUM-T(I))**2
        END DO
        IF(ZBD.GT.Z) THEN
          Z=ZBD
c          write(*,*) 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                                  ! RULE #1 ADJACENCY TEST
        R2=Q(M-2)
        R4=Q(M-1)
        IDX1=ASUM(R2)-2*T(M-2)-2*A(R2,R3)
        IDX2=ASUM(R3)-2*T(M)
        IDX3=2*T(M-1)-ASUM(R4)
        IF(IDX1.GE.IDX2.AND.IDX2.GE.IDX3) GO TO 91
        D(R3)=0
        GO TO 2
 91     Z1=IDX2**2
        ZS=ZVAL+Z1

C
        R4=Q(M-1)                           ! RULE #2: DEFAYS SIMPLE BOUND
        ISUM2=ASUM(R4)-T(M-1)
        Z9=(ISUM2-T(M-1))**2
        Z6=Z1
        IF(Z9.GT.Z6) Z6=Z9
        Z6=(N-M)*Z6
        IF(ZS+Z6.LT.Z) THEN
          D(R3)=0
          GO TO 2
        END IF
C
        ZT=MAXB(M+1)                        ! RULE #3: MY SIMPLE BOUND
        IF(ZS+ZT.LT.Z) THEN
          D(R3)=0
          GO TO 2
        END IF
C
        IF(SIGN.EQ.1) THEN                  ! RULE #4: INTERCHANGE TEST
        DO L = M-4,1,-2            ! INTERCHANGING 'M' WITH OBJECTS ON LEFT
          R2=Q(L)
          DELTA=0
          DELTA=DELTA-(ASUM(R3)-2*T(M))**2    ! Loss from moving R3 to L
          DELTA=DELTA-(ASUM(R2)-2*T(L))**2    ! Loss from moving R2 to M
          JKL=T(M)
          DO LL = L,M-2,2
            R1=Q(LL)
            JKL=JKL-A(R3,R1)
          END DO
          DELTA=DELTA+(ASUM(R3)-2*JKL)**2     ! Gain from moving R3 to L
          JKL=T(L)
          DO LL = L+2,M,2
            R1=Q(LL)
            JKL=JKL+A(R2,R1)                  ! GAIN FROM MOVING R2 TO M
          END DO
          DELTA=DELTA+(ASUM(R2)-2*JKL)**2
          DO LL = L+2,M-2,2
            R1=Q(LL)
            DELTA=DELTA-(ASUM(R1)-2*T(LL))**2
            JKL=T(LL)+A(R3,R1)-A(R2,R1)
            DELTA=DELTA+(ASUM(R1)-2*JKL)**2
          END DO
          IF(DELTA.GT.0) THEN
            D(R3)=0
            GO TO 2
          END IF
        END DO
C
        DO L = M-1,2,-2     ! INTERCHANGING 'M' WITH OBJECTS ON RIGHT
          R2=Q(L)
          DELTA=0
          DELTA=DELTA-(ASUM(R3)-2*T(M))**2    ! Loss from moving R3 to L
          DELTA=DELTA-(ASUM(R2)-2*T(L))**2    ! Loss from moving R2 to M
          JKL=0
          DO LL = L-2,2,-2
            R1=Q(LL)
            JKL=JKL+A(R3,R1)
          END DO
          DELTA=DELTA+(ASUM(R3)-2*JKL)**2     ! Gain from moving R3 to L
          JKL=0
          DO LL = 1,M-2,2
            R1=Q(LL)
            JKL=JKL+A(R2,R1)                  ! GAIN FROM MOVING R2 TO M
          END DO
          DELTA=DELTA+(ASUM(R2)-2*JKL)**2
          DO LL = L,M-1,2
            R1=Q(LL)
            DELTA=DELTA-(ASUM(R1)-2*T(LL))**2
            JKL=T(LL)+A(R3,R1)-A(R2,R1)
            DELTA=DELTA+(ASUM(R1)-2*JKL)**2
          END DO
          IF(DELTA.GT.0) THEN
            D(R3)=0
            GO TO 2
          END IF
        END DO
C
        ELSE
C
        DO L = M-4,2,-2    ! INTERCHANGING 'M' WITH OBJETS ON RIGHT
          R2=Q(L)
          DELTA=0
          DELTA=DELTA-(ASUM(R3)-2*T(M))**2    ! Loss from moving R3 to L
          DELTA=DELTA-(ASUM(R2)-2*T(L))**2    ! Loss from moving R2 to M
          JKL=T(M)
          DO LL = L,M-2,2
            R1=Q(LL)
            JKL=JKL-A(R3,R1)
          END DO
          DELTA=DELTA+(ASUM(R3)-2*JKL)**2     ! Gain from moving R3 to L
          JKL=T(L)
          DO LL = L+2,M,2
            R1=Q(LL)
            JKL=JKL+A(R2,R1)
          END DO
          DELTA=DELTA+(ASUM(R2)-2*JKL)**2
          DO LL = L+2,M-2,2
            R1=Q(LL)
            DELTA=DELTA-(ASUM(R1)-2*T(LL))**2
            JKL=T(LL)+A(R3,R1)-A(R2,R1)
            DELTA=DELTA+(ASUM(R1)-2*JKL)**2
          END DO
          IF(DELTA.GT.0) THEN
            D(R3)=0
            GO TO 2
          END IF
        END DO
C
        DO L = M-1,1,-2     ! INTERCHANGING 'M' WITH OBJECTS ON LEFT
          R2=Q(L)
          DELTA=0
          DELTA=DELTA-(ASUM(R3)-2*T(M))**2    ! Loss from moving R3 to L
          DELTA=DELTA-(ASUM(R2)-2*T(L))**2    ! Loss from moving R2 to M
          JKL=0
          DO LL = L-2,1,-2
            R1=Q(LL)
            JKL=JKL+A(R3,R1)
          END DO
          DELTA=DELTA+(ASUM(R3)-2*JKL)**2     ! Gain from moving R3 to L
          JKL=0
          DO LL = 2,M-2,2
            R1=Q(LL)
            JKL=JKL+A(R2,R1)                  ! GAIN FROM MOVING R2 TO M
          END DO
          DELTA=DELTA+(ASUM(R2)-2*JKL)**2
          DO LL = L,M-1,2
            R1=Q(LL)
            DELTA=DELTA-(ASUM(R1)-2*T(LL))**2
            JKL=T(LL)+A(R3,R1)-A(R2,R1)
            DELTA=DELTA+(ASUM(R1)-2*JKL)**2
          END DO
          IF(DELTA.GT.0) THEN
            D(R3)=0
            GO TO 2
          END IF
        END DO
        END IF
C
        Z3=0                                  ! RULE #4 MY BOUND TEST
        DO 90 I = 1,N
          IF(D(I).EQ.1) GO TO 90
          ISUM1=SLEFT(I)
          ISUM2=SRIGHT(I)
          ISUM1=ISUM1+A(I,R3)*SIGN
          ISUM2=ISUM2+A(I,R3)*(1-SIGN)
          IBEST=ISUM1
          IF(ISUM2.LT.IBEST) IBEST=ISUM2
          IBEST2=ASUM(I)-IBEST
          Z3=Z3+(IBEST2-IBEST)**2
  90    CONTINUE
        IF(ZS+Z3.LT.Z) THEN
          D(R3)=0
          GO TO 2
        END IF
C
  92    ZVAL=ZS
        DO 225 I = 1,N
          IF(D(I).EQ.1) GO TO 225
          SLEFT(I)=SLEFT(I)+A(I,R3)*SIGN
          SRIGHT(I)=SRIGHT(I)+A(I,R3)*(1-SIGN)
 225    CONTINUE
        SLEFT(R3)=0
        SRIGHT(R3)=0
        GO TO 1
      END IF
C
   7  D(Q(M))=0                             ! ACTUAL RETRACTION
      Q(M)=0
      M=M-1
      SIGN=1-SIGN
      R2=Q(M)
      ISUM2=ASUM(R2)-T(M)
      ZVAL=ZVAL-(ISUM2-T(M))**2
      IF(SIGN.EQ.0) THEN
        DO 226 I = 1,N
          IF(D(I).EQ.1) GO TO 226
          SRIGHT(I)=SRIGHT(I)-A(I,R2)
 226    CONTINUE
        SRIGHT(R2)=T(M)
        SLEFT(R2)=0
        DO I = 1,M,2
          R1=Q(I)
          SLEFT(R2)=SLEFT(R2)+A(R1,R2)
        END DO
      ELSE
        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)
        SRIGHT(R2)=0
        DO I = 2,M,2
          R1=Q(I)
          SRIGHT(R2)=SRIGHT(R2)+A(R1,R2)
        END DO
      END IF
      D(R2)=0
      GO TO 2
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
      DOG=-1
      CAT=0
      ICT=0
      ICT2=N+1
 86   ICT=ICT+1
      ICT2=ICT2-1
      DOG=DOG+2
      CAT=CAT+2
      XX(ICT)=X(DOG)
      XX(ICT2)=X(CAT)
      RSET=FLOAT(N)/2.
      IF(FLOAT(ICT+1).LE.RSET+.01) GO TO 86
      ISET=N/2
      IF(RSET-FLOAT(ISET).GT..01) XX(ISET+1)=X(N)
C
      DO I = 1,N-1
        K1A=0
        K2A=0
        R2=XX(I)
        R3=XX(I+1)
        DO J = 1,I-1
          R1=XX(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,XX(I),COORD(I)
      END DO
      Z2=0.0D0
      DO I = 1,N-1
        R1=XX(I)
        DO J = I+1,N
          R2=XX(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 ',F18.5,' CPU TIME ',F8.2)
 70   FORMAT(2I3,F12.5)
 71   FORMAT(' MINIMUM LOSS FUNCTION VALUE ', F18.5)
C
      END