home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
High Voltage Shareware
/
high1.zip
/
high1
/
DIR9
/
ACMALG01.ZIP
/
ACM525.FOR
< prev
next >
Wrap
Text File
|
1991-02-15
|
90KB
|
2,256 lines
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