C DAPM 10
C ..................................................................DAPM 20
C DAPM 30
C SUBROUTINE DAPMM DAPM 40
C DAPM 50
C PURPOSE DAPM 60
C APPROXIMATE A FUNCTION TABULATED IN N POINTS BY ANY LINEAR DAPM 70
C COMBINATION OF M GIVEN CONTINUOUS FUNCTIONS IN THE SENSE DAPM 80
C OF CHEBYSHEV. DAPM 90
C DAPM 100
C USAGE DAPM 110
C CALL DAPMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER) DAPM 120
C PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT IN THE DAPM 130
C CALLING PROGRAM. DAPM 140
C DAPM 150
C DESCRIPTION OF PARAMETERS DAPM 160
C FCT - NAME OF SUBROUTINE TO BE SUPPLIED BY THE USER. DAPM 170
C IT COMPUTES VALUES OF M GIVEN FUNCTIONS FOR DAPM 180
C ARGUMENT VALUE X. DAPM 190
C USAGE DAPM 200
C CALL FCT(Y,X,K) DAPM 210
C DESCRIPTION OF PARAMETERS DAPM 220
C Y - DOUBLE PRECISION RESULT VECTOR OF DIMEN- DAPM 230
C SION M CONTAINING THE VALUES OF GIVEN DAPM 240
C CONTINUOUS FUNCTIONS FOR GIVEN ARGUMENT X DAPM 250
C X - DOUBLE PRECISON ARGUMENT VALUE DAPM 260
C K - AN INTEGER VALUE WHICH IS EQUAL TO M-1 DAPM 270
C REMARKS DAPM 280
C IF APPROXIMATION BY NORMAL CHEBYSHEV, SHIFTED DAPM 290
C CHEBYSHEV, LEGENDRE, LAGUERRE, HERMITE POLYNO- DAPM 300
C MIALS IS DESIRED SUBROUTINES DCNP,DCSP,DLEP, DAPM 310
C DLAP,DHEP, RESPECTIVELY FROM SSP COULD BE USED. DAPM 320
C N - NUMBER OF DATA POINTS DEFINING THE FUNCTION WHICH DAPM 330
C IS TO BE APPROXIMATED DAPM 340
C M - NUMBER OF GIVEN CONTINUOUS FUNCTIONS FROM WHICH DAPM 350
C THE APPROXIMATING FUNCTION IS CONSTRUCTED. DAPM 360
C TOP - DOUBLE PRECISION VECTOR OF DIMENSION 3*N. DAPM 370
C ON ENTRY IT MUST CONTAIN FROM TOP(1) UP TO TOP(N) DAPM 380
C THE GIVEN N FUNCTION VALUES AND FROM TOP(N+1) UP DAPM 390
C TO TOP(2*N) THE CORRESPONDING NODES DAPM 400
C ON RETURN TOP CONTAINS FROM TOP(1) UP TO TOP(N) DAPM 410
C THE ERRORS AT THOSE N NODES. DAPM 420
C OTHER VALUES OF TOP ARE SCRATCH. DAPM 430
C IHE - INTEGER VECTOR OF DIMENSION 3*M+4*N+6 DAPM 440
C PIV - DOUBLE PRECISION VECTOR OF DIMENSION 3*M+6. DAPM 450
C ON RETURN PIV CONTAINS AT PIV(1) UP TO PIV(M) THE DAPM 460
C RESULTING COEFFICIENTS OF LINEAR APPROXIMATION. DAPM 470
C T - DOUBLE PRECISION AUXILIARY VECTOR OF DIMENSION DAPM 480
C (M+2)*(M+2) DAPM 490
C ITER - RESULTANT INTEGER WHICH SPECIFIES THE NUMBER OF DAPM 500
C ITERATIONS NEEDED DAPM 510
C IER - RESULTANT ERROR PARAMETER CODED IN THE FOLLOWING DAPM 520
C FORM DAPM 530
C IER=0 - NO ERROR DAPM 540
C IER=1 - THE NUMBER OF ITERATIONS HAS REACHED DAPM 550
C THE INTERNAL MAXIMUM N+M DAPM 560
C IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARA- DAPM 570
C METER M OR N OR SINCE AT SOME ITERATION DAPM 580
C NO SUITABLE PIVOT COULD BE FOUND DAPM 590
C DAPM 600
C REMARKS DAPM 610
C NO ACTION BESIDES ERROR MESSAGE IN CASE M LESS THAN 1 OR DAPM 620
C N LESS THAN 2. DAPM 630
C DAPM 640
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DAPM 650
C THE EXTERNAL SUBROUTINE FCT MUST BE FURNISHED BY THE USER. DAPM 660
C DAPM 670
C METHOD DAPM 680
C THE PROBLEM OF APPROXIMATION A TABULATED FUNCTION BY ANY DAPM 690
C LINEAR COMBINATION OF GIVEN FUNCTIONS IN THE SENSE OF DAPM 700
C CHEBYSHEV (I.E. TO MINIMIZE THE MAXIMUM ERROR) IS TRANS- DAPM 710
C FORMED INTO A LINEAR PROGRAMMING PROBLEM. DAPMM USES A DAPM 720
C REVISED SIMPLEX METHOD TO SOLVE A CORRESPONDING DUAL DAPM 730
C PROBLEM. FOR REFERENCE, SEE DAPM 740
C I.BARRODALE/A.YOUNG, ALGORITHMS FOR BEST L-SUB-ONE AND DAPM 750
C L-SUB-INFINITY, LINEAR APPROXIMATIONS ON A DISCRETE SET, DAPM 760
C NUMERISCHE MATHEMATIK, VOL.8, ISS.3 (1966), PP.295-306. DAPM 770
C DAPM 780
C ..................................................................DAPM 790
C DAPM 800
SUBROUTINE DAPMM(FCT,N,M,TOP,IHE,PIV,T,ITER,IER) DAPM 810
C DAPM 820
C DAPM 830
DIMENSION TOP(1),IHE(1),PIV(1),T(1) DAPM 840
DOUBLE PRECISION DSUM,TOP,PIV,T,SAVE,HELP,REPI,TOL DAPM 850
C DAPM 860
C TEST ON WRONG INPUT PARAMETERS N AND M DAPM 870
IER=-1 DAPM 880
IF (N-1) 81,81,1 DAPM 890
1 IF(M) 81,81,2 DAPM 900
C DAPM 910
C INITIALIZE CHARACTERISTIC VECTORS FOR THE TABLEAU DAPM 920
2 IER=0 DAPM 930
C DAPM 940
C PREPARE TOP-ROW TOP DAPM 950
DO 3 I=1,N DAPM 960
K=I+N DAPM 970
J=K+N DAPM 980
TOP(J)=TOP(K) DAPM 990
3 TOP(K)=-TOP(I) DAPM1000
C DAPM1010
C PREPARE INVERSE TRANSFORMATION MATRIX T DAPM1020
L=M+2 DAPM1030
LL=L*L DAPM1040
DO 4 I=1,LL DAPM1050
4 T(I)=0.D0 DAPM1060
K=1 DAPM1070
J=L+1 DAPM1080
DO 5 I=1,L DAPM1090
T(K)=1.D0 DAPM1100
5 K=K+J DAPM1110
C DAPM1120
C PREPARE INDEX-VECTOR IHE DAPM1130
DO 6 I=1,L DAPM1140
K=I+L DAPM1150
J=K+L DAPM1160
IHE(I)=0 DAPM1170
IHE(K)=I DAPM1180
6 IHE(J)=1-I DAPM1190
NAN=N+N DAPM1200
K=L+L+L DAPM1210
J=K+NAN DAPM1220
DO 7 I=1,NAN DAPM1230
K=K+1 DAPM1240
IHE(K)=I DAPM1250
J=J+1 DAPM1260
7 IHE(J)=I DAPM1270
C DAPM1280
C SET COUNTER ITER FOR ITERATION-STEPS DAPM1290
ITER=-1 DAPM1300
8 ITER=ITER+1 DAPM1310
C DAPM1320
C TEST FOR MAXIMUM ITERATION-STEPS DAPM1330
IF(N+M-ITER) 9,9,10 DAPM1340
9 IER=1 DAPM1350
GO TO 69 DAPM1360
C DAPM1370
C DETERMINE THE COLUMN WITH THE MOST POSITIVE ELEMENT IN TOP DAPM1380
10 ISE=0 DAPM1390
IPIV=0 DAPM1400
K=L+L+L DAPM1410
SAVE=0.D0 DAPM1420
C DAPM1430
C START TOP-LOOP DAPM1440
DO 14 I=1,NAN DAPM1450
IDO=K+I DAPM1460
HELP=TOP(I) DAPM1470
IF(HELP-SAVE) 12,12,11 DAPM1480
11 SAVE=HELP DAPM1490
IPIV=I DAPM1500
12 IF(IHE(IDO)) 14,13,14 DAPM1510
13 ISE=I DAPM1520
14 CONTINUE DAPM1530
C END OF TOP-LOOP DAPM1540
C DAPM1550
C IS OPTIMAL TABLEAU REACHED DAPM1560
IF(IPIV) 69,69,15 DAPM1570
C DAPM1580
C DETERMINE THE PIVOT-ELEMENT FOR THE COLUMN CHOSEN UPOVE DAPM1590
15 ILAB=1 DAPM1600
IND=0 DAPM1610
J=ISE DAPM1620
IF(J) 21,21,34 DAPM1630
C DAPM1640
C TRANSFER K-TH COLUMN FROM T TO PIV DAPM1650
16 K=(K-1)*L DAPM1660
DO 17 I=1,L DAPM1670
J=L+I DAPM1680
K=K+1 DAPM1690
17 PIV(J)=T(K) DAPM1700
C DAPM1710
C IS ANOTHER COLUMN NEEDED FOR SEARCH FOR PIVOT-ELEMENT DAPM1720
18 IF(ISE) 22,22,19 DAPM1730
19 ISE=-ISE DAPM1740
C DAPM1750
C TRANSFER COLUMNS IN PIV DAPM1760
J=L+1 DAPM1770
IDO=L+L DAPM1780
DO 20 I=J,IDO DAPM1790
K=I+L DAPM1800
20 PIV(K)=PIV(I) DAPM1810
21 J=IPIV DAPM1820
GO TO 34 DAPM1830
C DAPM1840
C SEARCH PIVOT-ELEMENT PIV(IND) DAPM1850
22 SAVE=1.D38 DAPM1860
IDO=0 DAPM1870
K=L+1 DAPM1880
LL=L+L DAPM1890
IND=0 DAPM1900
C DAPM1910
C START PIVOT-LOOP DAPM1920
DO 29 I=K,LL DAPM1930
J=I+L DAPM1940
HELP=PIV(I) DAPM1950
IF(HELP) 29,29,23 DAPM1960
23 HELP=-HELP DAPM1970
IF(ISE) 26,24,26 DAPM1980
24 IF(IHE(J)) 27,25,27 DAPM1990
25 IDO=I DAPM2000
GO TO 29 DAPM2010
26 HELP=-PIV(J)/HELP DAPM2020
27 IF(HELP-SAVE) 28,29,29 DAPM2030
28 SAVE=HELP DAPM2040
IND=I DAPM2050
29 CONTINUE DAPM2060
C END OF PIVOT-LOOP DAPM2070
C DAPM2080
C TEST FOR SUITABLE PIVOT-ELEMENT DAPM2090
IF(IND) 30,30,32 DAPM2100
30 IF(IDO) 68,68,31 DAPM2110
31 IND=IDO DAPM2120
C PIVOT-ELEMENT IS STORED IN PIV(IND) DAPM2130
C DAPM2140
C COMPUTE THE RECIPROCAL OF THE PIVOT-ELEMENT REPI DAPM2150
32 REPI=1.D0/PIV(IND) DAPM2160
IND=IND-L DAPM2170
C DAPM2180
C UPDATE THE TOP-ROW TOP OF THE TABLEAU DAPM2190
ILAB=0 DAPM2200
SAVE=-TOP(IPIV)*REPI DAPM2210
TOP(IPIV)=SAVE DAPM2220
C DAPM2230
C INITIALIZE J AS COUNTER FOR TOP-LOOP DAPM2240
J=NAN DAPM2250
33 IF(J-IPIV) 34,53,34 DAPM2260
34 K=0 DAPM2270
C DAPM2280
C SEARCH COLUMN IN TRANSFORMATION-MATRIX T DAPM2290
DO 36 I=1,L DAPM2300
IF(IHE(I)-J) 36,35,36 DAPM2310
35 K=I DAPM2320
IF(ILAB) 50,50,16 DAPM2330
36 CONTINUE DAPM2340
C DAPM2350
C GENERATE COLUMN USING SUBROUTINE FCT AND TRANSFORMATION-MATRIX DAPM2360
I=L+L+L+NAN+J DAPM2370
I=IHE(I)-N DAPM2380
IF(I) 37,37,38 DAPM2390
37 I=I+N DAPM2400
K=1 DAPM2410
38 I=I+NAN DAPM2420
C DAPM2430
C CALL SUBROUTINE FCT DAPM2440
CALL FCT(PIV,TOP(I),M-1) DAPM2450
C DAPM2460
C PREPARE THE CALLED VECTOR PIV DAPM2470
DSUM=0.D0 DAPM2480
IDO=M DAPM2490
DO 41 I=1,M DAPM2500
HELP=PIV(IDO) DAPM2510
IF(K) 39,39,40 DAPM2520
39 HELP=-HELP DAPM2530
40 DSUM=DSUM+HELP DAPM2540
PIV(IDO+1)=HELP DAPM2550
41 IDO=IDO-1 DAPM2560
PIV(L)=-DSUM DAPM2570
PIV(1)=1.D0 DAPM2580
C DAPM2590
C TRANSFORM VECTOR PIV WITH ROWS OF MATRIX T DAPM2600
IDO=IND DAPM2610
IF(ILAB) 44,44,42 DAPM2620
42 K=1 DAPM2630
43 IDO=K DAPM2640
44 DSUM=0.D0 DAPM2650
HELP=0.D0 DAPM2660
C DAPM2670
C START MULTIPLICATION-LOOP DAPM2680
DO 46 I=1,L DAPM2690
DSUM=DSUM+PIV(I)*T(IDO) DAPM2700
TOL=DABS(DSUM) DAPM2710
IF(TOL-HELP) 46,46,45 DAPM2720
45 HELP=TOL DAPM2730
46 IDO=IDO+L DAPM2740
C END OF MULTIPLICATION-LOOP DAPM2750
C DAPM2760
TOL=1.D-14*HELP DAPM2770
IF(DABS(DSUM)-TOL) 47,47,48 DAPM2780
47 DSUM=0.D0 DAPM2790
48 IF(ILAB) 51,51,49 DAPM2800
49 I=K+L DAPM2810
PIV(I)=DSUM DAPM2820
C DAPM2830
C TEST FOR LAST COLUMN-TERM DAPM2840
K=K+1 DAPM2850
IF(K-L) 43,43,18 DAPM2860
50 I=(K-1)*L+IND DAPM2870
DSUM=T(I) DAPM2880
C DAPM2890
C COMPUTE NEW TOP-ELEMENT DAPM2900
51 DSUM=DSUM*SAVE DAPM2910
TOL=1.D-14*DABS(DSUM) DAPM2920
TOP(J)=TOP(J)+DSUM DAPM2930
IF(DABS(TOP(J))-TOL) 52,52,53 DAPM2940
52 TOP(J)=0.D0 DAPM2950
C DAPM2960
C TEST FOR LAST TOP-TERM DAPM2970
53 J=J-1 DAPM2980
IF(J) 54,54,33 DAPM2990
C END OF TOP-LOOP DAPM3000
C DAPM3010
C TRANSFORM PIVOT-COLUMN DAPM3020
54 I=IND+L DAPM3030
PIV(I)=-1.D0 DAPM3040
DO 55 I=1,L DAPM3050
J=I+L DAPM3060
55 PIV(I)=-PIV(J)*REPI DAPM3070
C DAPM3080
C UPDATE TRANSFORMATION-MATRIX T DAPM3090
J=0 DAPM3100
DO 57 I=1,L DAPM3110
IDO=J+IND DAPM3120
SAVE=T(IDO) DAPM3130
T(IDO)=0.D0 DAPM3140
DO 56 K=1,L DAPM3150
ISE=K+J DAPM3160
56 T(ISE)=T(ISE)+SAVE*PIV(K) DAPM3170
57 J=J+L DAPM3180
C DAPM3190
C UPDATE INDEX-VECTOR IHE DAPM3200
C INITIALIZE CHARACTERISTICS DAPM3210
J=0 DAPM3220
K=0 DAPM3230
ISE=0 DAPM3240
IDO=0 DAPM3250
C DAPM3260
C START QUESTION-LOOP DAPM3270
DO 61 I=1,L DAPM3280
LL=I+L DAPM3290
ILAB=IHE(LL) DAPM3300
IF(IHE(I)-IPIV) 59,58,59 DAPM3310
58 ISE=I DAPM3320
J=ILAB DAPM3330
59 IF(ILAB-IND) 61,60,61 DAPM3340
60 IDO=I DAPM3350
K=IHE(I) DAPM3360
61 CONTINUE DAPM3370
C END OF QUESTION-LOOP DAPM3380
C DAPM3390
C START MODIFICATION DAPM3400
IF(K) 62,62,63 DAPM3410
62 IHE(IDO)=IPIV DAPM3420
IF(ISE) 67,67,65 DAPM3430
63 IF(IND-J) 64,66,64 DAPM3440
64 LL=L+L+L+NAN DAPM3450
K=K+LL DAPM3460
I=IPIV+LL DAPM3470
ILAB=IHE(K) DAPM3480
IHE(K)=IHE(I) DAPM3490
IHE(I)=ILAB DAPM3500
IF(ISE) 67,67,65 DAPM3510
65 IDO=IDO+L DAPM3520
I=ISE+L DAPM3530
IHE(IDO)=J DAPM3540
IHE(I)=IND DAPM3550
66 IHE(ISE)=0 DAPM3560
67 LL=L+L DAPM3570
J=LL+IND DAPM3580
I=LL+L+IPIV DAPM3590
ILAB=IHE(I) DAPM3600
IHE(I)=IHE(J) DAPM3610
IHE(J)=ILAB DAPM3620
C END OF MODIFICATION DAPM3630
C DAPM3640
GO TO 8 DAPM3650
C DAPM3660
C SET ERROR PARAMETER IER=-1 SINCE NO SUITABLE PIVOT IS FOUND DAPM3670
68 IER=-1 DAPM3680
C DAPM3690
C EVALUATE FINAL TABLEAU DAPM3700
C COMPUTE SAVE AS MAXIMUM ERROR OF APPROXIMATION AND DAPM3710
C HELP AS ADDITIVE CONSTANCE FOR RESULTING COEFFICIENTS DAPM3720
69 SAVE=0.D0 DAPM3730
HELP=0.D0 DAPM3740
K=L+L+L DAPM3750
DO 73 I=1,NAN DAPM3760
IDO=K+I DAPM3770
J=IHE(IDO) DAPM3780
IF(J) 71,70,73 DAPM3790
70 SAVE=-TOP(I) DAPM3800
71 IF(M+J+1) 73,72,73 DAPM3810
72 HELP=TOP(I) DAPM3820
73 CONTINUE DAPM3830
C DAPM3840
C PREPARE T,TOP,PIV DAPM3850
T(1)=SAVE DAPM3860
IDO=NAN+1 DAPM3870
J=NAN+N DAPM3880
DO 74 I=IDO,J DAPM3890
74 TOP(I)=SAVE DAPM3900
DO 75 I=1,M DAPM3910
75 PIV(I)=HELP DAPM3920
C DAPM3930
C COMPUTE COEFFICIENTS OF RESULTING POLYNOMIAL IN PIV(1) UP TO PIDAPM3940
C AND CALCULATE ERRORS AT GIVEN NODES IN TOP(1) UP TO TOP(N) DAPM3950
DO 79 I=1,NAN DAPM3960
IDO=K+I DAPM3970
J=IHE(IDO) DAPM3980
IF(J) 76,79,77 DAPM3990
76 J=-J DAPM4000
PIV(J)=HELP-TOP(I) DAPM4010
GO TO 79 DAPM4020
77 IF(J-N) 78,78,79 DAPM4030
78 J=J+NAN DAPM4040
TOP(J)=SAVE+TOP(I) DAPM4050
79 CONTINUE DAPM4060
DO 80 I=1,N DAPM4070
IDO=NAN+I DAPM4080
80 TOP(I)=TOP(IDO) DAPM4090
81 RETURN DAPM4100
END DAPM4110