SUBROUTINE APPROX(M, N, K, X, Y, EPSH, REF, MAXIT, HMAX, H,       APP   10
     * A, EQUAL)
C THIS SUBROUTINE, IN CONJUNCTION WITH SUBROUTINE EXCH,
C COMPUTES THE MINIMAX (CHEBYCHEV) POLYNOMIAL WHICH FITS
C THE DATA IN X AND Y.
C THIS SUBPROGRAM IS A TRANSLATION FROM ALGOL OF
C H. SCHMITT*S ALGORITHM NUMBER 409 IN CACM, V14,
C PP. 355-356(1971).
C THE PRIMARY REFERENCE IS STIEFEL, E. L.  NUMERICAL
C METHODS OF CHEBYCHEV APPROXIMATION, IN *ON NUMERICAL
C APPROXIMATION*, R. LANGER, (EDITOR), UNIVERSITY OF
C WISCONSIN PRESS, 1958, PP.217-232.
C TRANSLATION BY JOSEPH C. SIMPSON, MARQUETTE UNIVERSITY.
C THE EXCHANGE ALGORITHM IS A FINITE ITERITIVE PROCESS
C REQUIRING, AT MOST, BICO(N,P) ITERATIONS, WHERE P =
C NUMBER OF POINTS FIT IN ONE ITERATION (SEE EXCH FOR THE
C VALUE OF P).  SINCE BICO(N,P) CAN BE VERY LARGE, IT IS
C POSSIBLE THAT THE ROUTINE WILL NOT CONVERGE WITHIN MAXIT
C ITERATIONS.  THE OTHER POSSIBILITY OF FAILURE OCCURS
C WHEN INSUFFICIENT FLOATING POINT PRECISION IS AVAILABLE
C FOR THE INPUT DATA CHOSEN.
C              DESCRIPTION OF INPUT-OUTPUT VARIABLES.
C M = DEGREE OF FITTING POLYNOMIAL.
C N = NUMBER OF DATA POINTS STORED IN X AND Y.
C X = INDEPENDENT VARIABLE VALUES STORED IN ASCENDING ORDER.
C Y = DEPENDENT VARIABLE VALUES CORRESPONDING TO ABOVE DATA.
C K CONTROLS THE CHARACTER OF THE POLYNOMIAL.
C K = 0 IMPLIES MIXED PARITY POLYNOMIAL.  X(1) MAY BE ANY
C FLOATING POINT VALUE.
C K = 1 IMPLIES ODD POLYNOMIAL.  X(1) MUST BE .GT.0.
C K = 2 IMPLIES EVEN POLYNOMIAL.  X(1) MUST BE .GE.0.
C EPSH.  ABS(EPSH) = TOLERENCE FOR LEVELING.  IF EPSH IS
C NEGATIVE, REF CONTAINS THE INITIAL POINT SET.  ABS(EPSH)
C SHOULD BE SIGNIFICANT WHEN COMPARED WITH ONE.  A USEFUL
C VALUE FOR 24 BIT MANTISSA IS EPSH=2.E-7.
C REF = AN ARRAY OF INITIAL POINTS IF EPSH.LT.0.  OTHERWISE
C IT IS NOT AN INPUT ARRAY.  UPON RETURN IT CONTAINS THE
C MAXIMUM DEVIATION POINTS.
C MAXIT = UPPER LIMIT FOR NUMBER OF EXCHANGE STEPS.
C HMAX CONTAINS THE MAXIMUM DEVIATION(OUTPUT).
C H CONTAINS THE APPROXIMATION ERRORS(OUTPUT).
C A CONTAINS THE POLYNOMIAL COEFFICIENTS IN ORDER OF
C INCREASING POWERS(OUTPUT). IF MAXIT IS EXCEEDED, REF
C CONTAINS A NEW REFERENCE SET.
C DIMENSION REF(REFMAX),DAA(M+1),D(M+2),Z(M+4),A(M+1),
C AA(M+1),C(M+2),X(N),Y(N),H(N)
C REFMAX(K=0)=M+2. REFMAX(K=1)=INT((M+3)/2).
C REFMAX(K=2)=2+M/2.
C INPUT VARIABLES ALTERED DURING EXECUTION CAN INCLUDE M,
C EPSH,REF.
C EQUAL RECORDS THE SUCCESS OR FAILURE OF THE ALGORITHM.
C EQUAL=1  SUCCESFUL.
C EQUAL=0  CONVERGENCE OF EXCHANGE PROCESS NOT ACHEIVED.
C EQUAL=-1  INPUT ERROR.
C EQUAL=-2  ALGORITHM FAILURE.  NOTE.. A MAY CONTAIN THE
C           COLLOCATION POLYNOMIAL.
C          DESCRIPTION OF SOME LOCAL VARIABLES.
C SMALLEST REPRESENTABLE NUMBER.  IT IS USED ONLY WHEN
C COLLOCATION IS ENCOUNTERED.
C THE FOLLOWING VARIABLES ARE DECLARED INTEGER, BUT ARE USED
C AS LOGICAL VARIABLES ONLY,  .TRUE.=1  AND  .FALSE.=0.
C K0,K1,B1,B2.
C Z(I+1) CONTAINS THE ITH POINT OF THE CURRENT REFERENCE
C SET.  DAA(I+1) CONTAINS THE CORRECTION TO THE COEFFICIENT
C OF X**(I).
C TO OBTAIN A DOUBLE PRECISION VERSION, CHANGE THE
C REAL DECLARATION TO DOUBLE PRECISION, ABS TO DABS, COS
C TO DCOS, AND FLOAT TO DFLOAT.
C THIS DIMENSION STATEMENT ALLOWS 50 POINTS AND A MIXED
C ORDER POLYNOMIAL OF DEGREE 30.
      DIMENSION XX(50), AA(31), DAA(31), D(32), Z(34), C(32)
      DIMENSION X(1), Y(1), H(1), REF(1), A(1)
      DOUBLE PRECISION SD, QD, C, T, D, DT, AA, DAA, XX, MAX
      REAL EPSH, HMAX, Y, H, A, PI, Q, S, X1, XA, XE, COS
      INTEGER N, M, K, MAXIT, REF, I, J, P, Q1, Q2, R, K0, K1, B1,
     * B2, Z, ITEMP, MOLD, HIGH, JJ, II, LOW, EQUAL, NTAPE, B3
C MACHINE DEPENDENT CONSTANTS.
C NTAPE IS THE PRINTER UNIT NUMBER.
      NTAPE = 6
      PI = 3.1415926535
      MOLD = M
C ZERO THE COEFFICIENT ARRAY.
      IF (K-1) 10, 20, 30
   10 K0 = 1
      K1 = 0
      Q1 = 0
      Q2 = 1
      GO TO 50
   20 K0 = 0
      K1 = 1
      Q1 = 1
      Q2 = 2
      GO TO 40
   30 K0 = 0
      K1 = 0
      Q1 = 0
      Q2 = 2
C ADJUST M ACCORDING TO THE VALUE OF K
   40 M = (M-Q1)/2
   50 P = M + 2
C SCREEN INPUT PARAMETERS.
      IF (P-N) 60, 60, 780
   60 IF (M) 780, 70, 70
   70 IF (K-1) 100, 80, 90
   80 IF (X(1)) 790, 790, 100
   90 IF (X(1)) 790, 100, 100
  100 DO 120 I=2,N
        ITEMP = I - 1
        IF (X(I)-X(ITEMP)) 800, 110, 110
  110   CONTINUE
  120 CONTINUE
      ITEMP = MOLD + 1
      DO 130 I=1,ITEMP
        A(I) = 0.
  130 CONTINUE
      Z(1) = 0
      ITEMP = P + 2
      Z(ITEMP) = N + 1
      IF (EPSH) 140, 170, 170
C IF EPSH.GT.0, BRANCH TO Z SETUP SECTION.
C IF EPSH IS NEGATIVE, TEST REF AND LOAD IT INTO Z.
  140 EPSH = -EPSH
      J = 0
      DO 160 I=1,P
        ITEMP = I + 1
        R = REF(I)
        Z(ITEMP) = R
C BRANCH TO ERROR EXIT UNLESS REF IS MONOTONICALLY
C INCREASING.
        IF (J-R) 150, 810, 810
  150   J = R
  160 CONTINUE
      IF (J-N) 300, 300, 810
C BRANCH AROUND Z SETUP SECTION.
C THIS SECTION LOADS Z WITH THE POINTS CLOSEST TO THE
C CHEBYCHEV ABSCISSAS.
  170 X1 = X(1)
      XE = X(N)
      IF (K0) 190, 190, 180
C IF THE POLYNOMIAL IS NOT MIXED, I.E.HAS A DEFINITE PARITY,
C BRANCH TO 20.
  180 XA = XE + X1
      XE = XE - X1
      Q = PI/FLOAT(M+1)
      GO TO 200
  190 XA = 0.
      XE = XE + XE
      ITEMP = 2*(M+1) + Q1
      Q = PI/FLOAT(ITEMP)
C CALCULATE THE JTH CHEBYCHEV ABSCISSA AND LOAD Z(J+1) WITH
C THE APPROPRIATE INDEX FROM THE DATA ABSCISSAS.
  200 DO 270 JJ=1,P
        J = P + 1 - JJ
        X1 = XA + XE*(COS(Q*FLOAT(P-J)))
        ITEMP = J + 2
        R = Z(ITEMP)
        HIGH = R - 1
        IF (HIGH-2) 230, 210, 210
  210   DO 220 II=2,HIGH
          I = HIGH + 2 - II
          ITEMP = I - 1
          IF (X(I)+X(ITEMP)-X1) 240, 240, 220
C WHEN THE CHEBYCHEV ABSCISSA IS BRACKETED BY TWO INPUT
C ABSCISSAS, BRANCH TO 240.
  220   CONTINUE
  230   I = 1
  240   ITEMP = J + 1
        IF (I-R) 250, 260, 260
  250   Z(ITEMP) = I
        GO TO 270
  260   Z(ITEMP) = R - 1
  270 CONTINUE
C IF THE LOWER CHEBYCHEV ABSCISSAS ARE LESS THAN X(1), LOAD
C THE LOWER ELEMENTS OF Z WITH THE LOWEST POINTS.
      J = 0
  280 J = J + 1
      ITEMP = J + 1
      IF (Z(ITEMP)-J) 290, 300, 300
  290 Z(ITEMP) = J
      GO TO 280
C M1 ENTRY.  INITIALIZE VARIABLES TO PREPARE FOR EXCHANGE
C ITERATION.
  300 ITEMP = M + 1
C ZERO THE AA ARRAY.
      DO 310 I=1,ITEMP
        AA(I) = 0.
  310 CONTINUE
C LOAD H WITH THE ORDINATES AND XX(I) WITH THE ABSCISSAS IF
C THE POLYNOMIAL IS MIXED.  IF THE POLYNOMIAL IS EVEN OR ODD
C LOAD XX WITH THE SQUARES OF THE ABSCISSAS.
      DO 340 I=1,N
        H(I) = Y(I)
        Q = X(I)
        IF (K0) 320, 320, 330
  320   XX(I) = Q*Q
        GO TO 340
  330   XX(I) = Q
  340 CONTINUE
      B1 = 0
      B2 = 0
      B3 = 0
      R = -1
      T = 0.
C ITERATION ENTRY.  R IS THE ITERATION INDEX.
  350 R = R + 1
      S = 1.
C COMPUTATION OF DIVIDED DIFFERENCES SCHEME.
      IF (K1) 380, 380, 360
C IF THE POLYNOMIAL IS ODD, BRANCH TO 380.
  360 DO 370 I=1,P
        S = -S
        ITEMP = I + 1
        J = Z(ITEMP)
        Q = X(J)
        C(I) = (H(J)+S*T)/Q
        D(I) = S/Q
  370 CONTINUE
      GO TO 400
  380 DO 390 I=1,P
        S = -S
        ITEMP = I + 1
        ITEMP = Z(ITEMP)
        C(I) = H(ITEMP) + S*T
        D(I) = S
  390 CONTINUE
  400 CONTINUE
      DO 420 I=2,P
        DO 410 JJ=I,P
          J = P + I - JJ
          ITEMP = J + 1
          ITEMP = Z(ITEMP)
          QD = XX(ITEMP)
          ITEMP = 2 + J - I
          ITEMP = Z(ITEMP)
          QD = QD - XX(ITEMP)
          ITEMP = J - 1
          C(J) = (C(J)-C(ITEMP))/QD
          D(J) = (D(J)-D(ITEMP))/QD
  410   CONTINUE
  420 CONTINUE
      DT = -C(P)/D(P)
      T = T + DT
C COMPUTATION OF POLYNOMIAL COEFFICIENTS.
      HIGH = M + 1
      DO 450 II=1,HIGH
        I = HIGH - II
        ITEMP = I + 1
        DAA(ITEMP) = C(ITEMP) + DT*D(ITEMP)
        ITEMP = I + 2
        ITEMP = Z(ITEMP)
        QD = XX(ITEMP)
        LOW = I + 1
        IF (M-LOW) 450, 430, 430
  430   DO 440 J=LOW,M
          JJ = J + 1
          DAA(J) = DAA(J) - QD*DAA(JJ)
  440   CONTINUE
  450 CONTINUE
      DO 460 I=1,HIGH
        AA(I) = AA(I) + DAA(I)
  460 CONTINUE
C EVALUATION OF THE POLYNOMIAL TO GET THE APPROXIMATION
C ERRORS.
      MAX = 0.
      DO 540 I=1,N
        SD = AA(HIGH)
        QD = XX(I)
        IF (-M) 470, 490, 490
  470   DO 480 JJ=1,M
          J = M - JJ + 1
          SD = SD*QD + AA(J)
  480   CONTINUE
  490   CONTINUE
        IF (K1) 510, 510, 500
C IF THE POLYNOMIAL IS  ODD MULTIPLY SD BY X(I).
  500   SD = SD*X(I)
  510   QD = Y(I) - SD
        H(I) = QD
        IF (DABS(QD)-MAX) 530, 530, 520
C LOAD MAX WITH THE LARGEST MAGNITUDE OF THE APPROXIMATION
C ARRAY.
  520   MAX = DABS(QD)
  530   CONTINUE
  540 CONTINUE
C TEST FOR ALTERNATING SIGNS.
      ITEMP = Z(2)
      IF (H(ITEMP)) 550, 560, 600
  550 J = 1
      GO TO 610
  560 J = 0
C THIS REPRESENTS A CASE WHERE THE POLYNOMIAL EXACTLY
C PREDICTS A DATA POINT.
      WRITE (NTAPE,99998)
      IF (B3) 570, 570, 690
  570 B3 = 1
      IF (EPSH-MAX) 580, 710, 710
  580 WRITE (NTAPE,99999)
      LOW = (N+1)/2 - (P+1)/2 + 1
      HIGH = LOW + P
      DO 590 I=LOW,HIGH
        Z(I) = I - 1
  590 CONTINUE
      GO TO 350
  600 J = -1
  610 CONTINUE
      DO 670 I=2,P
        ITEMP = I + 1
        ITEMP = Z(ITEMP)
C IF(H(ITEMP)) 339,340,341
        IF (H(ITEMP)) 620, 560, 630
  620   JJ = -1
        GO TO 640
  630   JJ = 1
  640   IF (J-JJ) 650, 660, 650
C ERROR ENTRY.  INSUFFICIENT ACCURACY FOR CALCULATION.
  650   B1 = 1
        GO TO 720
  660   J = -J
  670 CONTINUE
C SEARCH FOR ANOTHER REFERENCE.
      CALL EXCH(N, P, H, EPSH, Z, EQUAL)
      IF (EQUAL) 680, 680, 720
  680 IF (R-MAXIT) 350, 700, 700
  690 B3 = -1
      GO TO 720
  700 B2 = 1
      GO TO 720
  710 WRITE (NTAPE,99988) MAX
C END OF ITERATION LOOP.
C M2 ENTRY.  LOAD OUTPUT VARIABLES AND RETURN.
  720 HIGH = M + 1
C LOAD THE COEFFICIENTS INTO THE A ARRAY.
      DO 730 I=1,HIGH
        ITEMP = Q1 + Q2*(I-1) + 1
        A(ITEMP) = AA(I)
  730 CONTINUE
C LOAD REF WITH THE FINAL REFERENCE POINTS.
      DO 740 I=1,P
        ITEMP = I + 1
        REF(I) = Z(ITEMP)
  740 CONTINUE
      HMAX = MAX
      MAXIT = R
      IF (B3) 840, 750, 750
  750 IF (B1) 760, 760, 820
  760 IF (B2) 770, 770, 830
  770 M = MOLD
C NORMAL EXIT
      RETURN
C ERROR EXITS.
  780 WRITE (NTAPE,99990)
      WRITE (NTAPE,99996) N, MOLD, K
      EQUAL = -1
      GO TO 770
  790 WRITE (NTAPE,99990)
      WRITE (NTAPE,99995) K, X(1)
      EQUAL = -1
      GO TO 770
  800 WRITE (NTAPE,99990)
      WRITE (NTAPE,99994) (X(I),I=1,N)
      EQUAL = -1
      GO TO 770
  810 WRITE (NTAPE,99990)
      WRITE (NTAPE,99993) (REF(I),I=1,N)
      EQUAL = -1
      GO TO 770
  820 WRITE (NTAPE,99997)
      WRITE (NTAPE,99992) HMAX
      EQUAL = -2
      GO TO 770
  830 WRITE (NTAPE,99997)
      ITEMP = M + 2
      WRITE (NTAPE,99991) MAXIT, (REF(I),I=1,ITEMP)
      EQUAL = 0
      GO TO 770
  840 WRITE (NTAPE,99997)
      WRITE (NTAPE,99989)
      EQUAL = -2
      GO TO 770
99999 FORMAT (1H+, 26X, 35HONE MORE ATTEMPT WILL BE MADE USING,
     * 19H THE MIDDLE POINTS.)
99998 FORMAT (25H COLLOCATION HAS OCCURED.)
99997 FORMAT (26H ALGORITHM APPROX FAILURE.)
99996 FORMAT (10X, 28HINSUFFICIENT INFORMATION. N=, I3, 3H M=, I3,
     * 3H K=, I3)
99995 FORMAT (10X, 37HPOLYNOMIAL PARITY - ABCISSA CONFLICT.,
     * 6H  K = , I3, 7H  X(1)=, F20.10)
99994 FORMAT (10X, 40HABCISSAE OUT OF ORDER.  X ARRAY EQUALS
     * /(7E17.7))
99993 FORMAT (10X, 40HINITIAL POINT ARRAY NOT MONOTONIC INCREA,
     * 23HSING.  REF ARRAY EQUALS/(13I9))
99992 FORMAT (1X, 39HAPPROXIMATION ERRORS AT POINTS OF REFER,
     * 30HENCE DO NOT ALTERNATE IN SIGN./23H SUSPECT INSUFFICIENT W,
     * 11HORD LENGTH., 16H MAXIMUM ERROR= , E15.7/14H IF POLYNOMIAL,
     * 34H WAS ASSUMED TO BE OF MIXED PARITY, 18H, CHECK FOR EVEN O,
     * 13HR ODD PARITY.)
99991 FORMAT (31H MAXIMUM NUMBER OF ITERATIONS, , I6, 10H THE REFER,
     * 23HENCE ARRAY OBTAINED IS /(13I9))
99990 FORMAT (31H INPUT ERROR SUBPROGRAM APPROX.)
99989 FORMAT (1H+, 26X, 18HAPPROX TERMINATES.)
99988 FORMAT (1H+, 26X, 6H MAX= , E15.7, 23H .LT.HMIN.  NORMAL EXIT,
     * 1H.)
      END
      SUBROUTINE EXCH(N, P, H, EPSH, Z, EQUAL)                          EXC   10
C DIMENSION H(N),Z(P+2)
C N = NUMBER OF POINTS.
C P = NUMBER OF REFERENCE POINTS.
C EPSH = APPROXIMATION ERROR STANDARD.
C Z= INPUT   OLD REFERENCE INDICES.
C    OUTPUT  NEW REFERENCE INDICES.
C Z(1) SHOULD CONTAIN ZERO.
C Z(P+2) SHOULD CONTAIN (N+1).
C Z(I+1) CONTAINS THE ITH REFERENCE POINT INDEX.
C H = ARRAY OF APPROXIMATION ERRORS.
C EQUAL = 0 IMPLIES NORMAL EXCHANGE.
C EQUAL = 1 IMPLIES OLD AND NEW REFERENCE SETS ARE EQUAL.
C THE APPROXIMATION ERRORS ARE COMPARED RELATIVE TO EPSH.
C          DESCRIPTION OF SOME LOCAL VARIABLES.
C HZ1 CONTAINS THE LOW POINT APPROXIMATION ERROR.
C HZP CONTAINS THE HIGH POINT APPROXIMATION ERROR.
C MAXR CONTAINS THE LARGEST APPROXIMATION ERROR ABOVE THE
C REFERENCE POINT SET. INDR IS THE POINT INDEX FOR MAXR.
C MAXL CONTAINS THE LARGEST APPROXIMATION ERROR BELOW THE
C REFERENCE POINT SET. INDL IS THE POINT INDEX FOR MAXL.
C ITEMP IS A FIXED POINT VARIABLE USED FOR TEMPORARY STORAGE
C CODING FOR HIGHER LEVEL FORTRAN COULD ELIMINATE ITEMP.
C TO OBTAIN A DOUBLE PRECISION VERSION, CHANGE THE REAL
C DECLARATION TO DOUBLE PRECISION AND ABS TO DABS.
      DIMENSION Z(1), H(1)
      INTEGER Z, N, P, I, J, L, INDEX, INDL, INDR, ZE, ITEMP, LOW,
     * HIGH, EQUAL, II
      REAL EPSH, H, HZ1, HZP, MAX, MAXL, MAXR, SIG
      EQUAL = 0
      L = 0
      ITEMP = Z(2)
C ARBITRARILY CHOSEN EQUAL TO THE SIGN OF THE INPUT POINT.
C THIS WILL BE ADJUSTED LATER IF NECESSARY.
      IF (H(ITEMP)) 10, 10, 20
   10 SIG = 1.
      GO TO 30
   20 SIG = -1.
   30 DO 90 I=1,P
C DO 90 PRESCANS Z TO INSURE IT IS A PROPER CHOICE, I. E.
C RESETS Z IF NECESSARY SO THAT MAXIMUM ERROR POINTS ARE
C CHOSEN, GIVEN THE SIGN CONVENTION MENTIONED ABOVE.
C IN ORDER TO WORK PROPERLY, THIS SECTION REQUIRES Z(1)=0
C AND Z(P+2)=N+1.
        MAX = 0.
        SIG = -SIG
        ITEMP = I + 2
        ZE = Z(ITEMP) - 1
        LOW = Z(I) + 1
C SCAN THE OPEN POINT INTERVAL CONTAINING ONLY THE ITH
C INITIAL REFERENCE POINT. IN THE INTERVAL PICK THE
C POINT WITH LARGEST MAGNITUDE AND CORRECT SIGN.  MOST OF
C THE SORTING OCCURS IN THIS SECTION.  SIG CONTAINS THE
C SIGN ASSUMED FOR H(I).
        DO 60 J=LOW,ZE
          IF (SIG*(H(J)-MAX)) 50, 50, 40
   40     MAX = H(J)
          INDEX = J
   50     CONTINUE
   60   CONTINUE
        ITEMP = I + 1
        ITEMP = Z(ITEMP)
        MAXL = ABS(MAX)
        IF (ABS(MAX-H(ITEMP))/MAXL-EPSH) 80, 80, 70
C IF THE MAX ERROR IS SIGNIFICANTLY GREATER THAN THE INPUT
C POINT SWITCH TO THIS POINT.
   70   ITEMP = I + 1
        Z(ITEMP) = INDEX
        L = 1
   80   CONTINUE
   90 CONTINUE
      MAXL = 0.
      MAXR = 0.
      ITEMP = P + 1
      LOW = Z(ITEMP) + 1
      IF (LOW-N) 100, 100, 140
  100 CONTINUE
C FIND THE ERROR WITH LARGEST ABSOLUTE VALUE AND PROPER SIGN
C FROM AMONG THE POINTS ABOVE THE LAST REFERENCE POINT.
C THIS SECTION IS NECESSARY BECAUSE THE SET OF POINTS CHOSEN
C MAY BEGIN WITH THE WRONG SIGN ALTERNATION.
      DO 130 J=LOW,N
        IF (SIG*(MAXR-H(J))) 120, 120, 110
  110   MAXR = H(J)
        INDR = J
  120   CONTINUE
  130 CONTINUE
  140 CONTINUE
C FIND THE ERROR WITH LARGEST ABSOLUTE VALUE AND PROPER SIGN
C FROM AMONG THE POINTS BELOW THE FIRST REFERENCE POINT.
C THIS SECTION IS NECESSARY BECAUSE THE SET OF POINTS CHOSEN
C MAY BEGIN WITH THE WRONG SIGN ALTERNATION.
      ITEMP = Z(2)
      HZ1 = H(ITEMP)
      HIGH = ITEMP - 1
      IF (HIGH) 230, 230, 150
  150 CONTINUE
      IF (HZ1) 160, 170, 180
  160 SIG = -1.
      GO TO 190
  170 SIG = 0.
      GO TO 190
  180 SIG = 1.
  190 CONTINUE
      DO 220 J=1,HIGH
        IF (SIG*(MAXL-H(J))) 210, 210, 200
  200   MAXL = H(J)
        INDL = J
  210   CONTINUE
  220 CONTINUE
  230 CONTINUE
      MAXL = ABS(MAXL)
      MAXR = ABS(MAXR)
      HZ1 = ABS(HZ1)
      ITEMP = P + 1
      ITEMP = Z(ITEMP)
      HZP = ABS(H(ITEMP))
C MAXL AND MAXR CONTAIN THE MAGNITUDE OF THE SIGNIFICANT
C ERRORS OUTSIDE THE REFERENCE POINT SET.  IF EITHER IS
C ZERO, THE REFERENCE POINT SET EXTENDS TO THE END POINT
C ON THAT SIDE OF THE INTERVAL.
      IF (L) 290, 240, 290
C L=0 IMPLIES THAT THE PRESCAN OF DO 90 DID NOT CHANGE ANY
C POINTS.IF L=0 AND MAXL AND MAXR ARE NOT SIGNIFICANT WHEN
C COMPARED WITH UPPER AND LOWER REFERENCE POINT ERRORS,
C RESPECTIVELY, USE THE EQUAL EXIT.
  240 IF (MAXL) 250, 260, 250
  250 IF (EPSH-(MAXL-HZP)/MAXL) 290, 260, 260
  260 IF (MAXR) 270, 280, 270
  270 IF (EPSH-(MAXR-HZ1)/MAXR) 290, 280, 280
C EQUAL BRANCH POINT.
  280 EQUAL = 1
      RETURN
C ENTER HERE IF ANY CHANGES HAVE BEEN MADE.  THEN TEST TO
C SEE IF A POINT OUTSIDE THE PRESET POINT SET SHOULD BE
C INCLUDED.  IF NOT, RETURN.
  290 IF (MAXL) 320, 300, 320
  300 IF (MAXR) 320, 310, 320
C RETURN BRANCH (END).
  310 CONTINUE
      RETURN
C IF A POINT OUTSIDE THE PRESENT REFERENCE SET MUST BE
C INCLUDED, (I.E.  THE SIGN OF THE FIRST POINT, ASSUMED IN
C THE DO 90 SECTION IS WRONG) SHIFT TO THE SIDE OF LARGEST
C RELATIVE ERROR FIRST.
C CHECK THE OTHER SIDE.
  320 IF (MAXL-MAXR) 350, 350, 330
  330 IF (MAXL-HZP) 340, 340, 420
  340 IF (MAXR-HZ1) 310, 370, 370
  350 IF (MAXR-HZ1) 360, 360, 370
  360 IF (MAXL-HZP) 310, 420, 420
C SHR ENTRY.
C THIS SECTION INSERTS A POINT FROM ABOVE THE PRESCAN
C POINT SET.
  370 INDEX = Z(2)
      DO 380 I=2,P
C SHIFT POINT SET DOWN, DROPPING THE LOWEST POINT.
        ITEMP = I + 1
        Z(I) = Z(ITEMP)
  380 CONTINUE
      ITEMP = P + 1
C ADD THE NEW HIGH POINT.
      Z(ITEMP) = INDR
      IF (MAXL) 310, 310, 390
C IF MAXL=0, RETURN.
C IF MAXL.GT.0  REPLACE REFERENCE POINTS FROM THE LEFT,
C STOPPING THE FIRST TIME THE CANDIDATE FOR REPLACEMENT IS
C GREATER IN MAGNITUDE THAN THE PROSPECTIVE REPLACEE.
C ALTERNATION OF SIGNS IS PRESERVED BECAUSE NON-REPLACEMENT
C IMMEDIATLY TERMINATES THE PROCESS.
  390 DO 410 I=2,P
        ITEMP = Z(I)
        IF (ABS(H(INDL))-ABS(H(ITEMP))) 310, 400, 400
  400   J = Z(I)
        Z(I) = INDL
        INDL = INDEX
        INDEX = J
  410 CONTINUE
      GO TO 310
C ENTRY SHL.  THIS SECTION INSERTS A POINT FROM BELOW THE
C PRESCAN POINT SET.
  420 ITEMP = P + 1
      INDEX = Z(ITEMP)
C SHIFT REFERENCE POINT SET UP BY ONE.
      DO 430 II=2,P
        I = P + 2 - II
        ITEMP = I + 1
        Z(ITEMP) = Z(I)
  430 CONTINUE
C STORE LOWEST POINT IN Z(2)
      Z(2) = INDL
      IF (MAXR) 310, 310, 440
C IF MAXR=0, RETURN.
C IF MAXR.GT.0, START REPLACING REFERENCE POINTS FROM RIGHT,
C STOPPING THE FIRST TIME THE CANDIDATE FOR REPLACEMENT
C IS GREATER IN MAGNITUDE THAN THE PROSPECTIVE REPLACEE.
  440 DO 460 II=2,P
        I = P + 2 - II
        ITEMP = I + 1
        HIGH = Z(ITEMP)
        IF (ABS(H(INDR)-ABS(H(HIGH)))) 310, 450, 450
  450   J = Z(ITEMP)
        Z(ITEMP) = INDR
        INDR = INDEX
        INDEX = J
  460 CONTINUE
      GO TO 310
      END

.