C
C-----------------------------------------------------------------------
C  MAIN PROGRAM:  FWIIR
C  AUTHORS:       KENNETH STEIGLITZ AND BRUCE D. LADENDORF
C                 PRINCETON UNIVERSITY, PRINCETON, NJ 08540
C                 VERSION OCTOBER 15, 1978
C
C  INPUT:         A GRID OF FREQUENCY POINTS AND A SET OF
C                 (ESSENTIALLY) INFINITE PRECISION COEFFICIENTS.
C                 SEE SUBROUTINE INDATA FOR DETAILED INPUT INFORMATION
C
C  THIS PROGRAM DESIGNS FINITE WORD-LENGTH IIR DIGITAL FILTERS.
C  THE METHOD USES RANDOMIZED PATTERN SEARCH OF HOOKE AND JEEVES,
C  AND IS DESCRIBED IN "DESIGNING SHORT-WORD RECURSIVE DIGITAL FILTERS,"
C  BY KENNETH STEIGLITZ, IN PROC. 9TH ANNUAL ALLERTON CONF. ON CIRCUIT
C  AND SYSTEM TH., PP. 778-788; OCT. 1971. REPRINTED IN DIGITAL SIGNAL
C  PROCESSING II, EDITED BY THE DIGITAL SIGNAL PROCESSING COMMITTEE
C  OF THE GROUP ON ASSP, IEEE PRESS, N. Y., 1976.
C-----------------------------------------------------------------------
C
      DOUBLE PRECISION W(100), Y(100), WEIGHT(100), X(36)
      DOUBLE PRECISION SDELTA, DELTAZ, DUMMY, RHO, EST, F
      LOGICAL PRINT, TWOPT, FIXCOF
      COMMON /ABC/ X, NSET, NTWOPT, NBIT1, NBIT2, NSECT, N
      COMMON /RAW/ W, Y, WEIGHT, M, KOUNT
      COMMON /STEER/ PRINT, TWOPT, FIXCOF
      COMMON /MACH/ IND, IOUTD
C
C SET UP MACHINE CONSTANTS
C
      IND = I1MACH(1)
      IOUTD = I1MACH(2)
C
C READ IN PROBLEM DATA
C
      CALL INDATA
C
C SET SMALL AND LARGE DELTA
C
      SDELTA = (.5D0)**NBIT1
      DELTAZ = (.5D0)**NBIT2
C
C SET RANDOM NUMBER GENERATOR BY CALLING IT NSET TIMES
C
      DO 10 J=1,NSET
        DUMMY = UNI(0.)
  10  CONTINUE
C
C TURN OFF THE PRINTING OPTION FOR FUNCT
C
      PRINT = .FALSE.
C
C ROUND OFF THE COEFFICIENTS TO NBIT1 BITS
C
      DO 20 J=1,N
        X(J) = DSIGN(SDELTA*FLOAT(IDINT(DABS(X(J))/SDELTA+.5D0)),X(J))
  20  CONTINUE
      WRITE (IOUTD,9999) NBIT1
9999  FORMAT (38H1THE INITIAL COEFFICIENTS, ROUNDED TO , I3, 8H BITS, A,
     *    2HRE)
      DO 30 J=1,NSECT
        WRITE (IOUTD,9998) J, X(4*J-3), X(4*J-2), X(4*J-1), X(4*J)
  30  CONTINUE
9998  FORMAT (8H SECTION, I3/1H , 2D20.10/1H , 2D20.10)
C
C SET PARAMETER VALUES FOR H-AND-J; THESE MAY BE CHANGED IF DESIRED
C H-AND-J MULTIPLIES THE STEP SIZE BY RHO; TAKE RHO = .5
C
      RHO = .5D0
C
C H-AND-J STOPS IF THE FUNCTION VALUE FALLS BELOW EST; TAKE EST = 0.
C
      EST = 0.
C
C H-AND-J STOPS IF WE EXCEED "LIMIT" FUNCT CALLS; TAKE LIMIT = 10,000
C
      LIMIT = 10000
      WRITE (IOUTD,9997)
9997  FORMAT (48H0NEXT FOLLOWS A REPORT FROM THE SEARCH ALGORITHM)
      CALL HANDJ(N, DELTAZ, SDELTA, RHO, EST, LIMIT, X)
      WRITE (IOUTD,9996) NBIT1
9996  FORMAT (31H0THE FINAL COEFFICIENTS, HAVING, I3, 10H BITS, ARE)
      DO 40 J=1,NSECT
        WRITE (IOUTD,9998) J, X(4*J-3), X(4*J-2), X(4*J-1), X(4*J)
  40  CONTINUE
C
C PRINT OUT THE FINAL FINE GRID AND QUIT
C
      PRINT = .TRUE.
      CALL FUNCT(N, X, F)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:    INDATA
C THIS SUBROUTINE READS IN THE PROBLEM DATA:
C CARD 1 HAS THE NUMBER OF COEFFICIENTS PER STAGE ( 3 OR 4), THE
C INITIAL SETTING OF THE RANDOM NUMBER GENERATOR, WHETHER TWO-OPT
C SHOULD BE USED ( 0 IS NO, 1 IS YES), THE DESIRED BIT LENGTH OF THE
C COEFFICIENTS, AND THE INITIAL SEARCH DELTA ( IN BITS) FOR THE
C HOOKE AND JEEVES ALGORITHM.
C
C THE NEXT SET OF CARDS SPECIFIES THE GRID, ONE GRID POINT PER CARD:
C EACH CARD HAS FIRST THE FREQUENCY IN FRACTIONS OF THE NYQUIST
C FREQUENCY, NEXT THE DESIRED TRANSFER FUNCTION MAGNITUDE, AND THIRD
C THE TOLERANCE. TO INDICATE THE END OF THE GRID POINTS,
C SPECIFY A FREQUENCY OF 1.  THE M CARDS ARE COUNTED, M.LE.100 .
C THE PRESENT SUBROUTINE FUNCT ASSUMES FOR THE A CALCULATION THAT THE
C DESIRED TRANSFER FUNCTION TAKES ON ONLY THE VALUES 0 OR 1. THAT
C IS, THE FILTER HAS ONLY PASS AND STOP BANDS, AND EVERY PASS BAND
C HAS DESIRED VALUE 1, AND THE SAME TOLERANCE.
C
C THE NEXT CARD HAS THE NUMBER OF SECOND-ORDER SECTIONS, NSECT.
C
C THE NEXT CARDS HAVE THE "INFINITE PRECISION" COEFFICIENTS, ONE PER
C CARD. THE COEFFICIENTS FOR SECTION 1 ARE GIVEN FIRST, NEXT
C SECTION 2, ETC. WITHIN EACH SECTION THEY ARE IN THE ORDER: NUMERATOR
C COEFFICIENT OF Z**(-1), NUMERATOR COEFFICIENT OF Z**(-2) ( IF 4
C PER SECTION), DENOMINATOR COEFFICIENT OF Z**(-1), DENOMINATOR
C COEFFICIENT OF Z**(-2). THERE ARE 3*NSECT OR 4*NSECT COEFFICIENTS
C READ IN, DEPENDING ON WHETHER THERE ARE 3 OR 4 COEFFICIENTS PER
C SECTION, NSECT.LE.9 . IN THE CASE OF 3 COEFFICIENTS PER SECTION,
C THE SECOND NUMERATOR COEFFICIENT IS SET TO 1 AND FROZEN.
C-----------------------------------------------------------------------
C
      SUBROUTINE INDATA
      DOUBLE PRECISION W(100), Y(100), WEIGHT(100), X(36)
      LOGICAL PRINT, TWOPT, FIXCOF
      COMMON /ABC/ X, NSET, NTWOPT, NBIT1, NBIT2, NSECT, N
      COMMON /RAW/ W, Y, WEIGHT, M, KOUNT
      COMMON /STEER/ PRINT, TWOPT, FIXCOF
      COMMON /MACH/ IND, IOUTD
C
C INITIALIZE KOUNT, THE NUMBER OF FUNCTION EVALUATIONS
C
      KOUNT = 0
C
C READ AND WRITE CARD 1 PARAMETERS
C
      READ (IND,9999) NCOEF, NSET, NTWOPT, NBIT1, NBIT2
9999  FORMAT (5I3)
      WRITE (IOUTD,9998) NCOEF, NSET, NTWOPT, NBIT1, NBIT2
9998  FORMAT (24H1 ***** INPUT DATA *****/26H NUMBER OF COEFFICIENTS PE,
     *    10HR STAGE IS, I3/26H RANDOM INITIALIZATION IS , I3/7H TWOPT ,
     *    2HIS, I3/44H FINAL ( AND INITIAL ROUNDING ) PRECISION IS, I3,
     *    5H BITS/24H INITIAL SEARCH DELTA IS, I3, 5H BITS)
C
C SET LOGICAL VARIABLE TWOPT IF TWO-OPT IS CALLED FOR
C
      TWOPT = .FALSE.
      IF (NTWOPT.EQ.1) TWOPT = .TRUE.
C
C READ AND WRITE THE GRID
C
      WRITE (IOUTD,9997)
9997  FORMAT (20H0GRID SPECIFICATIONS/30H POINT NO.           FREQUENCY,
     *    40H           DESIRED Y           TOLERANCE)
      M = 0
  10  M = M + 1
      READ (IND,9996) W(M), Y(M), WEIGHT(M)
9996  FORMAT (3F10.5)
      WRITE (IOUTD,9995) M, W(M), Y(M), WEIGHT(M)
9995  FORMAT (1H , I9, 3D20.10)
      WEIGHT(M) = 1.D0/WEIGHT(M)
      IF (W(M).LT.1.D0) GO TO 10
C
C READ AND WRITE THE COEFFICIENTS
C
      READ (IND,9999) NSECT
      WRITE (IOUTD,9994) NSECT
9994  FORMAT (40H0INITIAL HIGH PRECISION COEFFICIENTS FOR, I3, 6H SECTI,
     *    7HONS ARE)
      N = 4*NSECT
      IF (NCOEF.EQ.3) GO TO 20
      FIXCOF = .FALSE.
      READ (IND,9993) (X(J),J=1,N)
9993  FORMAT (F20.16)
      GO TO 40
  20  FIXCOF = .TRUE.
      READ (IND,9993) (X(4*J-3),X(4*J-1),X(4*J),J=1,NSECT)
      DO 30 J=1,NSECT
        X(4*J-2) = 1.D0
  30  CONTINUE
  40  CONTINUE
      DO 50 J=1,NSECT
        WRITE (IOUTD,9992) J, X(4*J-3), X(4*J-2), X(4*J-1), X(4*J)
  50  CONTINUE
9992  FORMAT (8H SECTION, I3/1H , 2D20.10/1H , 2D20.10)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:    FUNCT
C THIS COMPUTES THE VALUE OF THE ERROR FUNCTION, USING THE
C GRID. IF PRINT.EQ. .TRUE., IT COMPUTES THE ERROR FUNCTION
C AT A FINER GRID ( 10 POINTS FOR EACH GRID POINT), AND PRINTS
C THE RESULTS; THIS IS USED TO PRINT OUT THE FINAL TRANSFER FUNCTION.
C
C NOTE THAT THE CALCULATION OF THE CONSTANT MULTIPLIER A ASSUMES
C THAT ALL THE PASSBANDS HAVE THE SAME DESIRED VALUE OF 1. THIS
C MUST BE CHANGED TO A BAND-BY-BAND CALCULATION IF THERE IS MORE
C THAN ONE PASSBAND DESIRED VALUE.
C-----------------------------------------------------------------------
C
      SUBROUTINE FUNCT(N, X, F)
      DOUBLE PRECISION W(100), Y(100), WEIGHT(100)
      DOUBLE PRECISION COST2(100), COS2T4(100)
      DOUBLE PRECISION SUB1(100), SUB2(100), SUB3(100)
      DOUBLE PRECISION SUB6(100), SUB7(100), SUB8(100)
      DOUBLE PRECISION X(36), YHT(100), E(100)
      DOUBLE PRECISION YHTHT(100,11)
      DOUBLE PRECISION A, F, PI, NUM, DEN, Y1, Y2, FREQ, CCC, SSS, A1,
     *    ERROR
      LOGICAL PRINT, TWOPT, FIXCOF
      COMMON /RAW/ W, Y, WEIGHT, M, KOUNT
      COMMON /STEER/ PRINT, TWOPT, FIXCOF
      COMMON /MACH/ IND, IOUTD
C
C CHECK THAT THE POLES OF THE FILTER ARE INSIDE OR ON THE
C UNIT CIRCLE; NO CHECKS ARE MADE ON NUMERATOR COEFFICIENTS
C
      K = N/4
      DO 10 J=1,K
        J4 = J*4
        IF (X(J4).LE.1.D0.AND.1.D0+X(J4).GE.DABS(X(J4-1))) GO TO 10
        F = 1.D12
        RETURN
  10  CONTINUE
      PI = 4.D0*DATAN(1.D0)
      IF (KOUNT.NE.0) GO TO 30
C
C COMPUTE AND SAVE THE TRIG FUNCTIONS AT THE GRID POINTS
C THE FIRST TIME THROUGH
C
      DO 20 I=1,M
        COST2(I) = 2.D0*DCOS(PI*W(I))
        COS2T4(I) = COST2(I)**2
  20  CONTINUE
C
C MOVE THE CALCULATION OF CONSTANTS OUTSIDE THE INNER LOOP
C
  30  DO 40 J=1,K
        J4 = J*4
        SUB1(J) = X(J4-3)**2 + (X(J4-2)-1.D0)**2
        SUB2(J) = X(J4-3)*(X(J4-2)+1.D0)
        SUB3(J) = X(J4-2)
        SUB6(J) = X(J4-1)**2 + (X(J4)-1.D0)**2
        SUB7(J) = X(J4-1)*(X(J4)+1.D0)
        SUB8(J) = X(J4)
  40  CONTINUE
C
C EVALUATE THE MAGNITUDE OF THE TRANSFER FUNCTION, YHT, AT EACH
C OF THE M GRID POINTS
C
      DO 60 I=1,M
        NUM = 1.D0
        DEN = 1.D0
        DO 50 J=1,K
          NUM = NUM*(SUB1(J)+SUB2(J)*COST2(I)+SUB3(J)*COS2T4(I))
          DEN = DEN*(SUB6(J)+SUB7(J)*COST2(I)+SUB8(J)*COS2T4(I))
  50    CONTINUE
        YHT(I) = DSQRT(NUM/DEN)
  60  CONTINUE
C
C FIND THE LARGEST AND SMALLEST VALUES OF YHT IN THE PASSBANDS
C
      Y1 = 1.D12
      Y2 = 0.D0
      DO 70 I=1,M
        IF (Y(I).EQ.0.D0) GO TO 70
        Y1 = DMIN1(Y1,YHT(I))
        Y2 = DMAX1(Y2,YHT(I))
  70  CONTINUE
C
C DEFINE THE CONSTANT MULTIPLIER A TO BE THE RECIPROCAL OF THE AVERAGE
C OF THE SMALLEST AND LARGEST YHT"S IN THE PASSBANDS. THIS MAKES SENSE
C IF ALL THE PASSBANDS HAVE THE SAME DESIRED VALUE OF 1 AND THE SAME
C TOLERANCE; OTHERWISE THIS CALCULATION OF A MUST BE MODIFIED.
C
      A = 2.D0/(Y1+Y2)
C
C CALCULATE THE WEIGHTED ERROR AT THE GRID POINTS
C
      F = 0.D0
      DO 80 I=1,M
        YHT(I) = A*YHT(I)
        E(I) = WEIGHT(I)*DABS(YHT(I)-Y(I))
        F = DMAX1(F,E(I))
  80  CONTINUE
      KOUNT = KOUNT + 1
      IF (.NOT.PRINT) RETURN
C
C IF REQUIRED, PRINT OUT THE RESULTS ON THE COARSE GRID
C
      WRITE (IOUTD,9999) A, F
9999  FORMAT (26H1***** FINAL RESULTS *****/18H THE CONSTANT A IS,
     *    D20.10/33H FUNCTION VALUE ON COARSE GRID IS, D20.10)
      WRITE (IOUTD,9998)
9998  FORMAT (1H , 7X, 9HFREQUENCY, 7X, 9HDESIRED Y, 8X, 8HACTUAL Y,
     *    11X, 5HERROR)
      DO 90 I=1,M
        WRITE (IOUTD,9997) W(I), Y(I), YHT(I), E(I)
  90  CONTINUE
9997  FORMAT (1H , 4D16.8)
C
C TO RE-CALCULATE A, CALCULATE N/D=YHTHT ON THE FINE GRID
C THE SAME RESERVATION ABOUT THIS CALCULATION OF A APPLIES AS ABOVE
C
      Y1 = 1.D12
      Y2 = 0.D0
      MM = M - 1
      DO 130 I=1,MM
C
C SKIP TRANSITION BANDS, DEFINED AS THOSE WHERE Y(I) CHANGES
C
        IF (Y(I).NE.Y(I+1)) GO TO 130
        DO 120 II=1,11
          IF (II.NE.11) GO TO 100
C
C FIND CASES WHERE THE LAST POINT OF THE BAND WILL NOT BE COMPUTED LATER
C
          IF (I.EQ.MM) GO TO 100
          IF (Y(I+1).NE.Y(I+2)) GO TO 100
          GO TO 120
 100      FREQ = W(I) + (W(I+1)-W(I))*.1D0*FLOAT(II-1)
          CCC = 2.D0*DCOS(PI*FREQ)
          SSS = CCC**2
          NUM = 1.D0
          DEN = 1.D0
          DO 110 J=1,K
            J4 = J*4
            NUM = NUM*(X(J4-3)**2+(X(J4-2)-1.D0)**2+X(J4-3)*(X(J4-2)
     *          +1.D0)*CCC+X(J4-2)*SSS)
            DEN = DEN*(X(J4-1)**2+(X(J4)-1.D0)**2+X(J4-1)*(X(J4)+1.D0)*
     *          CCC+X(J4)*SSS)
 110      CONTINUE
          YHTHT(I,II) = DSQRT(NUM/DEN)
          IF (Y(I).EQ.0.D0) GO TO 120
          Y1 = DMIN1(Y1,YHTHT(I,II))
          Y2 = DMAX1(Y2,YHTHT(I,II))
 120    CONTINUE
 130  CONTINUE
      A = 2.D0/(Y1+Y2)
C
C NOW REPEAT THE SAME LOOPS, FINDING THE ERROR AND PRINTING THE RESULTS
C
      WRITE (IOUTD,9996)
9996  FORMAT (53H1FINAL RESULTS ON FINE GRID- 10 POINTS PER GRID POINT)
      WRITE (IOUTD,9998)
      F = 0.D0
      DO 160 I=1,MM
        IF (Y(I).NE.Y(I+1)) GO TO 160
        IP = I + 1
        WRITE (IOUTD,9995) I, IP
9995    FORMAT (27H ***** BAND FROM GRID POINT, I5, 3H TO, I5, 6H *****)
        DO 150 II=1,11
          IF (II.NE.11) GO TO 140
          IF (I.EQ.MM) GO TO 140
          IF (Y(I+1).NE.Y(I+2)) GO TO 140
          GO TO 150
 140      FREQ = W(I) + (W(I+1)-W(I))*.1D0*FLOAT(II-1)
          A1 = A*YHTHT(I,II)
          ERROR = WEIGHT(I)*DABS(A1-Y(I))
          F = DMAX1(F,ERROR)
          WRITE (IOUTD,9997) FREQ, Y(I), A1, ERROR
 150    CONTINUE
 160  CONTINUE
      WRITE (IOUTD,9994) A, F
9994  FORMAT (37H0THE VALUE OF A FROM THE FINE GRID IS, D20.10/6H0THE F,
     *    43HINAL VALUE OF THE ERROR ON THE FINE GRID IS, D20.10)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:    HANDJ
C THIS IS AN IMPLEMENTATION OF A RANDOMIZED HOOKE AND JEEVES
C SEARCH ALGORITHM. IT FINDS LOCAL OPTIMA; THE COORDINATES ARE
C RANDOMLY ORDERED EACH TIME EXPLOR IS INVOKED. FLOW CHARTS FOR
C THIS SUBROUTINE AND THE NEXT ARE INCLUDED IN THE REFERENCED PAPER.
C-----------------------------------------------------------------------
C
      SUBROUTINE HANDJ(N, DELTAZ, SDELTA, RHO, FMIN, LIMIT, PSI)
      DOUBLE PRECISION PSI(36), THETA(36), PHI(36)
      DOUBLE PRECISION W(100), Y(100), WEIGHT(100)
      DOUBLE PRECISION FPSI, FPHI, DELTA, RHO, DELTAZ, SDELTA, FMIN
      COMMON /RAW/ W, Y, WEIGHT, M, KOUNT
      COMMON /MACH/ IND, IOUTD
      CALL FUNCT(N, PSI, FPSI)
      DELTA = DELTAZ
      WRITE (IOUTD,9999)
9999  FORMAT (10H NO. CALLS, 10H PAT. SIZE, 11X, 5HDELTA, 7X, 7HOLD VAL,
     *    2HUE, 7X, 9HNEW VALUE)
  10  IF (DELTA.LT.SDELTA .OR. FPSI.LT.FMIN .OR. KOUNT.GT.LIMIT) GO TO
     *    50
      CALL SET(PHI, PSI, N)
      FPHI = FPSI
      CALL EXPLOR(PHI, FPHI, N, DELTA)
      I = 0
  20  IF (FPHI.GE.FPSI .OR. FPSI.LT.FMIN .OR. KOUNT.GT.LIMIT) GO TO 40
      I = I + 1
      WRITE (IOUTD,9998) KOUNT, I, DELTA, FPSI, FPHI
9998  FORMAT (4X, I6, 4X, I6, 3D16.8)
      CALL SET(THETA, PSI, N)
      CALL SET(PSI, PHI, N)
      FPSI = FPHI
      DO 30 J=1,N
        PHI(J) = 2.*PHI(J) - THETA(J)
  30  CONTINUE
      CALL FUNCT(N, PHI, FPHI)
      CALL EXPLOR(PHI, FPHI, N, DELTA)
      GO TO 20
  40  IF (I.GT.0) GO TO 10
      DELTA = RHO*DELTA
      GO TO 10
  50  WRITE (IOUTD,9997) KOUNT, FPSI
9997  FORMAT (34H0FINAL VALUES- NO. OF FUNCT CALLS=, I6, 7H FUNCT=,
     *    D20.10)
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:    EXPLOR
C THIS IS A LOCAL EXPLORATION SUBROUTINE USED BY H-AND-J
C-----------------------------------------------------------------------
C
      SUBROUTINE EXPLOR(PHI, FPHI, N, DELTA)
      DOUBLE PRECISION PHI(36), FPHI, DELTA, SAVE, FNEW, SAVEI, SAVEJ
      INTEGER LIST(36)
      COMMON /STEER/ PRINT, TWOPT, FIXCOF
      LOGICAL PRINT, TWOPT, FIXCOF
C
C RANDOMIZE THE ORDER OF THE COORDINATES
C
      CALL SHUFF(N, LIST)
      ICOUNT = 1
  10  IF (ICOUNT.GT.N) RETURN
      I = LIST(ICOUNT)
C
C IF ONLY THREE COEFFICIENTS WERE READ IN PER STAGE
C THEN DO NOT CHANGE THE SECOND NUMERATOR COEFFICIENT
C
      IF (FIXCOF .AND. I-(I/4)*4.EQ.2) GO TO 100
      SAVE = PHI(I)
      PHI(I) = PHI(I) + DELTA
      CALL FUNCT(N, PHI, FNEW)
      IF (FNEW.GE.FPHI) GO TO 20
      FPHI = FNEW
      GO TO 100
  20  PHI(I) = PHI(I) - 2.*DELTA
      CALL FUNCT(N, PHI, FNEW)
      IF (FNEW.GE.FPHI) GO TO 30
      FPHI = FNEW
      GO TO 100
  30  PHI(I) = SAVE
      IF (.NOT.TWOPT) GO TO 100
C
C BEGINNING OF 2-OPT
C
      JCOUNT = ICOUNT + 1
  40  IF (JCOUNT.GT.N) GO TO 100
      J = LIST(JCOUNT)
C
C SIMILARLY FOR TWOPT, IF THREE COEFFICIENTS PER STAGE WERE
C READ IN, DO NOT CHANGE THE SECOND NUMERATOR COEFFICIENT
C
      IF (FIXCOF .AND. J-(J/4)*4.EQ.2) GO TO 90
      SAVEI = PHI(I)
      SAVEJ = PHI(J)
      PHI(I) = PHI(I) + DELTA
      PHI(J) = PHI(J) + DELTA
      CALL FUNCT(N, PHI, FNEW)
      IF (FNEW.GE.FPHI) GO TO 50
      FPHI = FNEW
      GO TO 90
  50  PHI(I) = PHI(I) - 2.*DELTA
      CALL FUNCT(N, PHI, FNEW)
      IF (FNEW.GE.FPHI) GO TO 60
      FPHI = FNEW
      GO TO 90
  60  PHI(J) = PHI(J) - 2.*DELTA
      CALL FUNCT(N, PHI, FNEW)
      IF (FNEW.GE.FPHI) GO TO 70
      FPHI = FNEW
      GO TO 90
  70  PHI(I) = PHI(I) + 2.*DELTA
      CALL FUNCT(N, PHI, FNEW)
      IF (FNEW.GE.FPHI) GO TO 80
      FPHI = FNEW
      GO TO 90
  80  PHI(I) = SAVEI
      PHI(J) = SAVEJ
  90  JCOUNT = JCOUNT + 1
      GO TO 40
C
C END OF 2-OPT
C
 100  ICOUNT = ICOUNT + 1
      GO TO 10
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:    SET
C THIS SETS VECTOR A EQUAL TO VECTOR B
C-----------------------------------------------------------------------
C
      SUBROUTINE SET(A, B, N)
      DOUBLE PRECISION A(36), B(36)
      DO 10 I=1,N
        A(I) = B(I)
  10  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:    SHUFF
C THIS RANDOMLY ORDERS LIST
C-----------------------------------------------------------------------
C
      SUBROUTINE SHUFF(N, LIST)
      INTEGER LIST(36)
      DO 10 I=1,N
        LIST(I) = I
  10  CONTINUE
      DO 20 LL=1,N
        L = N - LL + 1
        J = INT(FLOAT(L)*UNI(0.)) + 1
        K = LIST(L)
        LIST(L) = LIST(J)
        LIST(J) = K
  20  CONTINUE
      RETURN
      END