      PROGRAM SANNEAL
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION AA(250,250),A(250,250),ASUM(250),SLEFT(250),
     1       X(250),COORD(250),SOLS(100,3),SOLS2(100,3),RMED(100,3),
     1       RES(250,250),LAMB1(250,250),LAMB2(250,250),CB(250,3)
      REAL S1,RCRIT
      INTEGER U(250),S(250),UNSEL,T(100,250),SB(250),Q,GB(250,3)
      CHARACTER JNK$
C
C ####################################################################
C  DHA.FOR IS A DISCRETE SIMULATED ANNEALING CODE FOR LEAST-SQAURES
C      UNIDIMENSIONAL SCALING BASED ON THE DEFAYS STRATEGY. THE CODE
C      USE PAIRWISE INTERCHANGE AND SINGLE-OBJECT INSERTIONS TO
C      GENERATE NEIGHBORING SOLUTIONS.
C   - 20 REPLICATIONS FROM RANDOM STARTS
C   - SUMMARY STATISTICS COMPUTED BY PROGRAM
C INPUT FILE: "MATFULL" AN N X N SYMMETRIC DISSIMILARITY MATRIX WHERE
C             THE FIRST ELEMENT IN THE FILE IS N.
C OUTPUT FILES: "OUTPUT" IS A SOLUTION SUMMARY
C               "STARTS" CONTAINS THE PERMUTATIONS FOR EACH REPLICATE,
C               WHICH CAN BE READ BY A POSTPROCESSOR SUCH AS THE "HAM"
C               STRATEGY.
C ####################################################################

      OPEN(1, FILE = 'DATASET')
      OPEN(2, FILE = 'OUTPUT')
      OPEN(3, FILE = 'PERMUT')
      OPEN(4, FILE = 'SUMMARY')
      CALL GETTIM (IHR, IMIN, ISEC, I100)
      CALL GETDAT (IYR, IMON, IDAY)  
      TIMEA = DFLOAT(86400*IDAY+3600*IHR+60*IMIN+ISEC)+DFLOAT(I100)/100.  
      RULE = .5D0    ! PROB OF SELECTING PAIRWISE INTERCHANGE
      READ(1,*) N
      READ(1,*) JNK$
      READ(1,*) ((AA(I,J),J=1,N),I=1,N)
      SSTOT = 0.
      SSY = 0.
      YBAR = 0.
      DO I = 1,N
        DO J = 1,N
          A(I,J) = AA(I,J)
          IF(I.LT.J) THEN
            SSTOT = SSTOT + AA(I,J)**2
            YBAR = YBAR + AA(I,J)
          END IF
        END DO
      END DO
      YBAR = YBAR/DFLOAT((N*(N-1))/2)
      DO I = 1,N-1
        DO J = I+1,N
          SSY = SSY+(AA(I,J)-YBAR)**2
        END DO
      END DO
      WRITE(*,177)
      WRITE(2,177)
C ###############################################################
C   BEGIN OUTER LOOP FOR THREE-SCALE EXTRACTION
C ###############################################################
C
      DO 3500 IQ = 1,3

      DO I = 1,N
        ASUM(I)=0.
      END DO
      EPS = 1.0D-08
      ZNORM = 0.0D0
      DO I = 1,N
        DO J = 1,N
          ASUM(I) = ASUM(I) + A(I,J)
          IF(I.LT.J) ZNORM = ZNORM + A(I,J)**2
        END DO
      END DO
      NREPS = 20
      DO 999 III = 1,NREPS
        DO I = 1,N
          U(I) = I
          T(III,I) = 0
        END DO
        UNSEL = N
        DO 1 I = 1,N
          CALL RANDOM(S1)
          ISET = S1 * FLOAT(UNSEL) + 1.
          IF(ISET.GT.UNSEL) ISET = UNSEL
          T(III,I) = U(ISET)
          DO J = ISET,UNSEL
            U(J) = U(J+1)
          END DO
          UNSEL = UNSEL - 1
  1     CONTINUE
 999  CONTINUE
C
      ZMIN = 9.9D+12
      ZAVG = 0.
      ZMAX = 0.
      DO 1000 III = 1,NREPS
        DO I = 1,N
          S(I) = T(III,I)
          SLEFT(I) = 0.0D0
        END DO
        Z = 0.0D0
        DO I = 1,N
          K = S(I)
          DO J = 1,I-1
            L = S(J)
            SLEFT(K) = SLEFT(K) + A(K,L)
          END DO
          Z = Z + (ASUM(K)-2*SLEFT(K))**2
        END DO
        ZBEST = Z
        TMAX = 0.0D0
        DO LLL = 1,5000
          CALL RANDOM(S1)
          I1 = S1 * FLOAT(N) + 1.
          IF(I1.GT.N) I1 = N
 199      CALL RANDOM(S1)
          J1 = S1 * FLOAT(N) + 1.
          IF(J1.GT.N) J1 = N
          IF(I1.EQ.J1) GO TO 199
          IF(I1.GT.J1) THEN
            JDUM = J1
            J1 = I1
            I1 = JDUM
          END IF
          K = S(I1)
          M = S(J1)
          DELTA = 0.0D0
          XK = 0.0D0
          XM = 0.0D0
          DO L = I1+1,J1
            Q = S(L)
            XK = XK + A(K,Q)
          END DO
          DO L = I1, J1-1
            Q = S(L)
            XM = XM - A(M,Q)
          END DO
          XK2 = XK * XK
          XM2 = XM * XM
          DELTA = DELTA + 4*XM2 + 8*XM*SLEFT(M) - 4*XM*ASUM(M)
          DELTA = DELTA + 4*XK2 + 8*XK*SLEFT(K) - 4*XK*ASUM(K)
          DO L = I1+1,J1-1
            Q = S(L)
            X(L) = A(M,Q) - A(K,Q)
            XL2 = X(L) * X(L)
            DELTA = DELTA + 4*XL2 + 8*X(L)*SLEFT(Q) - 4*X(L)*ASUM(Q)
          END DO
          IF(DELTA.LT.0) THEN
            IF(ABS(DELTA).GT.TMAX) TMAX = ABS(DELTA)
          END IF
        END DO
        COOL = .95
        TMIN = .0001d0
C        TMAX = Z
        ILOOP = 100*N
        NLOOP = (LOG(TMIN)-LOG(TMAX))/LOG(COOL)
C        WRITE(*,21) TMIN,TMAX,NLOOP
C  21    FORMAT(2F14.5,I6)
C        GO TO 889
        TEMP = TMAX
        DO I = 1,N
          SB(I) = S(I)
        END DO
C
        DO 2000 IJK = 1,NLOOP
          DO 2001 KKK = 1,ILOOP
            CALL RANDOM(S1)
            IF(S1.LE.RULE) THEN         ! INTERCHANGE
            CALL RANDOM(S1)
            I1 = S1 * FLOAT(N) + 1.
            IF(I1.GT.N) I1 = N
 99         CALL RANDOM(S1)
            J1 = S1 * FLOAT(N) + 1.
            IF(J1.GT.N) J1 = N
            IF(I1.EQ.J1) GO TO 99
            IF(I1.GT.J1) THEN
              JDUM = J1
              J1 = I1
              I1 = JDUM
            END IF
            K = S(I1)
            M = S(J1)
            DELTA = 0.0D0
            XK = 0.0D0
            XM = 0.0D0
            DO L = I1+1,J1
              Q = S(L)
              XK = XK + A(K,Q)
            END DO
            DO L = I1, J1-1
              Q = S(L)
              XM = XM - A(M,Q)
            END DO
            XK2 = XK * XK
            XM2 = XM * XM
            DELTA = DELTA + 4*XM2 + 8*XM*SLEFT(M) - 4*XM*ASUM(M)
            DELTA = DELTA + 4*XK2 + 8*XK*SLEFT(K) - 4*XK*ASUM(K)
            DO L = I1+1,J1-1
              Q = S(L)
              X(L) = A(M,Q) - A(K,Q)
              XL2 = X(L) * X(L)
              DELTA = DELTA + 4*XL2 + 8*X(L)*SLEFT(Q) - 4*X(L)*ASUM(Q)
            END DO
            IF(DELTA.GT.-EPS) THEN
              Z = Z + DELTA
              SLEFT(K) = SLEFT(K)+XK
              SLEFT(M) = SLEFT(M)+XM
              DO L = I1+1, J1-1
                Q = S(L)
                SLEFT(Q) = SLEFT(Q) + X(L)
              END DO
              S(I1) = M
              S(J1) = K
              IF(Z.GT.ZBEST) THEN
                ZBEST = Z
                DO I = 1,N
                  SB(I) = S(I)
                END DO
              END IF
            ELSE
              CALL RANDOM(S1)
              RCRIT = EXP(DELTA/TEMP)
              IF(S1.LE.RCRIT) THEN
                Z = Z + DELTA
                SLEFT(K) = SLEFT(K)+XK
                SLEFT(M) = SLEFT(M)+XM
                DO L = I1+1, J1-1
                  Q = S(L)
                  SLEFT(Q) = SLEFT(Q) + X(L)
                END DO
                S(I1) = M
                S(J1) = K
              END IF
            END IF

            ELSE                ! INSERTION

            CALL RANDOM(S1)
            I1 = S1 * FLOAT(N) + 1.      ! OBJECT POSITION IS I1
            IF(I1.GT.N) I1 = N
 599         CALL RANDOM(S1)
            J1 = S1 * FLOAT(N) + 1.
            IF(J1.GT.N) J1 = N
            IF(I1.EQ.J1) GO TO 599
            K = S(I1)
            DELTA = 0.0D0
            XK = 0.0D0
            IF(J1.GT.I1) THEN
              DO L = I1+1,J1
                Q = S(L)
                XK = XK + A(K,Q)
              END DO
            ELSE
              DO L = I1-1,J1,-1
                Q = S(L)
                XK = XK - A(K,Q)
              END DO
            END IF  
            XK2 = XK * XK
            DELTA = DELTA + 4*XK2 + 8*XK*SLEFT(K) - 4*XK*ASUM(K)
            IF(J1.GT.I1) THEN
              DO L = I1+1,J1
                Q = S(L)
                X(L) = -A(K,Q)
                XL2 = X(L) * X(L)
                DELTA=DELTA + 4*XL2 + 8*X(L)*SLEFT(Q) - 4*X(L)*ASUM(Q)
              END DO
            ELSE
              DO L = I1-1,J1,-1
                Q = S(L)
                X(L) = A(K,Q)
                XL2 = X(L) * X(L)
                DELTA=DELTA + 4*XL2 + 8*X(L)*SLEFT(Q) - 4*X(L)*ASUM(Q)
              END DO
            END IF
            IF(DELTA.GT.-EPS) THEN
              Z = Z + DELTA
              SLEFT(K) = SLEFT(K)+XK
              IF(J1.GT.I1) THEN
                DO L = I1+1, J1
                  Q = S(L)
                  SLEFT(Q) = SLEFT(Q) + X(L)
                END DO
                DO L = I1,J1-1
                  S(L)=S(L+1)
                END DO
                S(J1) = K
              ELSE
                DO L = I1-1,J1,-1
                  Q = S(L)
                  SLEFT(Q) = SLEFT(Q) + X(L)
                END DO
                DO L = I1,J1+1,-1
                  S(L)=S(L-1)
                END DO
                S(J1) = K
              END IF

              IF(Z.GT.ZBEST) THEN
                ZBEST = Z
                DO I = 1,N
                  SB(I) = S(I)
                END DO
              END IF
            ELSE
              CALL RANDOM(S1)
              RCRIT = EXP(DELTA/TEMP)
              IF(S1.LE.RCRIT) THEN
                Z = Z + DELTA
                SLEFT(K) = SLEFT(K)+XK
                IF(J1.GT.I1) THEN
                  DO L = I1+1, J1
                    Q = S(L)
                    SLEFT(Q) = SLEFT(Q) + X(L)
                  END DO
                  DO L = I1,J1-1
                    S(L)=S(L+1)
                  END DO
                  S(J1) = K
                ELSE
                  DO L = I1-1,J1,-1
                    Q = S(L)
                    SLEFT(Q) = SLEFT(Q) + X(L)
                  END DO
                  DO L = I1,J1+1,-1
                    S(L)=S(L-1)
                  END DO
                  S(J1) = K
                END IF
              END IF
            END IF

            END IF
 2001     CONTINUE
          TEMP = TEMP*COOL
 2000   CONTINUE
C
        ZLOSS = 0.0D0
        DO I = 1,N
          COORD(I) = 0.0D0
        END DO
        K = SB(1)
        COORD(K) = 0.0D0
        DO I = 2,N
          K = SB(I)
          K2 = SB(I-1)
          CRIT = 0.0D0
          DO L = I,N
            K3 = SB(L)
            CRIT = CRIT + A(K2,K3)
          END DO
          DO L = 1,I-2
            K3 = SB(L)
            CRIT = CRIT - A(K2,K3)
          END DO
          DO L = I+1,N
            K3 = SB(L)
            CRIT = CRIT - A(K,K3)
          END DO
          DO L = 1,I-1
            K3 = SB(L)
            CRIT = CRIT + A(K,K3)
          END DO
          CRIT = CRIT / DFLOAT(N)
          COORD(K) = COORD(K2) + CRIT
        END DO
        DO I = 1,N-1
          DO J = I+1,N
            DIST = ABS(COORD(I) - COORD(J))
            ZLOSS = ZLOSS + (A(I,J) - DIST)**2
          END DO
        END DO
        SOLS(III,IQ) = ZBEST
        SOLS2(III,IQ) = ZLOSS
        RMED(III,IQ) = ZLOSS
        ZSUM = ZSUM + ZLOSS
        IF(ZLOSS.LT.ZMIN) THEN
          ZMIN = ZLOSS
          DO I = 1,N
            GB(I,IQ) = SB(I)
            CB(I,IQ) = COORD(I)
          END DO
          DO I = 1,N-1
            DO J = I+1,N
              DIST = ABS(COORD(I) - COORD(J))
              IF(IQ.EQ.1) THEN
                LAMB1(I,J) = 1
                IF(A(I,J).LT.DIST) LAMB1(I,J) = -1
                LAMB1(J,I) = LAMB1(I,J)
              ELSEIF(IQ.EQ.2) THEN
                LAMB2(I,J) = 1
                DIST2 = LAMB1(I,J)*DIST+ABS(CB(I,1)-CB(J,1))
                IF(AA(I,J).LT.DIST2) LAMB2(I,J) = -1
                LAMB2(J,I) = LAMB2(I,J)
              END IF
              RES(I,J) = DIST-A(I,J)
              IF(RES(I,J).LT.0) RES(I,J) = -RES(I,J)
              RES(J,I) = RES(I,J)
            END DO
          END DO
        END IF
        IF(ZLOSS.GT.ZMAX) ZMAX = ZLOSS
        WRITE(*,77) iii,zbest,zloss,zloss/SSTOT
        WRITE(2,77) iii,zbest,zloss,zloss/SSTOT
 1000 CONTINUE
      ZSUM = ZSUM/DFLOAT(NREPS)
      IBEST = 0
      DO I = 1,NREPS
        IF(SOLS2(I,IQ).LT.ZMIN+EPS) IBEST = IBEST + 1
      END DO
      ATTR = 100.*(DFLOAT(IBEST) / DFLOAT(NREPS))
      DO I = 1,NREPS-1
        DO J = I + 1,NREPS
          IF(RMED(J,IQ).LT.RMED(I,IQ)) THEN
            RST = RMED(J,IQ)
            RMED(J,IQ) = RMED(I,IQ)
            RMED(I,IQ) = RST
          END IF
        END DO
      END DO
      NR = NREPS/2
      RMN = (RMED(NR,IQ)+RMED(NR+1,IQ))/2.0
      WRITE(2,*)
      WRITE(*,*)
      WRITE(2,201)
      WRITE(*,201)
      WRITE(2,202) ZMIN, ZSUM, RMN, ZMAX,ATTR
      WRITE(*,202) ZMIN, ZSUM, RMN, ZMAX,ATTR
      WRITE(2,203) ZMIN/SSTOT, ZSUM/SSTOT, RMN/SSTOT, ZMAX/SSTOT
      WRITE(4,204) ZMIN/SSTOT, ZSUM/SSTOT, RMN/SSTOT, ZMAX/SSTOT, ATTR
      WRITE(*,203) ZMIN/SSTOT, ZSUM/SSTOT, RMN/SSTOT, ZMAX/SSTOT
      DO I = 1,N
        DO J = 1,N
          A(I,J) = RES(I,J)
        END DO
      END DO
 3500 CONTINUE
C ############################################################
C  COMPLETED GENERATION OF THREE SCALES
C ################################################
C
      WRITE(3,*) N
      DO I = 1,N
        WRITE(3,207) I,GB(I,1),GB(I,2),GB(I,3),CB(GB(I,1),1),
     1               CB(GB(I,2),2),CB(GB(I,3),3)
      END DO
C
C DOUBLE CHECK RESULTS AFTER SECOND SCALE AND FIND CONSTANT
C
      RLOSS = 0.
      RMAX2 = -999999.
      RMIN2 = 999999.
      DO I = 1,N-1
        DO J = I+1,N
          DIST = ABS(CB(I,1)-CB(J,1))+LAMB1(I,J)*ABS(CB(I,2)-CB(J,2))
          IF(AA(I,J)-DIST.LT.RMIN2) RMIN2 = AA(I,J) - DIST
          IF(AA(I,J)-DIST.GT.RMAX2) RMAX2 = AA(I,J) - DIST
          RLOSS = RLOSS + (AA(I,J)-DIST)**2
        END DO
      END DO
      WRITE(*,*) RLOSS
      WRITE(*,209) RLOSS/SSTOT
      WRITE(2,209) RLOSS/SSTOT
      RSTEP = (RMAX2-RMIN2)/1000
      C2 = RMIN2-RSTEP
      RBEST = 999999999.
      DO K = 1,1000
        C2 = C2 + RSTEP
        RLOSS = 0.
        DO I = 1,N-1
          DO J = I+1,N
            DIST = ABS(CB(I,1)-CB(J,1))+LAMB1(I,J)*ABS(CB(I,2)-CB(J,2))
            DIST = DIST + C2
            RLOSS = RLOSS + (AA(I,J)-DIST)**2
          END DO
        END DO
        IF(RLOSS.LT.RBEST) THEN
          RBEST = RLOSS
          C2BEST = C2
        END IF
      END DO
C      WRITE(*,*) RBEST,SSY
      WRITE(*,210) C2BEST,1-RBEST/SSY,RBEST/SSTOT
      WRITE(2,210) C2BEST,1-RBEST/SSY,RBEST/SSTOT
C
C DOUBLE CHECK RESULTS AFTER THIRD SCALE AND FIND CONSTANT
C
      RLOSS = 0.
      RMAX3 = -9999999.
      RMIN3 = 9999999.
      DO I = 1,N-1
        DO J = I+1,N
          DIST = ABS(CB(I,1)-CB(J,1))+LAMB1(I,J)*ABS(CB(I,2)-CB(J,2))
          DIST = DIST+LAMB2(I,J)*ABS(CB(I,3)-CB(J,3))
          IF(AA(I,J)-DIST.LT.RMIN3) RMIN3 = AA(I,J) - DIST
          IF(AA(I,J)-DIST.GT.RMAX3) RMAX3 = AA(I,J) - DIST
          RLOSS = RLOSS + (AA(I,J)-DIST)**2
        END DO
      END DO
      write(*,*) rloss
      WRITE(*,211) RLOSS/SSTOT
      WRITE(2,211) RLOSS/SSTOT
      RSTEP = (RMAX3-RMIN3)/1000
      C3 = RMIN3-RSTEP
      RBEST = 999999999.
      DO K = 1,1000
        C3 = C3 + RSTEP
        RLOSS = 0.
        DO I = 1,N-1
          DO J = I+1,N
            DIST = ABS(CB(I,1)-CB(J,1))+LAMB1(I,J)*ABS(CB(I,2)-CB(J,2))
            DIST = DIST+LAMB2(I,J)*ABS(CB(I,3)-CB(J,3))
            DIST = DIST + C3
            RLOSS = RLOSS + (AA(I,J)-DIST)**2
          END DO
        END DO
        IF(RLOSS.LT.RBEST) THEN
          RBEST = RLOSS
          C3BEST = C3
        END IF
      END DO
      WRITE(*,212) C3BEST,1-RBEST/SSY,RBEST/SSTOT
      WRITE(2,212) C3BEST,1-RBEST/SSY,RBEST/SSTOT
      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(*,208) TIMTOT
  177 FORMAT(2X,'REP',4X,'DEFAYS CRIT',2X,'L2 LOSS VALUE')
   77 format(I5,2f15.4,f12.8)
  200 FORMAT(I3,2F20.6)
  201 FORMAT(9X,'MIN',11X,'AVG',11X,'MED',10X,'MAX',5X,'ATT',3X,'TIME')
  202 FORMAT(4F14.4,F5.0)
  203 FORMAT(4F14.8)
  204 FORMAT(4F14.8,F5.0,F7.2)
  207 FORMAT(4I4,3F10.5)
  208 FORMAT(' TOTAL TIME = ',F8.2)
  209 FORMAT(' DOUBLE CHECK AFTER SECOND SCALE ',F12.8)
  210 FORMAT(' CONSTANT 2 ',F12.6,' VAF = ',F12.8,' STRESS = ',F12.8)
  211 FORMAT(' DOUBLE CHECK AFTER THIRD SCALE ',F12.8)
  212 FORMAT(' CONSTANT 3 ',F12.6,' VAF = ',F12.8,' STRESS = ',F12.8)
  889 continue
      END
