C     ALGORITHM 525, COLLECTED ALGORITHMS FROM ACM. THIS WORK
C     PUBLISHED IN TRANS. MATH. SOFTWARE, 4(1), PP.82-94.
 
C
C  **** THIS PROGRAM RUNS ADAPT WITH DATA READ IN FROM CARDS
C       IT IS DESIGNED TO TEST ADAPT WITH THE 20 FUNCTIONS IN F(X).
C       IT READS DATA FOR EVERYTHING EXCEPT EDIST(=ATYPE) AND PRINTS
C       A HEADING WITH THE FUNCTIONS NAME. IT CALLS AND TIMES ADAPT,
C       THEN INDEPENDENTLY CHECKS THE ERROR IN ADAPT AND PLOTS THE
C       FUNCTION AND ERROR CURVE.
C
C NOTES -- USER MUST SUPPLY SYSTEM ROUTINES SECOND AND GRAPH
C          INTEGRATION USES ERRINT FROM ADAPT.
C          PRESENT DIMENSIONS (ON ZKNOTS + ZCOEFS) PREVENT USING
C             ERRINT TO CHECK ACCURACY FOR MORE THAN 100 KNOTS OR
C             POLYNOMIAL DEGREE 6.
C          QPOLY USED BY ERRINT IS A SINGLE ARGUMENT VERSION OF PPOLY.
C          SEE COMMENTS IN ADAPT ABOUT DIMENSION CONSTRAINTS AND
C             PORTABILITY.
C          IF ADAPT AND THIS DRIVER ARE CHANGED TO DOUBLE PRECISION
C             ONE PROBABLY WANTS TO LEAVE TIME,TSTART,TSTOP,XVALU,
C             GRAF1,GRAF2 AS SINGLE PRECISION
C             THEY ARE SEPARATELY DECLARED HERE.
C
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      DIMENSION XKNOTS(140), COEFS(140,12)
      REAL COEFS, ECHECK, F, XKNOTS
      REAL TIME, TSTART
      DIMENSION NAMES(25,10)
      EXTERNAL F
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
      COMMON /SAVEIT/ IKNOT
C
      COMMON /FDATA/ JFUNK, KOUNT, MAXK
      COMMON /INFO/ ECHECK, KFUNK, TIME, TSTART
      DATA NAMES(1,1), NAMES(1,2), NAMES(1,3), NAMES(1,4), NAMES(1,5),
     * NAMES(1,6), NAMES(1,7), NAMES(1,8), NAMES(1,9), NAMES(1,10) /
     * 4HSIMP,4HLE A,4HLGEB,4HRAIC,4H SIN,4HGULA,4HRITY,4H AT ,4H0.0 ,
     * 4H    /
      DATA NAMES(2,1), NAMES(2,2), NAMES(2,3), NAMES(2,4), NAMES(2,5),
     * NAMES(2,6), NAMES(2,7), NAMES(2,8), NAMES(2,9), NAMES(2,10) /
     * 4HINTE,4HRIOR,4H SIN,4HGULA,4HRITY,4H AT ,4HX = ,4H.307,4H7070,
     * 4H7707/
      DATA NAMES(3,1), NAMES(3,2), NAMES(3,3), NAMES(3,4), NAMES(3,5),
     * NAMES(3,6), NAMES(3,7), NAMES(3,8), NAMES(3,9), NAMES(3,10) /
     * 4HBROK,4HEN L,4HINE ,4H- BR,4HEAKS,4H AT ,4H1,2,,4H3.5 ,4HAND ,
     * 4H4   /
      DATA NAMES(4,1), NAMES(4,2), NAMES(4,3), NAMES(4,4), NAMES(4,5),
     * NAMES(4,6), NAMES(4,7), NAMES(4,8), NAMES(4,9), NAMES(4,10) /
     * 4HHUMP,4H AT ,4H2.5 ,4HRECI,4HPROC,4HAL O,4HF QU,4HARTI,4HC   ,
     * 4H    /
      DATA NAMES(5,1), NAMES(5,2), NAMES(5,3), NAMES(5,4), NAMES(5,5),
     * NAMES(5,6), NAMES(5,7), NAMES(5,8), NAMES(5,9), NAMES(5,10) /
     * 4HSIMP,4HLE C,4HUBIC,4H = (,4H4X-2,4H)**3,4H - (,4H4X-2,4H)   ,
     * 4H    /
      DATA NAMES(6,1), NAMES(6,2), NAMES(6,3), NAMES(6,4), NAMES(6,5),
     * NAMES(6,6), NAMES(6,7), NAMES(6,8), NAMES(6,9), NAMES(6,10) /
     * 4HVARI,4HABLE,4H OSC,4HILLA,4HTION,4H FRO,4HM SI,4HN(X*,4H*2-1,
     * 4H)   /
      DATA NAMES(7,1), NAMES(7,2), NAMES(7,3), NAMES(7,4), NAMES(7,5),
     * NAMES(7,6), NAMES(7,7), NAMES(7,8), NAMES(7,9), NAMES(7,10) /
     * 4HEXPO,4HNENT,4HIAL ,4H    ,4H    ,4H    ,4H    ,4H    ,4H    ,
     * 4H    /
      DATA NAMES(8,1), NAMES(8,2), NAMES(8,3), NAMES(8,4), NAMES(8,5),
     * NAMES(8,6), NAMES(8,7), NAMES(8,8), NAMES(8,9), NAMES(8,10) /
     * 4HFIFT,4HH DE,4HGREE,4H POL,4HY PL,4HUS F,4HINE ,4HSAW ,4HTOOT,
     * 4HH   /
      DATA NAMES(9,1), NAMES(9,2), NAMES(9,3), NAMES(9,4), NAMES(9,5),
     * NAMES(9,6), NAMES(9,7), NAMES(9,8), NAMES(9,9), NAMES(9,10) /
     * 4HSEVE,4HNTH ,4HDEGR,4HEE P,4HOLYN,4HOMIA,4HL   ,4H    ,4H    ,
     * 4H    /
      DATA NAMES(10,1), NAMES(10,2), NAMES(10,3), NAMES(10,4),
     * NAMES(10,5), NAMES(10,6), NAMES(10,7), NAMES(10,8), NAMES(10,9),
     * NAMES(10,10) /4HWILD,4H VAR,4HIATI,4HON B,4HETWE,4HEN .,4H47 A,
     * 4HND .,4H48  ,4H    /
      DATA NAMES(11,1), NAMES(11,2), NAMES(11,3), NAMES(11,4),
     * NAMES(11,5), NAMES(11,6), NAMES(11,7), NAMES(11,8), NAMES(11,9),
     * NAMES(11,10) /4HLOTS,4H OF ,4HOSCI,4HLLAT,4HION ,4H - F,4HIRST,
     * 4H EXA,4HMPLE,4H    /
      DATA NAMES(12,1), NAMES(12,2), NAMES(12,3), NAMES(12,4),
     * NAMES(12,5), NAMES(12,6), NAMES(12,7), NAMES(12,8), NAMES(12,9),
     * NAMES(12,10) /4HLOTS,4H OF ,4HOSCI,4HLLAT,4HION ,4H - S,4HECON,
     * 4HD EX,4HAMPL,4HE   /
      DATA NAMES(13,1), NAMES(13,2), NAMES(13,3), NAMES(13,4),
     * NAMES(13,5), NAMES(13,6), NAMES(13,7), NAMES(13,8), NAMES(13,9),
     * NAMES(13,10) /4HSING,4HULAR,4HITIE,4HS AT,4H BOT,4HH EN,4HD PO,
     * 4HINTS,4H    ,4H    /
      DATA NAMES(14,1), NAMES(14,2), NAMES(14,3), NAMES(14,4),
     * NAMES(14,5), NAMES(14,6), NAMES(14,7), NAMES(14,8), NAMES(14,9),
     * NAMES(14,10) /4HRAND,4HOM N,4HUMBE,4HRS A,4HDDED,4H TO ,4HSIN(,
     * 4HX)  ,4H    ,4H    /
      DATA NAMES(15,1), NAMES(15,2), NAMES(15,3), NAMES(15,4),
     * NAMES(15,5), NAMES(15,6), NAMES(15,7), NAMES(15,8), NAMES(15,9),
     * NAMES(15,10) /4HLARG,4HE FU,4HNCTI,4HON -,4H ADJ,4HUSTA,4HBLE ,
     * 4HSIZE,4H    ,4H    /
      DATA NAMES(16,1), NAMES(16,2), NAMES(16,3), NAMES(16,4),
     * NAMES(16,5), NAMES(16,6), NAMES(16,7), NAMES(16,8), NAMES(16,9),
     * NAMES(16,10) /4HSHIF,4HTED ,4HORIG,4HIN -,4H ADJ,4HUSTA,4HBLE ,
     * 4HSHIF,4HT   ,4H    /
      DATA NAMES(17,1), NAMES(17,2), NAMES(17,3), NAMES(17,4),
     * NAMES(17,5), NAMES(17,6), NAMES(17,7), NAMES(17,8), NAMES(17,9),
     * NAMES(17,10) /4HNUME,4HRICA,4HL DI,4HFFER,4HENTI,4HATIO,4HN OF,
     * 4H F(X,4H)   ,4H    /
      DATA NAMES(18,1), NAMES(18,2), NAMES(18,3), NAMES(18,4),
     * NAMES(18,5), NAMES(18,6), NAMES(18,7), NAMES(18,8), NAMES(18,9),
     * NAMES(18,10) /4HCOMP,4HLICA,4HTED ,4HFUNC,4HTION,4H WIT,4HH 7 ,
     * 4HPIEC,4HES  ,4H    /
      DATA NAMES(19,1), NAMES(19,2), NAMES(19,3), NAMES(19,4),
     * NAMES(19,5), NAMES(19,6), NAMES(19,7), NAMES(19,8), NAMES(19,9),
     * NAMES(19,10) /4HJUMP,4H DIS,4HCONT,4HINUI,4HTY A,4HT VA,4HRIAB,
     * 4HLE P,4HOINT,4H    /
      DATA NAMES(20,1), NAMES(20,2), NAMES(20,3), NAMES(20,4),
     * NAMES(20,5), NAMES(20,6), NAMES(20,7), NAMES(20,8), NAMES(20,9),
     * NAMES(20,10) /4HINFI,4HNITE,4H SIN,4HGULA,4HRITY,4H AT ,4HP10A,
     * 4H    ,4H    ,4H    /
C
   10 READ (5,99999) JFUNK, SMOOTH, DEGREE, LEVEL, NORM, A, B, ACCUR,
     * CHARF, NBREAK
      IF (JFUNK.EQ.0) STOP
      IF (CHARF.EQ.0.) CHARF = (B-A)*1.001
      IF (NBREAK.EQ.0) GO TO 20
      READ (5,99998) (DBREAK(K),XBREAK(K),BLEFT(K),BRIGHT(K),K=1,NBREAK)
   20 CONTINUE
      MAXK = 3200
      EDIST = 0
      WRITE (6,99997) JFUNK, (NAMES(JFUNK,K),K=1,10)
      KOUNT = 0
C
C     *****  CALL ON CLOCK TO INITIALIZE  TSTART  IN SECONDS  *****
C
      TSTART = 0.0
      CALL ADAPT(F, XKNOTS, COEFS, 140, 12)
      CALL SUMRY2(XKNOTS, COEFS, 140, 12)
      GO TO 10
99999 FORMAT (4I3, F4.3, 2F10.3, E12.4, F10.3, 2I3)
99998 FORMAT (2(I5, 3F11.5))
99997 FORMAT (//5(3H+++), 18HADAPT FOR FUNCTION, I3, 3H = , 10A4, 4X,
     * 5(3H+++))
      END
      BLOCK DATA
C         SET PARAMETERS FOR FUNCTIONS IN F
      REAL PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A, P6B, P7,
     * P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B, P10C
      COMMON /FPARS/ PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A,
     * P6B, P7, P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B,
     * P10C
      DATA PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B
     * /.02,.4,.6,.001,6.,8000.,2037./
      DATA P6A, P6B, P7, P8A, P8B, P8C, P8D, P8E /4000.,3998.,.0001,-2.,
     * 1.5,7.,5.,9.6/
      DATA P8F, P8G, P8H, P9, P10A, P10B, P10C /.6,14.,16.,1.0,.25,.5,
     * .6667/
      END
      REAL FUNCTION F(X, FDERV)
C
C **  AN ARRAY OF TWENTY TEST FUNCTIONS WITH 2 TO 5 DERIVATIVES.
C     THE DERIVATIVE VALUES ARE PLACED IN FDERV - DIMENSIONED AT 5
C     THE FUNCTIONS OF THE SECOND GROUP ARE PARAMETERIZED IN VARIOUS
C     WAYS.  THE PARAMETERS ARE SET IN BLOCK DATA AND AVAILABLE
C     THROUGH THE COMMON BLOCK / FPARS /
C         REQUIRED INPUT INFORMATION IN COMMON BLOCK / FDATA /
C            JFUNK  = INDEX OF THE FUNCTION SELECTED
C            KOUNT  = NUMBER OF F-EVALUATIONS
C            MAXK   = MAXIMUN ALLOWED VALUE OF KOUNT
C         REQUIRED CONTROL INFORMATION IN COMMON BLOCK / KONTROL /
C            FINISH = SWITCH THAT STOPS ADAPT
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL C, CP, CS, C1, DDR, DDT1, DDT2, DDT3, DDT4, DDT5, DDT6,
     * DDT7, DDXEP, DL, DR, DS, DT1, DT2, DT3, DT4, DT5, DT6, DT7,
     * DXEP, D2L, D2R, D2S, D3L, D3R, D3S, EPS, EX, FDERV, FD1, FD2,
     * FK, FL, FLL, FR, FRR, F27, PAR2, PAR3A, PAR3B, PAR4, PAR4A,
     * P10A, P10B, P10C, P5A, P5B, P6A, P6B, P7, P8A, P8B, P8C, P8D,
     * P8E, P8F, P8G, P8H, P9, R, RAND, R0, R4, S, SA, SAX, SAXX, SP,
     * S2, T, T1, T2, T3, T4, T5, T6, T7, X, XC, XE, XEP, XG, XL, XR,
     * X1, X2, X3, X8, Y, YY, YY3, Y2, Y3, BUFFER, SQRTBF, V, SING3,
     * VV, F26
      DIMENSION FDERV(5)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      COMMON /FDATA/ JFUNK, KOUNT, MAXK
      COMMON /FPARS/ PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A,
     * P6B, P7, P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B,
     * P10C
C
      DATA SING3 /.30770707707077/
      DATA BUFFER /1.E-12/
      DATA SQRTBF /1.E-7/
C        DEFINITION OF FUNCTION AT 2700
      F26(VV) = ALOG(2.+VV*VV/(1.+VV))
      F27(V) = SIN(V**2-3.1*EXP(V/(1.+V**2)))*(F26(V)+V)
C
C         FACILITY TO TERMINATE ADAPT TESTING BECAUSE OF EXCESSIVE
C         NUMBER OF FUNCTION EVALUATIONS.
      KOUNT = KOUNT + 1
      IF (KOUNT.NE.MAXK) GO TO 10
      WRITE (6,99999) FMESGE, MAXK
      IF (.NOT.FINISH) FATAL = .TRUE.
      FINISH = .TRUE.
   10 CONTINUE
C
C         PROTECT ARGUMENT X, INITIALIZE DERVS 3,4,5 = 0
      Y = X
      FDERV(3) = 0.0
      FDERV(4) = 0.0
      FDERV(5) = 0.0
C
C         SELECT FUNCTION NUMBER JFUNK
      GO TO (20, 40, 80, 140, 150, 160, 170, 180, 200, 210, 260, 270,
     * 280, 290, 300, 310, 320, 330, 360, 380), JFUNK
C
C        SIMPLE ALBEGRAIC SINGULARITY AT X= 0.0
   20 F = ABS(Y)**.4
      IF (ABS(Y).LE.BUFFER) Y = BUFFER
      FDERV(1) = .4*ABS(Y)**(-.6)*SIGN(1.,Y)
      DO 30 K=2,5
        FK = K
        FDERV(K) = (1.4-FK)*FDERV(K-1)/Y
   30 CONTINUE
      RETURN
C
C         INTERIOR SINGULARITY AT .307707077070770707707077...
C         DIFFERENT EXPONENTS (.481,.63) ON EACH SIDE
   40 T = Y - SING3
      IF (T.LT.0.0) GO TO 60
      F = T**.61
      IF (T.LE.BUFFER) T = BUFFER
      FDERV(1) = .61*F/T
      DO 50 K=2,5
        FDERV(K) = (1.61-FLOAT(K))*FDERV(K-1)/T
   50 CONTINUE
      RETURN
   60 F = (-T)**.481
      FDERV(1) = .481*F/T
      DO 70 K=2,5
        FK = K
        FDERV(K) = (1.481-FK)*FDERV(K-1)/T
   70 CONTINUE
      RETURN
C
C             BROKEN LINE - BREAKS AT 1,2,3.5,4
   80 CONTINUE
      FDERV(2) = 0.0
      IF (Y.LE.1.) GO TO 90
      IF (Y.GE.5.) GO TO 130
      IY = Y + 1.
      GO TO (90, 100, 110, 120, 130), IY
   90 F = 6.*Y + 1.
      FDERV(1) = 6.
      RETURN
  100 F = 7. - (Y-1.)*.5
      FDERV(1) = -.5
      RETURN
  110 F = 6.5 - 1.5*(Y-2.)
      FDERV(1) = -1.5
      RETURN
  120 IF (Y.LT.3.5) GO TO 110
      F = 4.25 - 2.5*(Y-3.5)
      FDERV(1) = -2.5
      RETURN
  130 F = 3. - 3.*(Y-4.)
      FDERV(1) = -3.0
      RETURN
C
C         HUMP AT 2.5 - ONLY 1ST + 2ND DERIVATIVES GIVEN
  140 YY = Y - 2.5
      F = 1./(1.+YY**4)
      YY3 = -4.*YY**3
      FDERV(1) = YY3*F**2
      FDERV(2) = (2.*YY3**2)*F**3 - 12.*YY**2*F**2
      RETURN
C
C           SIMPLE CUBIC
  150 YY = 4.*Y - 2.
      F = YY**3 - YY
      FDERV(1) = 12.*YY**2 - 4.
      FDERV(2) = 96.*YY
      FDERV(3) = 382.
      RETURN
C
C           VARIABLE OSCILLATION - ONLY CORRECT DERIVATIVES ARE 1,2,3
  160 YY = Y**2 - 1.
      F = SIN(YY)
      CS = COS(YY)
      FDERV(1) = 2.*Y*CS
      FDERV(2) = 2.*CS - 4.*(YY+1.)*F
      FDERV(3) = CS*(YY+1.)**3 - 12.*Y*F
      RETURN
C
C           EXPONENTIAL
  170 F = EXP(Y)
      FDERV(1) = F
      FDERV(2) = F
      FDERV(3) = F
      FDERV(4) = F
      FDERV(5) = F
      RETURN
C
C         FIFTH DEGREE POLYNOMIAL PLUS HIGH FREQUENCY SAWTOOTH
C                   ONLY 2 DERIVATIVES GIVEN
  180 F = Y**2*(Y**3-2.) + 6.37
      FD1 = 5.*Y**4 - 4.*Y
      FD2 = 20.*Y**3 - 4.
C         EPS IS FRACTIONAL PART OF 1000*Y
      L = 1000.*Y
      EPS = Y - .001*FLOAT(L)
      L = L - 2*(L/2)
      IF (L.EQ.1) GO TO 190
C             INCREASING PART OF THE SAWTOOTH
      F = F + 500.*EPS**2
      FDERV(1) = FD1 + 1000.*EPS
      FDERV(2) = FD2 + 1000.
      RETURN
C             DECREASING PART OF THE SAWTOOTH
  190 CONTINUE
      F = F + 500.*(1.E-6-EPS**2)
      FDERV(1) = FD1 + 1. - 1000.*EPS
      FDERV(2) = FD2 - 1000.
      RETURN
C
C           SEVENTH DEGREE POLYNOMIAL
  200 F = Y**7 - 2.*Y**6 + 3.1*Y**4 - 6.2*Y
      FDERV(1) = 7.*Y**6 - 12.*Y**5 + 12.4*Y**3 - 6.2
      FDERV(2) = 42.*Y**5 - 60.*Y**4 + 37.2*Y**2
      FDERV(3) = 210.*Y**4 - 240.*Y**3 + 74.4*Y
      FDERV(4) = 840.*Y**3 - 720.*Y**2 + 74.4
      FDERV(5) = 2520.*Y**2 - 1440.*Y
      RETURN
C
C         WILD VARIATION BETWEEN .47 AND .48 - ONLY 3 DERIVATIVES HERE
C         THERE ARE JUMPS AND INFINITIES IN THE SLOPE.
  210 FDERV(2) = 0.
      IF (ABS(Y-.35).GT..13) GO TO 230
      IF (ABS(Y-.45).GT..03) GO TO 240
      IF (ABS(Y-.475).GT..005) GO TO 250
C                SMALL PART - PROB NOT EXACTLY CONTINOUS
      Y2 = -(Y**2-.95*Y+.2256)
      Y3 = ABS(Y2)**.3
      F = 547.5*Y - 259.39 + 26.7*Y3
      IF (Y2.LE.BUFFER) GO TO 220
      FDERV(1) = 547.5 + 8.01*Y3/Y2*(2.*Y-.95)
      FDERV(2) = (-11.214*(2.*Y-.95)**2/Y2+16.02)*Y3/Y2
      FDERV(3) = ((19.0638*(2.*Y-.95)/Y2-44.856)*(2.*Y-.95)-11.214)*Y3/
     * Y2**2
      RETURN
C             SET DERIVATIVES = 0 AT BREAKS
  220 FDERV(1) = 0.0
      FDERV(2) = 0.0
      FDERV(3) = 0.0
      RETURN
C                 MAIN LINEAR PIECE
  230 F = 6.*Y + .55
      FDERV(1) = 6.
      RETURN
C                 NEXT LINEAR PART
  240 F = 25.2*Y - 3.674
      FDERV(1) = 25.2
      RETURN
C                 SMALL LINEAR PART
  250 F = -179.05*Y + 82.111
      FDERV(1) = -179.05
      RETURN
C
C          LOTS OF OSCILLATIONS - FIRST EXAMPLE - ONLY 2 DERIVATIVES
  260 CONTINUE
      X1 = .5*Y**2 + Y
      X2 = .1*Y**3 - 1.
      S = SIN(X1)
      C1 = COS(X1)
      C = COS(X2)
      S2 = SIN(X2)
      R0 = 1./(1.+Y**6)
      R4 = 1./(1.+(Y-4.)**6)
      SP = C1*(Y+1.)
      CP = -S2*(.3*Y**2)
      F = (S+C)*(R0+R4)
      FDERV(1) = -6.*(S+C)*(R0**2*Y**5+R4**2*(Y-4.)**5) +
     * (R0+R4)*(SP+CP)
      FDERV(2) = (S+C)*(-30.*(Y**4*R0**2+(Y-4.)**4*R4**2)+72.*(Y**10*R0*
     * *3+(Y-4.)**10*R4**3)) - 12.*(R0**2*Y**5+(Y-4.)**5*R4**2)*(SP+CP)
     * + (R0+R4)*(C1-(Y+1.)**2*S-.09*Y**4*C-.6*Y*S2)
      RETURN
C
C          LOTS OF OSCILLATIONS - SECOND EXAMPLE - ONLY 3 DERIVATIVES
  270 CONTINUE
      EX = EXP(Y)
      R = 1./(Y+PAR2)
      S = SIN(R)
      C = COS(R)
      DS = -C*R**2
      D2S = -S*R**4 + 2.*C*R**3
      D3S = C*R**6 + 6.*S*R**5 - 6.*C*R**4
      F = EX*S
      FDERV(1) = EX*(S+DS)
      FDERV(2) = EX*(S+2.*DS+D2S)
      FDERV(3) = EX*(S+3.*DS+3.*D2S+D3S)
      RETURN
C
C         END POINT SINGURITIES OF PAR3A,PAR3B - ONLY 3 DERIVATIVES
  280 CONTINUE
      IF (Y.LE.BUFFER) Y = BUFFER
      IF (Y.GE.1.-BUFFER) Y = 1. - BUFFER
      FL = Y**PAR3A
      FR = (1.-Y)**PAR3B
      DL = PAR3A*FL/Y
      DR = -PAR3B*FR/(1.-Y)
      D2L = (PAR3A-1.)*DL/Y
      D2R = -(PAR3B-1.)*DR/(1.-Y)
      D3L = (PAR3A-2.)*D2L/Y
      D3R = -(PAR3B-2.)*D2R/(1.-Y)
      F = FL*FR
      FDERV(1) = FL*DR + DL*FR
      FDERV(2) = FL*D2R + 2.*DL*DR + D2L*FR
      FDERV(3) = FL*D3R + 3.*DL*D2R + 3.*D2L*DR + D3L*FR
      RETURN
C
C         RANDOM NUMBERS ADDED TO SIN(X) - ONLY TWO DERIVATIVES
  290 CONTINUE
      S = SIN(Y)
      C = COS(Y)
C                   QUICKIE PSEUDO-RANDOM NUMBER GENERATOR USED HERE
      RAND = AMOD(4829.*S+C,7.23)/7.23
      F = S + PAR4*(RAND-.5)
      RAND = AMOD(2993.*F+S,3.19)/3.19
      FDERV(1) = C + PAR4*(RAND-.5)
      RAND = AMOD(4901.*RAND+C,8.13)/8.13
      FDERV(2) = -S + PAR4*(RAND-.5)
      RAND = AMOD(6987.*RAND+F*C,4.27)/4.27
      FDERV(3) = -C + PAR4*(RAND-.5)
      RETURN
C
C         LARGE FUNCTION
  300 CONTINUE
      EX = EXP(Y)*P5A
      X8 = Y**8*P5B
      F = EX + X8
      FDERV(1) = EX + 8.*X8/Y
      FDERV(2) = EX + 56.*X8/Y**2
      FDERV(3) = EX + 336.*X8/Y**3
      FDERV(4) = EX + 1680.*Y**4*P5B
      FDERV(5) = EX + 6720.*Y**3*P5B
      RETURN
C
C          SHIFTED FUNCTION   - ONLY 2 DERIVATIVES
  310 CONTINUE
      EX = EXP(Y-P6A)
      X3 = (Y-P6B)**3
      S = SIN(Y)
      C = COS(Y)
      R = 1./(1.+.5*S)
      F = EX + X3*R
      DR = -.5*C*R**2
      D2R = .5*S*R**2 + .5*C**2*R**3
      FDERV(1) = EX + 3.*(Y-P6B)**2*R + X3*DR
      FDERV(2) = EX + 6.*(Y-P6B)*R + 6.*(Y-P6B)**2*DR + X3*D2R
      RETURN
C
C         NUMERICAL DIFFERENTIATON OF ARITHMETIC STATEMENT FUNCTION
C             ONLY  FIRST, 2ND AND THIRD DERIVATIVES
  320 CONTINUE
      F = F27(Y)
      FR = F27(Y+P7)
      FL = F27(Y-P7)
      FRR = F27(Y+2.*P7)
      FLL = F27(Y-2.*P7)
      FDERV(1) = (FR-FL)/P7*.5
      FDERV(2) = (FR-2.*F+FL)/P7**2
      FDERV(3) = (-FLL+2.*FL-2.*FR+FRR)/P7**3*.5
      RETURN
C
C         RATHER COMPLEX FUNCTION - SUM OF 7 DIFFERENT THINGS
C             ONLY TWO DERIVATIVES
C             INFINITE SLOPE AT P8B AND CUSP AT P8E POSSIBLE
  330 CONTINUE
C                   TERM 1
      SA = SQRT(ABS(Y-P8A))
      T1 = EXP(-SA)
      IF (SA.LE.SQRTBF) SA = SQRTBF
      DT1 = -T1*SIGN(1.,Y-P8A)*.5/SA
      DDT1 = T1*(1.+1./SA)*.25/(SA*SA)
C                   TERM 2
      SA = 4.*EXP(-ABS(Y-P8B))
      SAX = Y*SA
      SAXX = Y*SAX
      T2 = 2. - SAXX
      DT2 = -2.*SAX + SAXX*SIGN(1.,Y-P8B)
      DDT2 = -2.*SA + 4.*SAX*SIGN(1.,Y-P8B) - SAXX
C                   TERM 3
      XC = Y - P8C
      T3 = 1./(1.+150.*XC**6)
      DT3 = -900.*XC**5*T3**2
      DDT3 = 1620000.*XC**10*T3**3 - 4500.*XC**4*T3**2
C                   TERM 4
      T4 = 6./(1.+P8D*XC**4)
      DT4 = -2./3.*P8D*XC**3*T4**2
      DDT4 = 8./9.*P8D**2*XC**6*T4**3 - 2.*P8D*XC**2*T4**2
C
C                   TERM 5
      XE = Y - P8E
      XEP = ABS(XE)**P8F
      R = 10./(1.+4.*XE**4)
      T5 = R*XEP
      DR = -1.6*XE**3*R**2
      DDR = 5.12*XE**6*R**3 - 4.8*XE**2*R**2
      DXEP = 0.0
      DDXEP = 0.0
      IF (ABS(XE).LE.1.E-12) GO TO 340
      DXEP = P8F*XEP/XE
      DDXEP = (P8F-1.)*DXEP/XE
  340 DT5 = DR*XEP + R*DXEP
      DDT5 = DDR*XEP + 2.*DR*DXEP + R*DDXEP
C
C                   TERM 6
      XG = ABS(Y-P8G)**2.5
      T6 = -2.*EXP(-XG)
      DT6 = -2.5*T6*SIGN(1.,Y-P8G)*XG**.6
      DDT6 = T6*6.25*(XG**1.2-XG**.2*.6)
C
C                   TERM 7
      S = SIN(4.*Y)
      C = COS(4.*Y)
      R = 1./(.5+ABS(Y-P8H)*5.)
      EX = EXP(Y+1.-P8H)
C             ACTUALLY USING EXP(AMAX1(X-15,0)) - SPLIT CALCULATION
      T7 = S*R
      DT7 = 4.*C*R - S*R**2*SIGN(1.,Y-P8H)*5.
      DDT7 = -16.*S*R - 40.*C*R**2*SIGN(1.,Y-P8H) + 50.*S*R**3
      IF (Y.LT.P8H-1.) GO TO 350
C                        INCLUDE EX TERM
      DDT7 = EX*(DDT7+2.*DT7+T7)
      DT7 = EX*(DT7+T7)
      T7 = EX*T7
  350 CONTINUE
C
      F = T1 + T2 + T3 + T4 + T5 + T6 + T7
      FDERV(1) = DT1 + DT2 + DT3 + DT4 + DT5 + DT6 + DT7
      FDERV(2) = DDT1 + DDT2 + DDT3 + DDT4 + DDT5 + DDT6 + DDT7
      RETURN
C
C         JUMP DISCONTINUITY (EXCEPT FOR P9 = 0.) - ONLY 2 DERIVATIVES
  360 CONTINUE
      IF (Y.GE.P9) GO TO 370
C                     EXPONENTIAL
      F = EXP(Y)
      FDERV(1) = F
      FDERV(2) = F
      FDERV(3) = F
      RETURN
  370 CONTINUE
C                     RATIONAL
      F = 1./(1.+Y**2)
      FDERV(1) = -2.*Y*F**2
      FDERV(2) = 2.*F**2*(4.*Y*F-1.)
      RETURN
C
C                    INFINITE SINGULARITY
  380 CONTINUE
      IF (Y.GT.P10A+BUFFER) GO TO 390
      IF (Y.GT.P10A-BUFFER) GO TO 400
C              LEFT SIDE OF SINGULARITY
      XL = P10A - Y
      F = XL**P10B
      FDERV(1) = -P10B*F/XL
      FDERV(2) = -(P10B-1.)*FDERV(1)/XL
      FDERV(3) = -(P10B-2.)*FDERV(2)/XL
      RETURN
  390 CONTINUE
C              RIGHT SIDE OF SINGULARITY
      XR = Y - P10A
      F = XR**P10C
      FDERV(1) = P10C*F/XR
      FDERV(2) = (P10C-1.)*FDERV(1)/XR
      FDERV(3) = (P10C-2.)*FDERV(2)/XR
      RETURN
  400 CONTINUE
C              AT THE SINGULARITY
      F = 100./(BUFFER**AMAX1(P10B,P10C))
      FDERV(1) = 0.0
      FDERV(2) = F**3
      RETURN
99999 FORMAT (//2X, 6A4, 14HEXCEEDED LIMIT, I5, 14H ON F(X) EVALS)
      END
      REAL FUNCTION QPOLY(T)
C
C **  THIS IS FUNCTION EXACTLY LIKE PPOLY EXCEPT IT HAS JUST 1 ARGUMENT
C         SO IT CAN BE USED WITH ERRINT BY SUMRY2.
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL T, X, ZCOEF, ZKNOTS
      DIMENSION ZKNOTS(100), ZCOEF(100,7)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      COMMON /PARAMS/ ZKNOTS, ZCOEF
      COMMON /SAVEIT/ IKNOT
   10 IF (T.LT.ZKNOTS(IKNOT+1)) GO TO 20
      IKNOT = IKNOT + 1
      IF (IKNOT.LT.KNOTS) GO TO 10
      IKNOT = KNOTS - 1
      GO TO 30
   20 IF (T.GE.ZKNOTS(IKNOT)) GO TO 30
      IKNOT = IKNOT - 1
      IF (IKNOT.LT.1) GO TO 20
      IKNOT = 1
   30 X = T - ZKNOTS(IKNOT)
      QPOLY = ZCOEF(IKNOT,NPAR)
      DO 40 K=1,DEGREE
        KK = NPAR - K
        QPOLY = ZCOEF(IKNOT,KK) + X*QPOLY
   40 CONTINUE
      RETURN
      END
      SUBROUTINE SUMRY2(XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C        =============================================================
C
C **  THIS PROGRAM SUMMARIZES THE PERFORMANCE OF ADAPT
C     IT COMPUTES EXECUTION TIME, INDEPENDENTLY CHECKS THE ACCURACY ,
C         PLOTS F(X) PLUS THE ERROR FOR 100 POINTS.
C     IT USES THE SYSTEM DEPENDENT ROUTINES SECOND AND GRAPH
C         AND THE SUBPROGRAM ERRINT FROM ADAPT.
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      REAL ABSC(4), DX, DXG, ECHECK, ERRINT, ERRP, P, PPOLY, QPOLY,
     * WGTS(4), XL, XR, ZCOEF, ZKNOTS, F
      REAL XVALU(100), GRAF1(100), GRAF2(100), TIME, TSTART, TSTOP,
     * ZKNOTS(100), ZCOEF(100,7), DUMMY(6)
      EXTERNAL F, QPOLY
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      COMMON /PARAMS/ ZKNOTS, ZCOEF
      COMMON /FDATA/ JFUNK, KOUNT, MAXK
      COMMON /INFO/ ECHECK, KFUNK, TIME, TSTART
      DATA ABSC(1) /-.86113631159405/
      DATA ABSC(2) /-.33998104358486/
      DATA ABSC(3) /.33998104358486/
      DATA ABSC(4) /.86113631159405/
      DATA WGTS(1) /.34785484513745/
      DATA WGTS(2) /.65214515486255/
      DATA WGTS(3) /.65214515486255/
      DATA WGTS(4) /.34785484513745/
C
      IF (LEVEL.LE.1) RETURN
C
C                  LEVEL 2 OUTPUT
C
C     ***** CALL CLOCK FOR VALUE OF  TSTOP  IN SECONDS *****
C
      TSTOP = 0.0
      TIME = TSTOP - TSTART
      KFUNK = KOUNT
C
C            VERIFY ERROR BY INDEPENDENT CHECK
      ECHECK = 0.
      IF (KNOTS.GT.100 .OR. NPAR.GT.7) GO TO 40
C             SKIP THIS CHECK
C
C      PUT PARAMETERS IN THE COMMON BLOCK PARAMS FOR QPOLY TO USE
      KNOTS1 = KNOTS - 1
      DO 20 J=1,KNOTS1
        ZKNOTS(J) = XKNOTS(J)
        DO 10 K=1,NPAR
          ZCOEF(J,K) = COEFS(J,K)
   10   CONTINUE
   20 CONTINUE
C
C         BREAK (A,B) INTO HALF AGAIN AS MANY PARTS AS KNOTS
      IPART = (3*KNOTS)/2 + 1
      DX = (B-A)/FLOAT(IPART)
      XL = A
      XR = A + DX
      ECHECK = 0.0
C             COMPUTE ERROR ESTIMATE OVER PART
      P = ABS(NORM)
      DO 30 K=1,IPART
        ERRP = ERRINT(F,QPOLY,XL,XR,ABSC,WGTS)
        XL = XR
        XR = XR + DX
C             UPDATE ERROR
        IF (NORM.EQ.3.) ECHECK = AMAX1(ERRP,ECHECK)
        IF (NORM.NE.3.) ECHECK = (ECHECK**P+ERRP)**(1./P)
   30 CONTINUE
   40 WRITE (6,99999) KFUNK, ACCUR, ERROR, ECHECK
      WRITE (6,99998) TIME
C
C            GRAPHICAL OUTPUT
      RETURN
C         ***** IF GRAPH ROUTINE IS AVAILABLE REMOVE THE ABOVE RETURN
C               AND CALL GRAPH ROUTINE BELOW
      WRITE (6,99997)
      DXG = (B-A)/99.
      DO 50 K=1,100
        XVALU(K) = A + FLOAT(K-1)*DXG
C                  DUMMY IS A DUMMY ARRAY BELOW
        GRAF1(K) = F(XVALU(K),DUMMY)
        GRAF2(K) = GRAF1(K) - PPOLY(XVALU(K),XKNOTS,COEFS,KDIMEN,NDIMEN)
   50 CONTINUE
C     ***** INSERT GRAPH CALL HERE FOR 100 POINTS OF (XVALU,GRAF1,GRAF2)
      RETURN
99999 FORMAT (/10X, 10HADAPT USED, I5, 28H FUNCTION VALUES FOR ERRORS
     * /10X, 11HSPECIFIED =, E15.5, 21H ESTIMATED BY ADAPT =, E15.5,
     * 21H INDEPENDENT CHECK = , E15.5)
99998 FORMAT (20X, 11HTIME USED =, F8.5, 9H SECONDS )
99997 FORMAT (41H1       PLOT OF F(X) AND THE ERROR CURVE )
      END
      SUBROUTINE ADAPT(F, XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C     THIS ALGORITHM COMPUTES A PIECEWISE POLYNOMIAL APPROXIMATION
C  OF SPECIFIED SMOOTHNESS, ACCURACY AND DEGREE.  THE INPUT TO THE
C  COMPUTATION IS
C
C   F      - FUNCTION BEING APPROXIMATED. IT MUST PROVIDE VALUES OF
C            DERIVATIVES UP TO THE ORDER OF SMOOTHNESS SPECIFIED FOR
C            THE APPROXIMATION.  THE CALLING SEQUENCE IS F(X,FDERV) AND
C            FDERV CONTAINS THE DERIVATIVES( SEE CONSTRAINT BELOW)
C   A,B    - THE ENDPOINTS OF THE INTERVAL OF APPROXIMATION
C   ACCUR  - THE ACCURACY REQUIRED FOR THE APPROXIMATION
C   SMOOTH - THE SMOOTHNESS REQUIRED FOR THE APPROXIMATION
C              = 0  MEANS CONTINUOUS
C              = 1  MEANS CONTINUOUS SLOPE
C              = 2  MEANS CONTINUOUS SECOND DERIVATIVE, ETC.
C   DEGREE - THE DEGREE OF THE POLYNOMIAL PIECES.
C            MUST HAVE DEGREE GT 2*SMOOTH
C
C     *****  ***** SECONDARY INPUT - ITEMS WITH DEFAULT VALUES POSSIBLE
C   CHARF  - CHARACTERISTIC LENGTH OF THE FUNCTION F(X). PIECES ARE NOT
C            LONGER THAN THIS LENGTH.
C                   DEFAULT=(B-A) IF DEGREE GT 1, ELSE (B-A)/3
C   NORM   - NORM TO MEASURE THE APPROXIMATION ERROR
C              = 1  L1 APPROXIMATION (LEAST DEVIATIONS)
C              = 2  L2 APPROXIMATION (LEAST SQUARES)
C              = 3  TCHEBYCHEFF (MINIMAX) APPROXIMATION
C              =-P  (NEGATIVE VALUE) GENERAL LP APPROXIMATION
C                   DEFAULT= 2
C   NBREAK - NUMBER OF SPECIAL BREAK POINTS IN THE APPROXIMATION.
C            ASSOCIATED INPUT VARIABLES ARE
C              XBREAK(J)  - LOCATION OF BREAK POINTS
C              DBREAK(J)  - DERIVATIVE BROKEN AT XBREAK
C              BLEFT (J)  - VALUE FROM LEFT FOR DBREAK DERIVATIVE
C              BRIGHT(J)  -   -    -   RIGHT -    -        -
C                   DEFAULT = 0
C   LEVEL  - LEVEL OF OUTPUT FROM ADAPT
C              = -1 NO OUTPUT AT ALL EXCEPT FOR FATAL INPUT ERRORS
C              =  0  ERROR CONDITIONS NOTED, FINAL SUMMARY
C              =  1  PRINT THE APPROXIMATION OUT
C              =  2  DETAILS OF THE COMPUTATION
C              =  3 DEBUG OUTPUT,  = 4 LOTS OF DEBUG OUTPUT
C                    DEFAULT = 0
C   EDIST  - SWITCH TO CHANGE FROM PROPORTIONAL ERROR DISTRIBUTION
C            TO FIXED DISTRIBUTION. THIS IS PRIMARILY OF USE IN
C            APPROXIMATION OF FUNCTIONS WITH SINGULARITIES. ONE SHOULD
C            USE NORM = 1. OR SO IN SUCH CASES
C              = 0  PROPORTIONAL DISTRIBUTION
C              = 1  APPROXIMATE FIXED ERROR DISTRIBUTION
C                   ATTEMPTS TO ACHIEVE SPECIFIED ACCURACY VALUE ACCUR
C              = 2  TRUE FIXED ERROR DISTRIBUTION
C
C     THE OUTPUT OF THE COMPUTATION CONSISTS OF 3 PARTS, EACH RETURNED
C     TO THE USER IN A DIFFERENT WAY. THEY ARE
C
C   XKNOTS,COEFS - ARRAYS DEFINING THE PIECEWISE POLYNOMIAL RESULT.
C   COEFS    XKNOTS(K)  = KNOTS OF THE APPROXIMATION ( K = 1 TO KNOTS)
C                         THE LAST ONE IS RIGHT END POINT OF INTERVAL
C            COEFS(K,N) = COEFFICIENT OF (X - XKNOT(K))**(N-1) IN THE
C                         INTERVAL XKNOT(K) TO XKNOT(K+1)
C                           K = 1 TO KNOTS-1  AND  N = 1 TO DEGREE+1
C            THESE ARRAYS ARE PASSED AS ARGUMENTS SO AS TO USE VARIABLE
C            DIMENSIONS.
C               KDIMEN - DIMENSION USED FOR  XKNOTS IN CALLING PROGRAM
C               NDIMEN - COEFS IS DECLARED  COEFS(KDIMEN,NDIMEN) IN THE
C                        CALLING PROGRAM.
C               ***** NOTE ***** SEVERAL SMALL ARRAYS HERE HAVE FIXED
C                        DIMENSIONS THAT LIMIT DEGREE AND THUS NDIMEN
C                        SHOULD NOT EXCEED THIS LIMIT (CURRENTLY = 12)
C
C   PPOLY  - THE PIECEWISE POLYNOMIAL APPROXIMATING FUNCTION.
C            THIS FUNCTION SUBPROGRAM IS AVAILABLE TO THE USER AT THE
C            COMPLETION OF ADAPT.
C
C   RESULZ - A LABELED COMMON BLOCK WITH  ERROR,KNOTS  IN IT
C              KNOTS - NUMBER OF KNOTS OF PPOLY
C              ERROR - ACCURACY OF PPOLY AS ESTIMATED BY ADAPT
C
C  ********** DIMENSION CONSTRAINTS **********
C     MAXKNT - MAX NUMBER OF KNOTS TAKEN FROM USER VIA KDIMEN
C              ARRAYS WITH THIS DIMENSION
C                   COEFS   XKNOTS
C     MAXPAR - MAX NUMBER OF PARAMETERS PER INTERVAL ( = 12 CURRENTLY )
C                   USER PROVIDED NDIMEN MUST HAVE NDIMEN LE MAXPAR
C                   MUST HAVE  DEGREE + 1 LE MAXPAR
C              ARRAYS WITH THIS DIMENSION (OR RELATED VALUES )
C                D       DDTEMP  FDERVL  FDERVR  FDUMB   FACTOR
C                FINTRP  FLEFT   FRIGHT  POWERS  XTEMP   XINTRP  XDD
C              ***** NOTE ***** MAXPAR ALSO AFFECTS ARGUMENT FDERV
C              OF FUNCTION F.  FDERVL, FDERVR ARE ALSO INVOLVED.
C              SHOULD DECLARE FDERV OF SIZE 6 IN F TO BE SAFE.
C     MAXAUX - MAXIMUM NUMBBER OF AUXILIARY INPUT( = 20 NOW ). ARRAYS
C                 XBREAK  DBREAK  BLEFT   BRIGHT
C     MAXSTK - MAX SIZE OF ACTIVE INTERVAL STACK
C              MIN INTERVAL LENGTH IS 2**(-MAXSTK)*(B-A). ARRAYS
C                 XLEFT   XRIGHT
C
C  **********  PORTABILITY CONSIDERATIONS  **********
C
C     THIS PROGRAM IS IN ANSI STANDARD FORTRAN.  IN ADDITION, IT MEETS
C     ALL THE REQUIREMENTS OF THE BELL LABS PORTABLE FORTRAN -PFORT-
C     EXCEPT ONE.  HOLLERITH CHARACTERS ARE PACKED 4/WORD RATHER THAN
C     1/WORD AS SPECIFIED BY PFORT.
C     NEVERTHELESS, THIS PROGRAM IS AFFECTED IN SEVERAL WAYS BY A
C     CHANGE IN MACHINE WORD LENGTH AND CHANGING TO DOUBLE PRECISION.
C     ***** THIS VERSION IS FOR THE MACHINE WITH THE LONGEST SINGLE
C           PRECISION WORD (CDC).  THE LENGTH OF SOME CONSTANTS IN
C           THE SUBPROGRAM COMPUT EXCEEDS THE CAPACITY OF SOME FORTRAN
C           COMPILERS AND CAN PREVENT COMPILATION.
C  INPUT-OUTPUT -- IS OF THREE TYPES.
C         FATAL ERROR MESSAGES - OCCUR IN SETUP,TAKE,PUT AND TERMIN
C                                THEY CANNOT BE SWITCHED OFF
C         USER AND DEBUGGING OUTPUT - CAN BE SWITCHED OFF
C     THESE INVOLVE MANY FORMATS LIKE  E15.8, F12.8, E22.13, ETC.
C     SOME FORTRAN COMPILERS REQUIRE D-FORMAT FOR DOUBLE PRECISION.
C     SOME DO NOT HANDLE E22.13 ON MACHINES OF SHORTER WORD LENGTH.
C
C     SUMARY PRINTS COEFFICIENTS 5 PER LINE, 6 PER LINE IS BETTER
C     FOR SHORTER WORD LENGTHS.  DOUBLE PRECISION ON MANY MACHINES
C     LIMITS ONE TO 4 PER LINE.
C
C  CONSTANTS -- THE GAUSS WEIGHTS AND ABSCISSA IN COMPUTE ARE GIVEN
C               TO 15 DIGITS IN THE COMMENTS
C
C  DOUBLE PRECISION CONVERSION -- REQUIRES FOUR STEPS
C     1. ALL REAL VARIABLES ARE DECLARED DOUBLE PRECISION. THIS IS
C        DONE BY CHANGING REAL TO DOUBLE PRECISION AS ALL REALS ARE
C        EXPLICITLY DECLARED AND ROOM IS LEFT FOR THIS CHANGE.  REAL
C        VARIABLES APPEAR BEFORE INTEGERS IN ALL COMMON BLOCKS.
C
C     2. ADD D TO CONSTANTS( E.G. 1.0 = 1.0D0 ). ADJUST LENGTH OF
C        GAUSS WEIGHTS IN COMPUTE.
C
C     3. CHANGE ABS,AMAX1,FLOAT  AT MANY PLACES
C
C     4. ADJUST THE INTERVAL STACK SIZE = DIMENSIONS OF XLEFT, XRIGHT
C        AND VALUE OF MAXSTK.  ADJUST THE VALUE BUFFER IN PUT TO BE
C        ABOUT .1E-K FOR A MACHINE WITH K+2 DECIMAL DIGITS.
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      EXTERNAL F
      REAL F
C
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
C                      KNTDIM - KDIMEN, NAME CHANGED TO PUT IN COMMON
C                      NPARDM - NDIMEN, NAME CHANGED TO PUT IN COMMON
      COMMON /RESULZ/ ERROR, KNOTS
C                      KNOTS = FINAL NO. OF KNOTS, INCLUDES B AS ONE.
C                      ERROR = ESTIMATE OF ERROR ACTUALLY ACHIEVED.
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
C             KONTRL CONTAINS GENERALLY USEFUL VARIABLES
C                      FATAL  - SWITCH FOR DETECTION OF FATAL ERROR
C                               INCLUDING EXCESSIVE INTERVAL SUBDIVISION
C                               WHICH DOES NOT TERMINATE THE COMPUTATION
C                      FINISH - SWITCH TO TERMINATE ALGORITHM
C                      MAXS   - SEE COMMENTS EARLIER
C                      NSTACK - COUNTER FOR INTERVAL STACK, CONSISTS OF
C                               (XLEFT(J),XRIGHT(J))  J = 1 TO NSTACK
C                      ERRORI - ERROR ESTIMATE FOR TOP INTERVAL
C                      DSCTOL - TOLERANCE TO CHECK DISCARDING INTERVALS
C                      DISCRD - SWITCH TO SIGNAL DISCARD OF TOP INTERVAL
C                      FACTOR - ARRAY OF FACTORIALS
C                      NPAR   - NUMBER OF PAREMETERS = DEGREE + 1
C                      FMESGE - STRING =  **** FATAL ERROR *****
C                      INTERP - NUMBER OF INTERIOR INTERPOLATION POINTS
C                               IN THE NORMAL INTERVAL
C                      IBREAK - COUNTER ON BREAK POINTS
C                      BREAK  - SWITCH FOR BREAK POINT IN TOP INTERVAL
C                               0     = NO BREAK PRESENT
C                               LEFT  = BREAK AT XLEFT(NSTACK)
C                               RIGHT = BREAK AT XRIGHT(NSTACK)
C                               BOTH  = BREAK AT BOTH ENDS
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C             COMDIF CONTAINS VARIABLES USED ONLY BY COMPUT AND FRIENDS.
C                      NINTRP - NUMBER OF INTERIOR INTERPOLATION POINTS
C                               FOR THE CURRENT INTERVAL
C                      XINTRP - INTERIOR INTERPOLATION POINTS
C                      FINTRP - F VALUES AT XINTRP POINTS
C                      LEFTX  - MULTIPLICITY OF INTERPOLATION AT XLEFT
C                               = NO. OF DERIVATIVES MATCHED AT XLEFT
C                      FLEFT  - VALUES OF F AND ITS DERIVATIVES AT XLEFT
C                      RIGHTX - MULTIPLICITY OF INTERPOLATION AT XRIGHT
C                      FRIGHT - VALUES OF F AND DERIVATIVES AT XRIGHT
C                      DDTEMP - THE ARRAY OF DIVIDED DIFFERENCES
C                      XDD    - THE X VALUES FOR DDTEMP WITH PROPER
C                               MULTIPLICITIES OF XLEFT AND XRIGHT
      COMMON /SAVEIT/ IKNOT
C
C------------------------ MAIN CONTROL PROGRAM -------------------------
C
C   ***** NOTE - ARGUMENTS BELOW ARE FOR READABILITY ONLY    *****
C   *****        EXCEPT FOR F AND XKNOTS,COEFS,KDIMEN,NDIMEN *****
C
C                   CHECK INPUT, INITIAL COMPUTATIONS, PRINT PROBLEM
C
      CALL SETUP(XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C         CHECK FOR FATAL ERROR IN PROBLEM SPECIFICATION
      IF (FATAL) RETURN
C
C                   LOOP OVER PROCESSING OF INTERVALS
   10 CALL TAKE(INTERV)
C
      CALL COMPUT(F, APPROX, INTERV)
C
C                   CHECK FOR DISCARDING INTERVALS
      CALL CHECK(FUNCT, CHARCT)
C
C            PUT NEW INTERVALS ON STACK OR DISCARD, UPDATE STATUS
      CALL PUT(INTERV, XKNOTS, COEFS, KDIMEN, NDIMEN)
C
      CALL TERMIN(TEST, AND, PRINT, XKNOTS, KDIMEN)
C
      IF (.NOT.FINISH) GO TO 10
C
      CALL SUMARY(XKNOTS, COEFS, KDIMEN, NDIMEN)
      RETURN
      END
      SUBROUTINE CHECK(FUNCT, CHAR)
C
C        ===============================================================
C
C **  THIS PROGRAM CHECKS FOR DISCARDING INTERVAL, APPLIES VARIOUS
C     TESTS ABOUT DISCARDING INVOLVING EDIST AND CHARF.
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL DTEST, DX
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C
      DISCRD = .FALSE.
      DX = XRIGHT(NSTACK) - XLEFT(NSTACK)
C
C         COMPUTE DTEST FOR IMPLEMENTING VARIOUS TYPES OF ADAPTIVE APPRX
      IF (EDIST.EQ.0) DTEST = DX*DSCTOL
C         FOR THE APPROXIMATE FIXED ERROR DISTRIBUTION TYPE WE ESTIMATE
C         THE FINAL NUMBER OF KNOTS BY( LIMITING IT A LITTLE )
C             (NSTACK+KNOTS+2)((B-A)/(XRIGHT-A))
      IF (EDIST.EQ.1) DTEST = DSCTOL/(FLOAT(NSTACK+KNOTS+2)*AMIN1((B-A)/
     * (XRIGHT(NSTACK)-A),5.))
      IF (EDIST.EQ.2 .OR. NORM.EQ.3.) DTEST = DSCTOL
C
C             CHECK FOR DISCARD OF INTERVAL
      IF (ERRORI.LE.DTEST) DISCRD = .TRUE.
C
C         CHECK CHARACTERISTIC LENGTH OF FUNCTION
      IF (DX.GE.CHARF) DISCRD = .FALSE.
      RETURN
      END
      SUBROUTINE COMPUT(F, APPROX, INTERV)
C
C        ===============================================================
C
C **  THIS PROGRAM COMPUTES THE PIECEWISE POLYNOMIAL APPROXIMATION ON
C     THE CURRENT INTERVAL. IT ALSO ESTIMATES THE ERROR
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL ABSC, DX, ERRINT, F, FDERVL, FDERVR, FDUMB, POLYDD, WGTS
      DIMENSION ABSC(4), WGTS(4), FDERVL(5), FDERVR(5), FDUMB(5)
      EXTERNAL F, POLYDD
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      EQUIVALENCE (FLEFT(2),FDERVL(1)), (FRIGHT(2),FDERVR(1))
C     FIFTEEN DIGIT VALUES FOR THESE GAUSS INTEGRATION CONSTANTS ARE
C  .861136311594053 .339981043584856 .347854845137454 .652145154862546
      DATA ABSC(1) /-.86113631159405/
      DATA ABSC(2) /-.33998104358486/
      DATA ABSC(3) /.33998104358486/
      DATA ABSC(4) /.86113631159405/
      DATA WGTS(1) /.34785484513745/
      DATA WGTS(2) /.65214515486255/
      DATA WGTS(3) /.65214515486255/
      DATA WGTS(4) /.34785484513745/
C
C             COMPUTE INTERPOLATION INFORMATION
      NINTRP = DEGREE - 2*SMOOTH - 1
C
C          INCREASE NUMBER OF INTERPOLATION POINTS IF BREAK POINTS ARE
C          SPECIFIED WITH FEWER DERIVATIVES THAN SMOOTH
      IF (BREAK.EQ.LEFT .OR. BREAK.EQ.RIGHT) NINTRP = NINTRP + SMOOTH -
     * DBREAK(IBREAK)
      IF (BREAK.EQ.BOTH) NINTRP = NINTRP + 2*SMOOTH - DBREAK(IBREAK) -
     * DBREAK(IBREAK+1)
      IF (NINTRP.EQ.0) GO TO 20
C             GENERATE EQUAL SPACED INTERPOLATION POINTS
      DX = (XRIGHT(NSTACK)-XLEFT(NSTACK))/FLOAT(NINTRP+1)
      DO 10 J=1,NINTRP
        XINTRP(J) = XLEFT(NSTACK) + FLOAT(J)*DX
   10 CONTINUE
C
C             GET LEFT AND RIGHT F-VALUES, PUT F-VALUE IN FIRST ELEMENT
C             OF ARRAYS FLEFT AND FRIGHT.  GET DERIVATIVES BACK AS
C             OTHER ELEMENTS VIA THE SUBARRAYS FDERVL AND FDERVR.
   20 FLEFT(1) = F(XLEFT(NSTACK),FDERVL)
      FRIGHT(1) = F(XRIGHT(NSTACK),FDERVR)
      LEFTX = SMOOTH + 1
      RIGHTX = LEFTX
C            GET F-VALUES AT OTHER INTERPOLATION POINTS, IF ANY
      IF (NINTRP.EQ.0) GO TO 40
      DO 30 J=1,NINTRP
        FINTRP(J) = F(XINTRP(J),FDUMB)
   30 CONTINUE
C
C          CHECK FOR BREAK POINTS, MODIFY VALUES IF NECESSARY
   40 CONTINUE
      IF (BREAK.NE.LEFT) GO TO 50
      LEFTX = DBREAK(IBREAK) + 1
      FLEFT(LEFTX) = BRIGHT(IBREAK)
   50 IF (BREAK.NE.RIGHT) GO TO 60
      RIGHTX = DBREAK(IBREAK) + 1
      FRIGHT(RIGHTX) = BLEFT(IBREAK)
   60 IF (BREAK.NE.BOTH) GO TO 70
      LEFTX = DBREAK(IBREAK) + 1
      RIGHTX = DBREAK(IBREAK+1) + 1
      FLEFT(LEFTX) = BRIGHT(IBREAK)
      FRIGHT(RIGHTX) = BLEFT(IBREAK+1)
   70 CONTINUE
C
C           COMPUTE DIVIDED DIFFERENCES, NEWTON FORM OF POLYNOMIAL
      CALL NEWTON(LEFTX, RIGHTX, NINTRP)
C
C         COMPUTE NORM OF ERROR OF THIS APPROMIMATION USING FOUR PTS
C         ADD 50 PERCENT FUDGE FACTOR
      ERRORI = ERRINT(F,POLYDD,XLEFT(NSTACK),XRIGHT(NSTACK),ABSC,WGTS)
      ERRORI = 1.5*ERRORI
      RETURN
      END
      REAL FUNCTION ERRINT(F, FIT, AAA, BBB, POINTS, WEIGHT)
C
C        ===============================================================
C
C **  THIS FUNCTION DOES A FOUR POINT INTEGRATION RULE FOR THE
C     ABSOLUTE VALUE OF THE DIFFERENCE OF TWO FUNCTIONS( F AND FIT )
C                   ABS( F(X) - FIT(X) )**NORM
C     THE INTEGRATION USES THE POINTS AND WEIGHTS GIVEN AND SCALED
C     FROM (-1,1) TO (AAA,BBB)
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL AAA, ABMID, BA, BBB, F, FDUMB, FIT, P, PJ, POINTS, WEIGHT
      REAL ER, F1, F2
      DIMENSION FDUMB(5), POINTS(4), WEIGHT(4)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C         COMPUTE MIDPOINT = ABMID AND HALF LENGTH = BA OF INTERVAL
      ABMID = (AAA+BBB)/2.
      BA = (BBB-AAA)/2.
      PJ = ABMID + BA*POINTS(1)
C
C         TEST FOR TCHEBYCHEFF (MINIMAX) NORM WHICH USES SPECIAL CODE
      IF (NORM.EQ.3.) GO TO 20
C
C         HAVE GENERAL LP NORM OR LEAST SQUARES OR LEAST DEVIATIONS
      P = ABS(NORM)
C              INITIALIZE THE QUADRATURE RULE
      ERRINT = ABS(F(PJ,FDUMB)-FIT(PJ))**P*WEIGHT(1)
C              LOOP THROUGH REMAINING POINTS
      DO 10 J=2,4
        PJ = ABMID + BA*POINTS(J)
C             DEBUG DEBUG DEBUG START OF THE COMPUTATION
        F1 = F(PJ,FDUMB)
        F2 = FIT(PJ)
        ER = ABS(F1-F2)**P
        IF (LEVEL.GE.4 .AND. (KNOTS.EQ.1 .OR. FINISH)) WRITE (6,99999)
     *   PJ, F1, F2, ER
        ERRINT = ERRINT + ABS(F(PJ,FDUMB)-FIT(PJ))**P*WEIGHT(J)
   10 CONTINUE
      ERRINT = ERRINT*BA
      GO TO 40
C
C         TCHEBYCHEFF NORM
   20 CONTINUE
C             FIND MAX ERROR ON POINTS
C               INITIALIZE
      ERRINT = ABS(F(PJ,FDUMB)-FIT(PJ))
C              LOOP THROUGH THE REMAINING POINTS
      DO 30 J=2,4
        PJ = ABMID + BA*POINTS(J)
        ERRINT = AMAX1(ERRINT,ABS(F(PJ,FDUMB)-FIT(PJ)))
   30 CONTINUE
   40 CONTINUE
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.GE.3) WRITE (6,99998) ERRINT, AAA, BBB
      RETURN
99999 FORMAT (5(3H --), 31HDEBUG ERROR CURVE,PJ,F1,F2,ER =, 4F15.8)
99998 FORMAT (15X, 9HERRINT = , F20.15, 4H ON , 2F15.8)
      END
      SUBROUTINE NEWTON(NL, NR, NI)
C
C        ===============================================================
C
C **  THIS PROGRAM COMPUTES THE DIVIDED DIFFERENCES ARRAY AS FOLLOWS
C         NL COALESCED POINTS ON LEFT   - DERIV VALUES IN FLEFT
C         NR COALESCED POINTS ON RIGHT  -   -      -    - FRIGHT
C         NI DISTINCT  POINTS INBETWEEN - FNCTN    -    - FINTRP
C
C     THE POINTS ARE ORDERED XL = XLEFT (NSTACK)
C                            XR = XRIGHT(NSTACK)
C                            XINTRP ARRAY
C
C         LAYOUT OF THE DDTEMP DIVIDED DIFFERENCE ARRAY
C
C     NL=6    LLLLLL****II
C     NR=4    LLLLL****II     L = FIRST TRIANGLE
C     NI=2    LLLL****II
C             LLL****II       R = SECOND TRIANGLE
C             LL****II
C             L****II         * = FILL BETWEEN TRIANGLES
C             RRRRII
C             RRRII           I = COMPLETION FOR INTERPOLATION POINTS
C             RRII
C             RII       IDIF = HORIZONTAL COORD. = DIFFERENCE ORDER
C             II        IPT  = VERTICAL COORD. ASSOCIATED WITH POINTS
C             I
C
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL DIFFF, DIFFX
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C                 MAIN CALCULATION OF DIVIDED DIFFERENCES
C         DEFINE A FEW SHORT CONSTANTS
      NL1 = NL - 1
      NL2 = NL + 1
      NR1 = NR - 1
      NR2 = NR + 1
      NRL = NR + NL
C
C         PUT X-VALUES IN A SINGLE ARRAY WITH NDDX = NL+NR+NI POINTS
      DO 10 NDDX=1,NL
        XDD(NDDX) = XLEFT(NSTACK)
   10 CONTINUE
      NDDX = NL
      DO 20 K=1,NR
        NDDX = NDDX + 1
        XDD(NDDX) = XRIGHT(NSTACK)
   20 CONTINUE
C            CHECK IF THERE ARE ANY INTERPOLATION POINTS TO ADD TO XDD
      IF (NI.EQ.0) GO TO 40
      DO 30 K=1,NI
        NDDX = NDDX + 1
        XDD(NDDX) = XINTRP(K)
   30 CONTINUE
C
C           FILL BORDER OF FIRST TRIANGLE - SIZE NL.
   40 CONTINUE
C         TOP BORDER
      DO 50 IDIF=1,NL
        DDTEMP(IDIF,1) = FLEFT(IDIF)/FACTOR(IDIF)
   50 CONTINUE
      IF (NL1.EQ.0) GO TO 70
C                   BOTTOM BORDER
      DO 60 IDIF=1,NL1
        IPT = NL2 - IDIF
        DDTEMP(IDIF,IPT) = DDTEMP(IDIF,1)
   60 CONTINUE
C
C          FILL BORDER OF SECOND TRIANGLE - SIZE NR
   70 CONTINUE
C         TOP BORDER
      DO 80 IDIF=1,NR
        DDTEMP(IDIF,NL2) = FRIGHT(IDIF)/FACTOR(IDIF)
   80 CONTINUE
      IF (NRL.EQ.NL2) GO TO 100
C                   BOTTOM BORDER
      DO 90 IDIF=1,NR1
        IPT = NRL + 1 - IDIF
        DDTEMP(IDIF,IPT) = DDTEMP(IDIF,NL2)
   90 CONTINUE
C
C           FILL PARALLOGRAM BETWEEN THE TWO TRIANGLES JUST FILLED
C          FILL ENTRIES PARALLEL TO BOTTOM OF FIRST TRIANGLE
  100 CONTINUE
C
C         LOOP STEPPING ALONG TOP SIDE OF SECOND TRIANGLE
      DO 120 J=2,NR2
        IDIF = J
C             LOOP STEPPING PARALLEL TO BOTTOM SIDE OF FIRST TRIANGLE
        DO 110 K=2,NL2
          IPT = NL + 2 - K
          DIFFF = DDTEMP(IDIF-1,IPT+1) - DDTEMP(IDIF-1,IPT)
          IPT2 = IPT - 1 + IDIF
          DIFFX = XDD(IPT2) - XDD(IPT)
          DDTEMP(IDIF,IPT) = DIFFF/DIFFX
          IDIF = IDIF + 1
  110   CONTINUE
  120 CONTINUE
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.GE.4 .AND. KNOTS.LE.1) WRITE (6,99999) NR2, NL2, IDIF,
     * IPT, DIFFF, DIFFX
C
C         FILL IN BOTTOM DIAGONALS FOR INTERPOLATION POINTS, IF ANY
      IF (NI.EQ.0) GO TO 150
C         LOOP THROUGH THE INTERPOLATATION POINTS
      DO 140 J=1,NI
        IDIF = 2
        NRLJ = NRL + J
        DDTEMP(1,NRLJ) = FINTRP(J)
C         LOOP THROUGH THE DIFFERENCES (IDIF INDEX)
        NRLJ1 = NRLJ - 1
        DO 130 K=1,NRLJ1
          IPT = NRLJ - K
          DIFFF = DDTEMP(IDIF-1,IPT+1) - DDTEMP(IDIF-1,IPT)
          DIFFX = XDD(NRLJ) - XDD(IPT)
          DDTEMP(IDIF,IPT) = DIFFF/DIFFX
          IDIF = IDIF + 1
  130   CONTINUE
  140 CONTINUE
  150 CONTINUE
      RETURN
99999 FORMAT (9X, 30HNR2,NL2,IDIF,IPT,DIFFF,DIFFX =, 4I3, 2F12.6)
      END
      REAL FUNCTION POLYDD(X)
C
C        ===============================================================
C
C **  THIS FUNCTION EVALUATES THE CURRENT POLYNOMIAL PIECE REPRESENTED
C     BY THE DIVIDED DIFFERENCES DDTEMP ON THE POINT SET XDD.
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL X
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C
      POLYDD = DDTEMP(DEGREE+1,1)
      DO 10 K=1,DEGREE
        J = DEGREE + 1 - K
        POLYDD = DDTEMP(J,1) + (X-XDD(J))*POLYDD
   10 CONTINUE
      RETURN
      END
      SUBROUTINE PPFIT4(F, XLFT, XRGT, EPSLN, NPIECE, ERREST, XKNOTS,
     * COEFS, KDIMEN, NDIMEN, NDEG, NSMTH, EMEAS, LPRNT, FOSCIL, ATYPE,
     * KBREAK, BRAKPT, KDERVB, VALLFT, VALRGT)
C
C        ===============================================================
C
C **  THIS SET OF FOUR CONTROL PROGRAMS SET VARYING NUMBERS OF DEFAULT
C     VALUES FOR ARGUMENTS.  IT USES ENTRY STATEMENTS WHICH ARE DONE
C     DIFFERENTLY BY VARIOUS FORTRANS.  FOR THIS REASON ENTRY STATEMENTS
C     ARE ONLY INDICATED BY COMMENT CARDS.  WRITING FOUR SEPARATE ROU-
C     TINES APPROXIMATELY TRIPLES THE LENGTH OF THE TOTAL CODE.
C
C     THE FOLLOWING TABULATES THE INTERNAL AND EXTERNAL NAMES OF THE
C     ARGUMENTS ALONG WITH THEIR DEFAULT VALUES FOR THE VARIOUS PPFIT.
C     NOTE THAT THIS ALLOWS THE ARGUMENTS TO BE PUT INTO A COMMON BLOCK
C     AND AVOIDS LONG ARGUMENT LISTS WITHIN ADAPT ITSELF.
C
C
C    INTERNAL   EXTERNAL     DEFAULT VALUE    SET BY PPFIT NUMBER
C     F          F              NONE
C     A,B        A,B            NONE
C     ACCUR      EPSLN          NONE
C     KNOTS      NPIECE         OUTPUT
C     ERROR      ERREST         OUTPUT
C     XKNOTS     XKNOTS         OUTPUT
C     COEFS      COEFS          OUTPUT
C     KDIMEN     KDIMEN         NONE
C     NDIMEN     NDIMEN         NONE
C     DEGREE     NDEG           3                1
C     SMOOTH     NSMTH          0                1
C     NORM       EMEAS          3                1 2
C     LEVEL      LPRNT          1                1 2
C     CHARF      FOSCIL         VARIABLE         1 2
C     EDIST      ATYPE          1                1 2 3
C     NBREAK     KBREAK         0                1 2 3
C     XBREAK     BRAKPT         -
C     DBREAK     KDERVB         -
C     BLEFT      VALLFT         -
C     BRIGHT     VALRGT         -
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      REAL BRAKPT, EMEAS, EPSLN, ERREST, FOSCIL, F, VALLFT, VALRGT,
     * XLFT, XRGT
      DIMENSION BRAKPT(20), KDERVB(20), VALLFT(20), VALRGT(20)
      INTEGER ATYPE
      EXTERNAL F
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
      COMMON /SAVEIT/ IKNOT
C
C
C             SKIP ALL DEFAULT SETTING FOR PPFIT4
      GO TO 10
C
C                 SET DEFAULTS FOR PPFIT1
C
C         ****** SINCE ENTRY IS NON-STANDARD  ******
C         ****** THE NATURAL WAY TO IMPLEMENT ******
C         ****** THE OTHER PPFITS IS ONLY     ******
C         ****** INDICATED IN THE COMMENTS    ******
C
C     ENTRY PPFIT1
C     ENTRY PPFIT1(F,XLFT,XRGT,EPSLN,NPIECE,ERREST,XKNOTS,COEFS,
C    A             KDIMEN,NDIMEN)
      NDEG = 3
      NSMTH = 0
C
C                 SET DEFAULTS FOR PPFIT2
C     ENTRY PPFIT2
C     ENTRY PPFIT2(F,XLFT,XRGT,EPSLN,NPIECE,ERREST,XKNOTS,COEFS,
C    A             KDIMEN,NDIMEN,NDEG,NSMTH)
      EMEAS = 3.0
      LPRNT = 1
      FOSCIL = B - A
      IF (NDEG.EQ.2) FOSCIL = .5*FOSCIL
      IF (NDEG.LE.1) FOSCIL = (B-A)/3.
C
C                 SET DEFAULTS FOR PPFIT3
C     ENTRY PPFIT3
C     ENTRY PPFIT3(F,XLFT,XRGT,EPSLN,NPIECE,ERREST,XKNOTS,COEFS,
C    A             KDIMEN,NDIMEN,NDEG,NSMTH,EMEAS,LPRNT,FOSCIL)
      ATYPE = 1
      KSING = 0
      KBREAK = 0
C
C         PUT INPUT INTO COMMON INPUTZ BY CHANGING TO INTERNAL NAMES
   10 A = XLFT
      B = XRGT
      ACCUR = EPSLN
      DEGREE = NDEG
      SMOOTH = NSMTH
      NORM = EMEAS
      LEVEL = LPRNT
      CHARF = FOSCIL
      EDIST = ATYPE
      NBREAK = KBREAK
      IF (NBREAK.LE.0 .OR. NBREAK.GE.21) GO TO 30
      DO 20 K=1,NBREAK
        XBREAK(K) = BRAKPT(K)
        DBREAK(K) = KDERVB(K)
        BLEFT(K) = VALLFT(K)
        BRIGHT(K) = VALRGT(K)
   20 CONTINUE
   30 CONTINUE
C
      CALL ADAPT(F, XKNOTS, COEFS, KDIMEN, NDIMEN)
      NPIECE = KNOTS
      ERREST = ERROR
      RETURN
      END
      REAL FUNCTION PPOLY(T, XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C        ===============================================================
C
C **  THIS FUNCTION EVALUATES THE PIECEWISE POLYNOMIAL APPROXIMATION
C     COMPUTED IN ADAPT
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      REAL T, X
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      COMMON /SAVEIT/ IKNOT
C           START SEARCH FOR RIGHT INTERVAL FROM POINT OF LAST
C           EVALUATION OF PPOLY
C
C           FORWARD SEARCHING LOOP
   10 IF (T.LT.XKNOTS(IKNOT+1)) GO TO 20
      IKNOT = IKNOT + 1
      IF (IKNOT.LT.KNOTS) GO TO 10
      IKNOT = KNOTS - 1
C
C                REACHED RIGHT END OF INTERVAL
      GO TO 30
C
C           BACKWARD SEARCHING LOOP
   20 IF (T.GE.XKNOTS(IKNOT)) GO TO 30
      IKNOT = IKNOT - 1
      IF (IKNOT.GT.1) GO TO 20
      IKNOT = 1
C
C                REACHED LEFT END OF INTERVAL
C
C                EVALUATE FROM POWERS BASED AT XKNOT(IKNOT)
C                USE NESTED MULTIPLICATION
   30 X = T - XKNOTS(IKNOT)
      PPOLY = COEFS(IKNOT,NPAR)
      DO 40 K=1,DEGREE
        KK = NPAR - K
        PPOLY = COEFS(IKNOT,KK) + X*PPOLY
   40 CONTINUE
      RETURN
      END
      SUBROUTINE PTRANS(D, POWERS)
C
C        ===============================================================
C
C **  THIS PROGRAM CONVERTS POLYNOMIAL REPRESENTATION FROM DIVIDED
C     DIFFERENCE TO POWER FORM.  THERE ARE COALESCED POINTS ON EACH
C     END OF THE INTERVAL (XL,XR) = (XLEFT(NSTACK),XRIGHT(NSTACK)).
C     THE NUMBER COALESCED AT EACH END IS LEFTX AND RIGHTX.
C     AND THERE ARE NINTRP OTHER PTS XINTRP(K)  INBETWEEN THEM.
C     SEE SUBROUTINE NEWTON FOR MORE DETAILS
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL D, POWERS, SHIFT, XL, XR, XTEMP
      DIMENSION D(12,12), POWERS(12), XTEMP(12)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C         SET SOME SHORT LOCAL VARIABLE NAMES
      XL = XLEFT(NSTACK)
      XR = XRIGHT(NSTACK)
      NL = LEFTX
      NL1 = NL + 1
      NR = RIGHTX
      NI = NINTRP
      NRL = NR + NL
      NRI = NR + NI
      NRI1 = NRI - 1
      NRLI = NRL + NI
C
C         STARTING REPRESENTATION IS (ASSUMING XL = 0 )
C
C    D(1) +D(2)X +D(3)X**2 + --- +D(NL)X**(NL-1)
C   +(X**NL)*( D(NL+1)(+D(NL+2)(X-XR)**2 + --- +D(NL+NR)*(X-XR)**(NR-1)
C        *((X-XR)**NR)*(D(NL+NR+1) + D(NL+NR+2)*(X-XINTRP(1))
C                      +D(NL+NR+3)*(X-XINTRP(1))(X-XINTRP(2)) + ---))
C
C         STRATEGY IS TO FIRST CONVERT THE PART FROM THE INTERP. PTS.
C         TO POLY IN (X-XR).  THIS POLY THEN HAS ORIGIN SHIFTED TO XL.
C
C     THE CONVERSION OF THE INTERP PART IS DONE EXPLICITLY FOR DEGREE
C     TWO OR LESS AND DONE BY SYNTHETIC DIVISION FOR HIGHER DEGREES
C
C   D1 + D2(X-X1) +D3(X**2-(X1+X2)X +X1*X2)
C
C     THE RESULTING COEFFICIENTS ARE PUT IN THE ARRAY POWERS
C
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.EQ.4) WRITE (6,99999) (D(K,1),K=1,NPAR)
      IF (NI.EQ.0) GO TO 100
C             BUILD UP THE POLYNOMIAL FOR THE INTERPOLATION POINTS
C
C         USE SPECIAL FORMULAS FOR NI LESS THAN 3
      IF (NI.EQ.1) GO TO 10
      IF (NI.EQ.2) GO TO 20
      GO TO 30
   10 POWERS(1) = D(NRL+1,1)
      GO TO 80
   20 POWERS(1) = D(NRL+1,1) + (XR-XINTRP(1))*D(NRL+2,1)
      POWERS(2) = D(NRL+2,1)
      GO TO 80
C
C         CONVERSION BY REPEATED SYNTHETIC DIVISION
   30 NI1 = NI - 1
C         INITIALIZE THE POWERS AND XTEMP ARRAYS
      DO 40 K=1,NI
        XTEMP(K) = XINTRP(K)
        NRLK = NRL + K
        POWERS(K) = D(NRLK,1)
   40 CONTINUE
C
C         DO THE REPEATED SYNTHETIC DIVISION TO REPLACE THE XTEMP
C         = XINTRP POINTS OF THE NEWTON EXPANSION BY THE XR POINTS
      DO 70 K=1,NI1
C              POWERS(NI) IS FIXED AND SET ABOVE
        DO 50 II=1,NI1
          I = NI - II
          POWERS(I) = POWERS(I) + (XR-XTEMP(I))*POWERS(I+1)
   50   CONTINUE
C              SHIFT THE NEWTON EXPANSION PTS. UP, PUT IN ONE MORE XR
        DO 60 II=1,NI1
          I = NI - II
          XTEMP(I+1) = XTEMP(I)
   60   CONTINUE
        XTEMP(1) = XR
   70 CONTINUE
      IF (LEVEL.EQ.4) WRITE (6,99998) (POWERS(K),K=1,NI)
   80 CONTINUE
C             SHIFT THE COEFFICIENTS TO THE TOP OF THE POWERS ARRAY
      DO 90 K=1,NI
        L = NI + 1 - K
        LTOP = L + NRL
        POWERS(LTOP) = POWERS(L)
   90 CONTINUE
C
C             HAVE THE INTERPOLATION PT. COEFS. IN THE ARRAY POWERS
  100 CONTINUE
C             PUT THE REMAINING DIVIDED DIFFS INTO THE POWERS ARRAY
      DO 110 J=1,NRL
        POWERS(J) = D(J,1)
  110 CONTINUE
C
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.EQ.4) WRITE (6,99997) (POWERS(K),K=NL1,NPAR)
C         TRANSFORM THE ORIGIN OF THE POLYNOMIAL FROM XR TO XL
C         WE USE REPEATED SYNTHETIC DIVISION
      IF (NRI.EQ.1) GO TO 140
      SHIFT = XR - XL
      KHI = NRI1
C             LOOP THROUGH THE COEFFICIENTS
      DO 130 J=2,NRI
C                   SYNTHETIC DIVISION LOOP
        DO 120 K=1,KHI
          KOEF = NRLI - K
          POWERS(KOEF) = POWERS(KOEF) - SHIFT*POWERS(KOEF+1)
  120   CONTINUE
        KHI = KHI - 1
  130 CONTINUE
  140 CONTINUE
C     THE COEFFICIENTS ARE NOW OF THE POWER FORM WITH ORIGIN XL
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.EQ.4) WRITE (6,99996) (POWERS(K),K=1,NPAR)
      RETURN
99999 FORMAT (15X, 26HDEBUG PTRANS, ORIG INPUT =/(5X, 8E15.5))
99998 FORMAT (10X, 17HINTERP PART COEFS/(5X, 8E15.5))
99997 FORMAT (10X, 25HRIGHT + INTERP PART COEFS/(5X, 8E15.5))
99996 FORMAT (10X, 18HFINAL COEFFICIENTS/(5X, 8E15.5))
      END
      SUBROUTINE PUT(INTERV, XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C        ===============================================================
C **  THIS PROGRAM PUTS INTERVALS ON THE STACK OR DISCARDS THEM.
C     WHEN AN INTERVAL IS DISCARDED A NEW KNOT IS FOUND. THEN THIS
C     PROGRAM UPDATES THE ERROR ESTIMATE, THE XKNOT ARRAY, TRANSFORMS
C     THE POLYNOMIAL TO THE POWER FORM AND PUT THE COEFFICIENTS INTO
C     THE ARRAY COEFS.  IT ALSO CHECKS FOR PASSING BREAK POINTS
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      REAL BUFFER, DX, POWERS, P, RATIO
      DIMENSION POWERS(12)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      DATA BUFFER /1.E-12/
C            CHECK FOR DISCARDING THE INTERVAL
      IF (DISCRD) GO TO 30
C         SUBDIVIDE INTERVAL AND PLACE ON STACK
      IF (NSTACK.LT.MAXSTK) GO TO 10
C             FATAL ERROR, EXCEEDED ACTIVE STACK SIZE
      FATAL = .TRUE.
      IF (LEVEL.GE.0) WRITE (6,99999) FMESGE, MAXSTK, XLEFT(NSTACK),
     * XRIGHT(NSTACK)
      DISCRD = .TRUE.
      GO TO 30
C
   10 DX = (XRIGHT(NSTACK)-XLEFT(NSTACK))*.5
C              CHECK FOR SMALL INTERVALS
      RATIO = DX/(ABS(A)+ABS(B))
      IF (RATIO.GT.BUFFER) GO TO 20
      DISCRD = .TRUE.
      FATAL = .TRUE.
      IF (LEVEL.GE.0) WRITE (6,99998) XLEFT(NSTACK), XRIGHT(NSTACK)
      GO TO 30
   20 CONTINUE
      NSTACK = NSTACK + 1
      XLEFT(NSTACK) = XLEFT(NSTACK-1)
      XLEFT(NSTACK-1) = XRIGHT(NSTACK-1) - DX
      XRIGHT(NSTACK) = XLEFT(NSTACK-1)
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.GE.3) WRITE (6,99997) NSTACK, XLEFT(NSTACK),
     * XRIGHT(NSTACK)
      RETURN
C
C            DISCARD INTERVAL, UPDATE GLOBAL ERROR, XKNOTS AND COEFS
   30 CONTINUE
C
      P = ABS(NORM)
      IF (NORM.EQ.3.) ERROR = AMAX1(ERROR,ERRORI)
      IF (NORM.NE.3.) ERROR = (ERROR**P+ERRORI)**(1./P)
C
C              CHECK FOR PASSING BREAK POINTS
      IF (BREAK.EQ.LEFT .OR. BREAK.EQ.BOTH) IBREAK = IBREAK + 1
C             DEBUG  DEBUG  DEBUG  DEBUG
      IF (LEVEL.GE.3) WRITE (6,99996) NSTACK
C
C             TRANSFORM REPRESENTATION OF POLYNOMIAL FROM DIVIDED
C             DIFFERENCES TO POWERS OF X WITH ORIGIN AT XKNOTS (KNOTS)
      CALL PTRANS(DDTEMP, POWERS)
C
C             PUT COEFS INTO THE MAIN ARRAY
      DO 40 K=1,NPAR
        COEFS(KNOTS,K) = POWERS(K)
   40 CONTINUE
C            PUT THE NEW KNOTS IN XKNOTS
      KNOTS = KNOTS + 1
      XKNOTS(KNOTS) = XRIGHT(NSTACK)
      NSTACK = NSTACK - 1
      RETURN
99999 FORMAT (//2X, 6A4, 36HINTERVAL DIVIDED TOO MUCH, EXCEEDED ,
     * 6HLIMIT , I3, 22H ON INTERVAL STACK AT /20X, 2E18.8/2X, 6HINTERV,
     * 38HAL DISCARDED AND COMPUTATION CONTINUED)
99998 FORMAT (25X, 23HGOT SHORT INTERVAL ****, 2E22.13, 12H **** DISCAR,
     * 4HD IT)
99997 FORMAT (15X, 13HPUT INTERVAL , I3, 10H ON STACK , 2F12.8)
99996 FORMAT (15X, 16HDISCARD INTERVAL, I5)
      END
      SUBROUTINE SETUP(XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C        ===============================================================
C
C **  THIS PROGRAM CHECKS THE INPUT DATA AND INITIALIZES THE COMPUTATION
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      INTEGER HMSGE(6), NMSGE(4,4)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C              IKNOT IS A COUNTER USED IN THE OUTPUT APPROX. PPOLY
      COMMON /SAVEIT/ IKNOT
      DATA HMSGE(1), HMSGE(2), HMSGE(3), HMSGE(4), HMSGE(5), HMSGE(6) /
     * 4H ***,4H* FA,4HTAL ,4HERRO,4HR **,4H**  /
      DATA NMSGE(1,1), NMSGE(1,2), NMSGE(1,3), NMSGE(1,4) /4HLEAS,
     * 4HT DE,4HVIAT,4HIONS/
      DATA NMSGE(2,1), NMSGE(2,2), NMSGE(2,3), NMSGE(2,4) /4HLEAS,
     * 4HT SQ,4HUARE,4HS   /
      DATA NMSGE(3,1), NMSGE(3,2), NMSGE(3,3), NMSGE(3,4) /4HMINI,
     * 4HMAX ,4HNORM,4H    /
      DATA NMSGE(4,1), NMSGE(4,2), NMSGE(4,3), NMSGE(4,4) /4HGENE,
     * 4HRAL ,4HL-P ,4HNORM/
      DATA KLEFT, KRIGHT, KBOTH /4HLEFT,4HRITE,4HBOTH/
C
C         PUT DATA STATEMENT ITEMS INTO COMMON VARIABLES
C         VARIOUS COMPILERS ALLOW MORE EFFICIENT WAYS
C
      LEFT = KLEFT
      RIGHT = KRIGHT
      BOTH = KBOTH
      DO 10 K=1,6
        FMESGE(K) = HMSGE(K)
   10 CONTINUE
C -------- SET CURRENT VALUES OF LIMITS ON DIMENSIONS ------------------
      KNTDIM = KDIMEN
      NPARDM = NDIMEN
      MAXKNT = KNTDIM
      MAXSTK = 50
      MAXPAR = MIN0(12,NPARDM)
      MAXAUX = 20
C -------- CHECK INPUT DATA --------------------------------------------
      FATAL = .FALSE.
      IF ((LEVEL+1)*LEVEL*(LEVEL-1)*(LEVEL-2).EQ.0) GO TO 20
C         FIXED ERROR ON PRINT CONTROL LEVEL
      WRITE (6,99999) LEVEL
C         FOR DEBUGGING DO NOT CHANGE OUTPUT LEVEL
C         FOR PRODUCTION VERSION SET LEVEL = 0
   20 IF (A.LT.B) GO TO 30
C         FATAL ERROR IN INTERVAL (A,B)
      FATAL = .TRUE.
      WRITE (6,99998) FMESGE, A, B
   30 IF (ACCUR.GT.0.0) GO TO 40
C         FATAL ERROR, NEGATIVE ACCURACY
      FATAL = .TRUE.
      WRITE (6,99997) FMESGE, ACCUR
   40 IF (DEGREE.LT.MAXPAR) GO TO 50
C         FATAL ERROR, DEGREE EXCEEDS MAXIMUM ALLOWED VALUE
      IF (LEVEL.GE.0) WRITE (6,99996) FMESGE, DEGREE
      FATAL = .TRUE.
   50 IF (2*SMOOTH.LT.DEGREE) GO TO 60
C         FATAL ERROR, SMOOTH AND DEGREE INCOMPATIBLE
      FATAL = .TRUE.
      WRITE (6,99995) FMESGE, SMOOTH, DEGREE
   60 IF (CHARF.GT.(B-A)/FLOAT(MAXKNT)) GO TO 70
C         FATAL ERROR, CHARF IS TOO SMALL
      FATAL = .TRUE.
      WRITE (6,99994) FMESGE, CHARF
   70 IF ((NORM-1.)*(NORM-2.)*(NORM-3.).EQ.0.0 .OR. NORM.LT.0.0) GO TO
     * 80
C         FATAL ERROR, UNDEFINED NORM REQUESTED
      IF (LEVEL.GE.0) WRITE (6,99993) FMESGE, NORM
      FATAL = .TRUE.
   80 IF (NBREAK.GE.0 .AND. NBREAK.LE.MAXAUX) GO TO 90
C         FATAL ERROR IN NUMBER OF BREAK POINTS
      FATAL = .TRUE.
      WRITE (6,99992) FMESGE, NBREAK
C
   90 IF (NBREAK.EQ.0) GO TO 150
C         CHECK THE BREAK POINT DATA, MONOTONICITY AND DEGREE
      J = 1
      IF (XBREAK(1).LT.A .OR. XBREAK(NBREAK).GT.B) GO TO 110
      IF (NBREAK.EQ.1) GO TO 120
      DO 100 J=2,NBREAK
        IF (XBREAK(J-1).GE.XBREAK(J)) GO TO 110
  100 CONTINUE
      GO TO 120
  110 FATAL = .TRUE.
C         BREAK POINTS ARE NOT MONOTONIC
      WRITE (6,99991) FMESGE, XBREAK(J)
  120 LIMSM = (DEGREE-1)/2
      DO 130 J=1,NBREAK
        IF (DBREAK(J).LT.0 .OR. DBREAK(J).GT.LIMSM) GO TO 140
  130 CONTINUE
      GO TO 150
C         BAD VALUE IN DERIVATIVE BREAKS
  140 FATAL = .TRUE.
      WRITE (6,99990) FMESGE, DBREAK(J)
C
  150 IF (EDIST*(EDIST-1)*(EDIST-2).EQ.0) GO TO 160
C         FATAL ERROR, ILLEGAL ERROR DISTRIBUTION REQUESTED
      IF (LEVEL.GE.0) WRITE (6,99989) FMESGE, EDIST
      FATAL = .TRUE.
  160 CONTINUE
C
C -------- INITIALIZATION OF VARIABLES ---------------------------------
C
C             ACTIVE INTERVAL STACK
      NSTACK = 1
      XLEFT(1) = A
      XRIGHT(1) = B
C             TERMINATION AND ERROR VALUES
      FINISH = .FALSE.
      ERROR = 0.
      DSCTOL = ACCUR**ABS(NORM)
      IF (EDIST.EQ.0) DSCTOL = DSCTOL/(B-A)
      IF (NORM.EQ.3.) DSCTOL = ACCUR
C             MISCELLANEOUS VARIABLES AND POINTERS
      IBREAK = 1
      KNOTS = 1
      IKNOT = 1
      INTERP = DEGREE + 2 - 2*SMOOTH
      XKNOTS(1) = A
      NPAR = DEGREE + 1
C
C             COMPUTE ARRAY OF NPAR FACTORIALS, FACTOR(K)= K-1 FACTORIAL
      FACTOR(1) = 1.
      FACTOR(2) = 1.
      DO 170 K=3,NPAR
        FACTOR(K) = FLOAT(K-1)*FACTOR(K-1)
  170 CONTINUE
C
C -------- PRINT OUT OF PROBLEM TO BE SOLVED ---------------------------
      IF (LEVEL.LE.0) GO TO 180
      NMES = NORM
      IF (NORM.LT.0.0) NMES = 4
      WRITE (6,99988) A, B, DEGREE, SMOOTH, ACCUR, (NMSGE(NMES,J),J=1,4)
      WRITE (6,99987) CHARF, NORM
      IF (EDIST.EQ.0) WRITE (6,99986)
      IF (EDIST.EQ.1) WRITE (6,99985)
      IF (EDIST.EQ.2) WRITE (6,99984)
      IF (EDIST.EQ.-1) WRITE (6,99983)
      IF (NBREAK.EQ.0) GO TO 180
      WRITE (6,99982) NBREAK
      WRITE (6,99981) (XBREAK(J),DBREAK(J),BLEFT(J),BRIGHT(J),J=1,
     * NBREAK)
  180 CONTINUE
      RETURN
99999 FORMAT (//2X, 29HILLEGAL OUTPUT LEVEL CONTROL , I2, 9H SET TO 0)
99998 FORMAT (//2X, 6A4, 20H INCORRECT INTERVAL , 2E15.5)
99997 FORMAT (//2X, 6A4, 28H ILLEGAL ACCURACY REQUESTED , E15.5)
99996 FORMAT (//2X, 6A4, 6HDEGREE, I3, 29H EXCEEDS MAXIMUM ALLOWED VALU,
     * 1HE)
99995 FORMAT (//2X, 6A4, 35H INCOMPATIBLE DEGREE AND SMOOTHNESS, 2I5)
99994 FORMAT (//2X, 6A4, 40H CHARACTERISTIC OSCILLATION LENGTH OF F ,
     * E15.5, 14H IS TOO SMALL )
99993 FORMAT (//2X, 6A4, 20H ERROR MEASURE NORM , E15.5, 11H NOT ALLOWE,
     * 1HD)
99992 FORMAT (//2X, 6A4, 11H IN NUMBER , I4, 17H OF BREAK POINTS )
99991 FORMAT (//2X, 6A4, 30H BREAK POINTS NOT IN ORDER AT , E15.5)
99990 FORMAT (//2X, 6A4, 20H ILLEGAL DERIVATIVE , I3, 10H AT BREAK )
99989 FORMAT (//2X, 6A4, 32HILLEGAL ERROR DISTRIBUTION TYPE , I3)
99988 FORMAT (//2X, 45H++++++ PIECEWISE POLYNOMIAL APPROXIMATION ON ,
     * 9HINTERVAL , 2E15.4, 11H OF DEGREE , I2, 6H WITH , I2, 7H CONTIN,
     * 17HUOUS DERIVATIVES /10X, 23H ACCURACY REQUESTED IS , E15.4,
     * 13H MEASURED BY , 4A4)
99987 FORMAT (10X, 43HOTHER INPUT/DEFAULT VARIABLES ARE FOSCIL = ,
     * E15.4, 10H, EMEAS = , E15.4)
99986 FORMAT (15X, 9(2H--), 33H PROPORTIONAL ERROR DISTRIBUTION )
99985 FORMAT (15X, 9(2H--), 36HAPPROXIMATE FIXED ERROR DISTRIBUTION)
99984 FORMAT (15X, 9(2H--), 25H FIXED ERROR DISTRIBUTION)
99983 FORMAT (15X, 9(2H--), 28HMODIFIED PROPORTIONAL ERROR , 8HDISTRIBU,
     * 4HTION)
99982 FORMAT (I12, 36H BREAK POINTS SPECIFIED WITH DATA = , 9H LOCATION,
     * 32H, DERIVATIVE, LEFT+RIGHT VALUES )
99981 FORMAT (5X, 2(E20.5, I4, E16.5, E15.5))
      END
      SUBROUTINE SUMARY(XKNOTS, COEFS, KDIMEN, NDIMEN)
C
C        ===============================================================
C
C **  THIS PROGRAM PRINTS OUT A SUMMARY OF RESULTS OF ADAPT
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN), COEFS(KDIMEN,NDIMEN)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      IF (LEVEL.EQ.-1) RETURN
C                   LEVEL 0 OUTPUT
      WRITE (6,99999) DEGREE, SMOOTH, KNOTS, ERROR
      IF (LEVEL.EQ.0) RETURN
C
C                   LEVEL 1 OUTPUT - KNOTS AND COEFFICIENTS
      KNOTS1 = KNOTS - 1
      IF (KNOTS.GT.15) GO TO 20
C
C               FORMAT FOR SMALL NO. OF KNOTS
      WRITE (6,99998)
      DO 10 K=1,KNOTS1
        WRITE (6,99997) K, XKNOTS(K), COEFS(K,1)
        WRITE (6,99996) (COEFS(K,J),J=2,NPAR)
   10 CONTINUE
      GO TO 40
C
C               FORMAT FOR LOTS OF KNOTS
   20 CONTINUE
      WRITE (6,99995)
      DO 30 K=1,KNOTS1
        WRITE (6,99994) K, XKNOTS(K), (COEFS(K,J),J=1,NPAR)
   30 CONTINUE
   40 RETURN
99999 FORMAT (///48H --- ADAPTIVE PIECEWISE POLYNOMIAL APPROXIMATION,
     * 11H OF DEGREE , I2, 6H WITH , I2, 12H CONTINUOUS , 10HDERIVATIVE,
     * 8HS NEEDED, I4, 17H KNOTS FOR ERROR , E10.4)
99998 FORMAT (8X, 13HKNOT LOCATION, 13X, 21HX-POWER COEFFICIENTS ,
     * 27HRELATIVE TO KNOT LOCATIONS )
99997 FORMAT (I10, 2E20.12)
99996 FORMAT (30X, E20.12)
99995 FORMAT (3X, 39HK   K-TH INTERIOR KNOT     POWERS OF X , 7HRELATIV,
     * 20HE TO KNOT LOCATIONS )
99994 FORMAT (I5, E25.12, 4E22.12/(E30.12, 4E22.12))
      END
      SUBROUTINE TAKE(INTERV)
C
C        ===============================================================
C
C **  THIS PROGRAM TAKES AN ACTIVE INTERVAL OFF THE TOP OF THE STACK
C     IT ALSO DOES MOST OF THE WORK OF LOCATING AND HANDLING
C     BREAK POINTS
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL DX, RATIO, BUFFER
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
      DATA BUFFER /1.E-12/
C
C         CHECK FOR BREAK POINT
      BREAK = 0
      IF (NBREAK.EQ.0 .OR. IBREAK.GT.NBREAK) GO TO 20
      IF (XBREAK(IBREAK).GT.XRIGHT(NSTACK)) GO TO 20
C
C               SET CONTROL VARIABLE BREAK, CHECK FOR LOCATION
      IF (XBREAK(IBREAK).GT.XLEFT(NSTACK)) GO TO 10
      BREAK = LEFT
      IF (IBREAK.EQ.NBREAK) GO TO 20
C             CHECK FOR SECOND BREAK POINT IN THIS INTERVAL
      IF (XBREAK(IBREAK+1).GE.XRIGHT(NSTACK)) GO TO 20
C                 NEXT BREAK IS INSIDE INTERVAL, SPLIT TOP INTERVAL
      BREAK = BOTH
C                  CHECK EXCEEDING STACK LIMIT. IF SO, STOP
      IF (NSTACK.EQ.MAXSTK) GO TO 30
C                     DONT SPLIT VERY SMALL INTERVALS, STOP INSTEAD
      DX = XBREAK(IBREAK+1) - XLEFT(NSTACK)
      RATIO = DX/(ABS(A)+ABS(B))
      IF (RATIO.LE.BUFFER) GO TO 40
      NSTACK = NSTACK + 1
      XLEFT(NSTACK) = XLEFT(NSTACK-1)
      XRIGHT(NSTACK) = XBREAK(IBREAK+1)
      XLEFT(NSTACK-1) = XRIGHT(NSTACK)
      IF (LEVEL.GE.2) WRITE (6,99999) NSTACK, XLEFT(NSTACK),
     * XRIGHT(NSTACK)
      GO TO 20
C
   10 BREAK = RIGHT
C                 CHECK TO SEE IF BREAK IS ALREADY AT RIGHT END POINT
      IF (XBREAK(IBREAK).GE.XRIGHT(NSTACK)) GO TO 20
C                 THE BREAK IS INSIDE INTERVAL, SPLIT TOP INTERVAL
C                  CHECK EXCEEDING STACK LIMIT. IF SO, STOP
      IF (NSTACK.EQ.MAXSTK) GO TO 30
C                     DONT SPLIT VERY SMALL INTERVALS, STOP INSTEAD
      DX = XBREAK(IBREAK) - XLEFT(NSTACK)
      RATIO = DX/(ABS(A)+ABS(B))
      IF (RATIO.LE.BUFFER) GO TO 40
      NSTACK = NSTACK + 1
      XLEFT(NSTACK) = XLEFT(NSTACK-1)
      XRIGHT(NSTACK) = XBREAK(IBREAK)
      XLEFT(NSTACK-1) = XRIGHT(NSTACK)
      IF (LEVEL.GE.2) WRITE (6,99998) NSTACK, XLEFT(NSTACK),
     * XRIGHT(NSTACK)
   20 CONTINUE
C             DEBUG  DEBUG  DEBUG
      IF (LEVEL.GE.3) WRITE (6,99997) XLEFT(NSTACK), XRIGHT(NSTACK),
     * BREAK
      RETURN
C
C        UNABLE TO PROCEED BECAUSE THE STACK LIMIT IS REACHED WITH
C        MULTIPLE BREAKPOINTS INSIDE THE SMALLEST ONE.
   30 FATAL = .TRUE.
      FINISH = .TRUE.
      IF (LEVEL.GE.0) WRITE (6,99996) FMESGE, MAXSTK, XLEFT(NSTACK),
     * XRIGHT(NSTACK), XBREAK(IBREAK)
      RETURN
C
C         SPLITTING INTERVALS TO ACCOMODATE BREAK POINTS HAS LED TO
C         AN EXCESSIVELY SMALL INTERVAL
   40 FATAL = .TRUE.
      FINISH = .TRUE.
      IF (LEVEL.GE.0) WRITE (6,99995) FMESGE, XLEFT(NSTACK),
     * XBREAK(IBREAK), XRIGHT(NSTACK)
      RETURN
99999 FORMAT (10X, 28H++++ SPLIT TOP INTERVAL, GET, I3, 2E18.8,
     * 17H FOR BREAK = BOTH)
99998 FORMAT (10X, 28H++++ SPLIT TOP INTERVAL, GET, I3, 2E18.8,
     * 19H FOR BREAK = RIGHT )
99997 FORMAT (15X, 14HTAKE INTERVAL , 2F12.8, 11H    BREAK =, A4)
99996 FORMAT (//2X, 6A4, 42H INTERVAL DIVIDED TOO MUCH, EXCEEDED LIMIT,
     * I3, 21H ON INTERVAL STACK AT/20X, 2E18.8/2X, 12HBREAK POINT ,
     * E18.8, 48H PRESENT REQUIRES SPLITTING INTERVAL, STOP ADAPT)
99995 FORMAT (//2X, 6A4, 43HBREAK POINT ADJUSTMENT HAS LED TO A FATALLY,
     * 1H , 40HSMALL INTERVAL, THE POINTS INVOLVED ARE /20X, 3E18.8)
      END
      SUBROUTINE TERMIN(TEST, AND, PRINT, XKNOTS, KDIMEN)
C
C        ===============================================================
C
C **  THIS PROGRAM TESTS FOR TERMINATION AND GIVES INTERMEDIATE OUTPUT
C
      REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR,
     * ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD,
     * XINTRP, XLEFT, XRIGHT
      DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20)
      DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6)
      DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6),
     * XDD(12), XINTRP(10)
      INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT,
     * RIGHTX, SMOOTH
      LOGICAL DISCRD, FATAL, FINISH
      REAL XKNOTS(KDIMEN)
      COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT,
     * DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM
      COMMON /RESULZ/ ERROR, KNOTS
      COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH,
     * FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR,
     * MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH
      COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP,
     * LEFTX, NINTRP, RIGHTX
C
C         INTERMEDIATE OUTPUT - ACCORDING TO LEVELS
C         NO INTERMEDIATE OUTPUT FOR LEVELS -1,0,1
      IF (LEVEL.LE.1) GO TO 40
C
C             LEVEL 2 OUTPUT - ONLY FOR DISCARDED INTERVAL
      IF (.NOT.DISCRD) GO TO 10
      WRITE (6,99999) KNOTS, XKNOTS(KNOTS), ERRORI, ERROR
      IF (BREAK.EQ.LEFT) WRITE (6,99998)
      IF (BREAK.EQ.RIGHT) WRITE (6,99997)
      IF (BREAK.EQ.BOTH) WRITE (6,99996)
C
      GO TO 20
C             DEBUG OUTPUT  -  LEVEL 3
   10 IF (LEVEL.EQ.2) GO TO 40
C
C                  INTERVAL SUMMARY
      WRITE (6,99995) NSTACK, XLEFT(NSTACK), XRIGHT(NSTACK), ERRORI
      IF (BREAK.NE.0) WRITE (6,99994) IBREAK, XBREAK(IBREAK)
C
C             DEBUG POLYNOMIAL OPERATIONS  -  LEVEL 4
   20 CONTINUE
      IF (LEVEL.LE.3) GO TO 40
      WRITE (6,99993) LEFTX, RIGHTX, NINTRP
      DO 30 K=1,NPAR
        KNPAR = NPAR + 1 - K
        WRITE (6,99992) K, (DDTEMP(I,K),I=1,KNPAR)
   30 CONTINUE
   40 CONTINUE
C             TEST FOR NORMAL TERMINATION
      IF (NSTACK.EQ.0) GO TO 50
C
C             TEST FOR ABNORMAL TERMINATION
      IF (KNOTS.LT.MAXKNT) RETURN
C
C            HAVE EXCEEDED LIMIT ON KNOTS
      WRITE (6,99991) FMESGE, MAXKNT, XKNOTS(KNOTS), ERROR
      FATAL = .TRUE.
C
C         TERMINATE COMPUTATION
   50 FINISH = .TRUE.
      RETURN
99999 FORMAT (2X, 9H**** KNOT, I4, 4H AT , E16.5, 17H, WITH LOCAL AND ,
     * 16HGLOBAL ERRORS = , 2E15.4)
99998 FORMAT (20X, 32HBREAK POINT ON LEFT OF INTERVAL )
99997 FORMAT (20X, 33HBREAK POINT ON RIGHT OF INTERVAL )
99996 FORMAT (20X, 37HBREAK POINT AT BOTH ENDS OF INTERVAL )
99995 FORMAT (10X, 19H++++ STACK ELEMENT , I3, 3H = , E20.5, E15.5,
     * 19H, WITH LOCAL ERROR , E15.4, 17H PLACED ON STACK )
99994 FORMAT (15X, 14HBBBREAK POINT , I3, 3H = , E15.4, 9H INVOLVED)
99993 FORMAT (5X, 44H---- DIVIDED DIFFERENCE INFO ---- NL,NR,NI =, 3I3)
99992 FORMAT (5X, 12HDIV DIFF ROW, I3, 9F12.5/93X, 3F12.5)
99991 FORMAT (//2X, 6A4, 14HEXCEEDED LIMIT, I3, 14H ON NUMBER OF ,
     * 9HKNOTS AT , E18.8, 13H WITH ERROR =, E14.4)
      END
 
