-
Notifications
You must be signed in to change notification settings - Fork 1
/
525.f
2252 lines (2252 loc) · 91.2 KB
/
525.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C MAN 10
C **** THIS PROGRAM RUNS ADAPT WITH DATA READ IN FROM CARDS MAN 20
C IT IS DESIGNED TO TEST ADAPT WITH THE 20 FUNCTIONS IN F(X). MAN 30
C IT READS DATA FOR EVERYTHING EXCEPT EDIST(=ATYPE) AND PRINTS MAN 40
C A HEADING WITH THE FUNCTIONS NAME. IT CALLS AND TIMES ADAPT, MAN 50
C THEN INDEPENDENTLY CHECKS THE ERROR IN ADAPT AND PLOTS THE MAN 60
C FUNCTION AND ERROR CURVE. MAN 70
C MAN 80
C NOTES -- USER MUST SUPPLY SYSTEM ROUTINES SECOND AND GRAPH MAN 90
C INTEGRATION USES ERRINT FROM ADAPT. MAN 100
C PRESENT DIMENSIONS (ON ZKNOTS + ZCOEFS) PREVENT USING MAN 110
C ERRINT TO CHECK ACCURACY FOR MORE THAN 100 KNOTS OR MAN 120
C POLYNOMIAL DEGREE 6. MAN 130
C QPOLY USED BY ERRINT IS A SINGLE ARGUMENT VERSION OF PPOLY. MAN 140
C SEE COMMENTS IN ADAPT ABOUT DIMENSION CONSTRAINTS AND MAN 150
C PORTABILITY. MAN 160
C IF ADAPT AND THIS DRIVER ARE CHANGED TO DOUBLE PRECISION MAN 170
C ONE PROBABLY WANTS TO LEAVE TIME,TSTART,TSTOP,XVALU, MAN 180
C GRAF1,GRAF2 AS SINGLE PRECISION MAN 190
C THEY ARE SEPARATELY DECLARED HERE. MAN 200
C MAN 210
C MAN 220
REAL A, ACCUR, B, BLEFT, BRIGHT, CHARF, DDTEMP, DSCTOL, ERROR, MAN 230
* ERRORI, FACTOR, FINTRP, FLEFT, FRIGHT, NORM, XBREAK, XDD, MAN 240
* XINTRP, XLEFT, XRIGHT MAN 250
DIMENSION XBREAK(20), DBREAK(20), BLEFT(20), BRIGHT(20) MAN 260
DIMENSION XLEFT(50), XRIGHT(50), FACTOR(12), FMESGE(6) MAN 270
DIMENSION DDTEMP(12,12), FINTRP(10), FLEFT(6), FRIGHT(6), MAN 280
* XDD(12), XINTRP(10) MAN 290
INTEGER BOTH, BREAK, DBREAK, DEGREE, EDIST, FMESGE, RIGHT, MAN 300
* RIGHTX, SMOOTH MAN 310
LOGICAL DISCRD, FATAL, FINISH MAN 320
DIMENSION XKNOTS(140), COEFS(140,12) MAN 330
REAL COEFS, ECHECK, F, XKNOTS MAN 340
REAL TIME, TSTART MAN 350
DIMENSION NAMES(25,10) MAN 360
EXTERNAL F MAN 370
COMMON /INPUTZ/ A, B, ACCUR, NORM, CHARF, XBREAK, BLEFT, BRIGHT, MAN 380
* DBREAK, DEGREE, SMOOTH, LEVEL, EDIST, NBREAK, KNTDIM, NPARDM MAN 390
COMMON /RESULZ/ ERROR, KNOTS MAN 400
COMMON /KONTRL/ DSCTOL, ERRORI, XLEFT, XRIGHT, BREAK, BOTH, MAN 410
* FACTOR, FMESGE, IBREAK, INTERP, LEFT, MAXAUX, MAXKNT, MAXPAR, MAN 420
* MAXSTK, NPAR, NSTACK, RIGHT, DISCRD, FATAL, FINISH MAN 430
COMMON /COMDIF/ DDTEMP, FINTRP, FLEFT, FRIGHT, XDD, XINTRP, MAN 440
* LEFTX, NINTRP, RIGHTX MAN 450
COMMON /SAVEIT/ IKNOT MAN 460
C MAN 470
COMMON /FDATA/ JFUNK, KOUNT, MAXK MAN 480
COMMON /INFO/ ECHECK, KFUNK, TIME, TSTART MAN 490
DATA NAMES(1,1), NAMES(1,2), NAMES(1,3), NAMES(1,4), NAMES(1,5), MAN 500
* NAMES(1,6), NAMES(1,7), NAMES(1,8), NAMES(1,9), NAMES(1,10) / MAN 510
* 4HSIMP,4HLE A,4HLGEB,4HRAIC,4H SIN,4HGULA,4HRITY,4H AT ,4H0.0 , MAN 520
* 4H / MAN 530
DATA NAMES(2,1), NAMES(2,2), NAMES(2,3), NAMES(2,4), NAMES(2,5), MAN 540
* NAMES(2,6), NAMES(2,7), NAMES(2,8), NAMES(2,9), NAMES(2,10) / MAN 550
* 4HINTE,4HRIOR,4H SIN,4HGULA,4HRITY,4H AT ,4HX = ,4H.307,4H7070, MAN 560
* 4H7707/ MAN 570
DATA NAMES(3,1), NAMES(3,2), NAMES(3,3), NAMES(3,4), NAMES(3,5), MAN 580
* NAMES(3,6), NAMES(3,7), NAMES(3,8), NAMES(3,9), NAMES(3,10) / MAN 590
* 4HBROK,4HEN L,4HINE ,4H- BR,4HEAKS,4H AT ,4H1,2,,4H3.5 ,4HAND , MAN 600
* 4H4 / MAN 610
DATA NAMES(4,1), NAMES(4,2), NAMES(4,3), NAMES(4,4), NAMES(4,5), MAN 620
* NAMES(4,6), NAMES(4,7), NAMES(4,8), NAMES(4,9), NAMES(4,10) / MAN 630
* 4HHUMP,4H AT ,4H2.5 ,4HRECI,4HPROC,4HAL O,4HF QU,4HARTI,4HC , MAN 640
* 4H / MAN 650
DATA NAMES(5,1), NAMES(5,2), NAMES(5,3), NAMES(5,4), NAMES(5,5), MAN 660
* NAMES(5,6), NAMES(5,7), NAMES(5,8), NAMES(5,9), NAMES(5,10) / MAN 670
* 4HSIMP,4HLE C,4HUBIC,4H = (,4H4X-2,4H)**3,4H - (,4H4X-2,4H) , MAN 680
* 4H / MAN 690
DATA NAMES(6,1), NAMES(6,2), NAMES(6,3), NAMES(6,4), NAMES(6,5), MAN 700
* NAMES(6,6), NAMES(6,7), NAMES(6,8), NAMES(6,9), NAMES(6,10) / MAN 710
* 4HVARI,4HABLE,4H OSC,4HILLA,4HTION,4H FRO,4HM SI,4HN(X*,4H*2-1, MAN 720
* 4H) / MAN 730
DATA NAMES(7,1), NAMES(7,2), NAMES(7,3), NAMES(7,4), NAMES(7,5), MAN 740
* NAMES(7,6), NAMES(7,7), NAMES(7,8), NAMES(7,9), NAMES(7,10) / MAN 750
* 4HEXPO,4HNENT,4HIAL ,4H ,4H ,4H ,4H ,4H ,4H , MAN 760
* 4H / MAN 770
DATA NAMES(8,1), NAMES(8,2), NAMES(8,3), NAMES(8,4), NAMES(8,5), MAN 780
* NAMES(8,6), NAMES(8,7), NAMES(8,8), NAMES(8,9), NAMES(8,10) / MAN 790
* 4HFIFT,4HH DE,4HGREE,4H POL,4HY PL,4HUS F,4HINE ,4HSAW ,4HTOOT, MAN 800
* 4HH / MAN 810
DATA NAMES(9,1), NAMES(9,2), NAMES(9,3), NAMES(9,4), NAMES(9,5), MAN 820
* NAMES(9,6), NAMES(9,7), NAMES(9,8), NAMES(9,9), NAMES(9,10) / MAN 830
* 4HSEVE,4HNTH ,4HDEGR,4HEE P,4HOLYN,4HOMIA,4HL ,4H ,4H , MAN 840
* 4H / MAN 850
DATA NAMES(10,1), NAMES(10,2), NAMES(10,3), NAMES(10,4), MAN 860
* NAMES(10,5), NAMES(10,6), NAMES(10,7), NAMES(10,8), NAMES(10,9), MAN 870
* NAMES(10,10) /4HWILD,4H VAR,4HIATI,4HON B,4HETWE,4HEN .,4H47 A, MAN 880
* 4HND .,4H48 ,4H / MAN 890
DATA NAMES(11,1), NAMES(11,2), NAMES(11,3), NAMES(11,4), MAN 900
* NAMES(11,5), NAMES(11,6), NAMES(11,7), NAMES(11,8), NAMES(11,9), MAN 910
* NAMES(11,10) /4HLOTS,4H OF ,4HOSCI,4HLLAT,4HION ,4H - F,4HIRST, MAN 920
* 4H EXA,4HMPLE,4H / MAN 930
DATA NAMES(12,1), NAMES(12,2), NAMES(12,3), NAMES(12,4), MAN 940
* NAMES(12,5), NAMES(12,6), NAMES(12,7), NAMES(12,8), NAMES(12,9), MAN 950
* NAMES(12,10) /4HLOTS,4H OF ,4HOSCI,4HLLAT,4HION ,4H - S,4HECON, MAN 960
* 4HD EX,4HAMPL,4HE / MAN 970
DATA NAMES(13,1), NAMES(13,2), NAMES(13,3), NAMES(13,4), MAN 980
* NAMES(13,5), NAMES(13,6), NAMES(13,7), NAMES(13,8), NAMES(13,9), MAN 990
* NAMES(13,10) /4HSING,4HULAR,4HITIE,4HS AT,4H BOT,4HH EN,4HD PO, MAN 1000
* 4HINTS,4H ,4H / MAN 1010
DATA NAMES(14,1), NAMES(14,2), NAMES(14,3), NAMES(14,4), MAN 1020
* NAMES(14,5), NAMES(14,6), NAMES(14,7), NAMES(14,8), NAMES(14,9), MAN 1030
* NAMES(14,10) /4HRAND,4HOM N,4HUMBE,4HRS A,4HDDED,4H TO ,4HSIN(, MAN 1040
* 4HX) ,4H ,4H / MAN 1050
DATA NAMES(15,1), NAMES(15,2), NAMES(15,3), NAMES(15,4), MAN 1060
* NAMES(15,5), NAMES(15,6), NAMES(15,7), NAMES(15,8), NAMES(15,9), MAN 1070
* NAMES(15,10) /4HLARG,4HE FU,4HNCTI,4HON -,4H ADJ,4HUSTA,4HBLE , MAN 1080
* 4HSIZE,4H ,4H / MAN 1090
DATA NAMES(16,1), NAMES(16,2), NAMES(16,3), NAMES(16,4), MAN 1100
* NAMES(16,5), NAMES(16,6), NAMES(16,7), NAMES(16,8), NAMES(16,9), MAN 1110
* NAMES(16,10) /4HSHIF,4HTED ,4HORIG,4HIN -,4H ADJ,4HUSTA,4HBLE , MAN 1120
* 4HSHIF,4HT ,4H / MAN 1130
DATA NAMES(17,1), NAMES(17,2), NAMES(17,3), NAMES(17,4), MAN 1140
* NAMES(17,5), NAMES(17,6), NAMES(17,7), NAMES(17,8), NAMES(17,9), MAN 1150
* NAMES(17,10) /4HNUME,4HRICA,4HL DI,4HFFER,4HENTI,4HATIO,4HN OF, MAN 1160
* 4H F(X,4H) ,4H / MAN 1170
DATA NAMES(18,1), NAMES(18,2), NAMES(18,3), NAMES(18,4), MAN 1180
* NAMES(18,5), NAMES(18,6), NAMES(18,7), NAMES(18,8), NAMES(18,9), MAN 1190
* NAMES(18,10) /4HCOMP,4HLICA,4HTED ,4HFUNC,4HTION,4H WIT,4HH 7 , MAN 1200
* 4HPIEC,4HES ,4H / MAN 1210
DATA NAMES(19,1), NAMES(19,2), NAMES(19,3), NAMES(19,4), MAN 1220
* NAMES(19,5), NAMES(19,6), NAMES(19,7), NAMES(19,8), NAMES(19,9), MAN 1230
* NAMES(19,10) /4HJUMP,4H DIS,4HCONT,4HINUI,4HTY A,4HT VA,4HRIAB, MAN 1240
* 4HLE P,4HOINT,4H / MAN 1250
DATA NAMES(20,1), NAMES(20,2), NAMES(20,3), NAMES(20,4), MAN 1260
* NAMES(20,5), NAMES(20,6), NAMES(20,7), NAMES(20,8), NAMES(20,9), MAN 1270
* NAMES(20,10) /4HINFI,4HNITE,4H SIN,4HGULA,4HRITY,4H AT ,4HP10A, MAN 1280
* 4H ,4H ,4H / MAN 1290
C MAN 1300
10 READ (5,99999) JFUNK, SMOOTH, DEGREE, LEVEL, NORM, A, B, ACCUR, MAN 1310
* CHARF, NBREAK MAN 1320
IF (JFUNK.EQ.0) STOP MAN 1330
IF (CHARF.EQ.0.) CHARF = (B-A)*1.001 MAN 1340
IF (NBREAK.EQ.0) GO TO 20 MAN 1350
READ (5,99998) (DBREAK(K),XBREAK(K),BLEFT(K),BRIGHT(K),K=1,NBREAK)MAN 1360
20 CONTINUE MAN 1370
MAXK = 3200 MAN 1380
EDIST = 0 MAN 1390
WRITE (6,99997) JFUNK, (NAMES(JFUNK,K),K=1,10) MAN 1400
KOUNT = 0 MAN 1410
C MAN 1420
C ***** CALL ON CLOCK TO INITIALIZE TSTART IN SECONDS ***** MAN 1430
C MAN 1440
TSTART = 0.0 MAN 1450
CALL ADAPT(F, XKNOTS, COEFS, 140, 12) MAN 1460
CALL SUMRY2(XKNOTS, COEFS, 140, 12) MAN 1470
GO TO 10 MAN 1480
99999 FORMAT (4I3, F4.3, 2F10.3, E12.4, F10.3, 2I3) MAN 1490
99998 FORMAT (2(I5, 3F11.5)) MAN 1500
99997 FORMAT (//5(3H+++), 18HADAPT FOR FUNCTION, I3, 3H = , 10A4, 4X, MAN 1510
* 5(3H+++)) MAN 1520
END MAN 1530
BLOCK DATA BLK 10
C SET PARAMETERS FOR FUNCTIONS IN F BLK 20
REAL PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A, P6B, P7, BLK 30
* P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B, P10C BLK 40
COMMON /FPARS/ PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B, P6A, BLK 50
* P6B, P7, P8A, P8B, P8C, P8D, P8E, P8F, P8G, P8H, P9, P10A, P10B, BLK 60
* P10C BLK 70
DATA PAR2, PAR3A, PAR3B, PAR4, PAR4A, P5A, P5B BLK 80
* /.02,.4,.6,.001,6.,8000.,2037./ BLK 90
DATA P6A, P6B, P7, P8A, P8B, P8C, P8D, P8E /4000.,3998.,.0001,-2.,BLK 100
* 1.5,7.,5.,9.6/ BLK 110
DATA P8F, P8G, P8H, P9, P10A, P10B, P10C /.6,14.,16.,1.0,.25,.5, BLK 120
* .6667/ BLK 130
END BLK 140
REAL FUNCTION F(X, FDERV) F 10
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) QPO 10
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) SUM 10
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(100), COEFS(100,100)
REAL ABSC(4), DX, DXG, ECHECK, ERRINT, ERRP, P, PPOLY, QPOLY,
* WGTS(4), XL, XR, F
REAL XVALU(100), GRAF1(100), GRAF2(100), TIME, TSTART, TSTOP,
* ZKNOTS(100), ZCOEF(100,7), DUMMY1(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
C GRAF1(K) = F(XVALU(K),DUMMY)
GRAF1(K) = (XVALU(K))
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) ADA 10
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