C DARA 10
C ..................................................................DARA 20
C DARA 30
C SUBROUTINE DARAT DARA 40
C DARA 50
C PURPOSE DARA 60
C CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE DARA 70
C FUNCTION IN THE LEAST SQUARES SENSE DARA 80
C DARA 90
C USAGE DARA 100
C CALL DARAT(DATI,N,WORK,P,IP,IQ,IER) DARA 110
C DARA 120
C DESCRIPTION OF PARAMETERS DARA 130
C DATI - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS DARA 140
C THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS, DARA 150
C THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND DARA 160
C THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY. DARA 170
C IF NO WEIGHTS ARE TO BE USED THEN THE THIRD DARA 180
C COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT DARA 190
C WHICH MUST CONTAIN A NONPOSITIVE VALUE DARA 200
C DATI MUST BE OF DOUBLE PRECISION DARA 210
C N - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION DARA 220
C WORK - WORKING STORAGE WHICH IS OF DIMENSION DARA 230
C (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST. DARA 240
C ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED DARA 250
C IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF DARA 260
C THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO DARA 270
C WORK(3*N) DARA 280
C WORK MUST BE OF DOUBLE PRECISION DARA 290
C P - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND DARA 300
C NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ DARA 310
C LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP DARA 320
C LOCATIONS. DARA 330
C COEFFICIENTS ARE ORDERED FROM LOW TO HIGH. DARA 340
C P MUST BE OF DOUBLE PRECISION DARA 350
C IP - DIMENSION OF THE NUMERATOR (INPUT VALUE) DARA 360
C IQ - DIMENSION OF THE DENOMINATOR (INPUT VALUE) DARA 370
C IER - RESULTANT ERROR PARAMETER DARA 380
C IER =-1 MEANS FORMAL ERRORS DARA 390
C IER = 0 MEANS NO ERRORS DARA 400
C IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION DARA 410
C IER IS ALSO USED AS INPUT VALUE DARA 420
C A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN DARA 430
C INITIAL APPROXIMATION STORED IN P DARA 440
C DARA 450
C REMARKS DARA 460
C THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR DARA 470
C OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P DARA 480
C STARTING WITH LOW POWERS (DENOMINATOR FIRST). DARA 490
C IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE. DARA 500
C SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL DARA 510
C FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL DARA 520
C (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEARDARA 530
C TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS. DARA 540
C IF A FIT IN OTHER FUNCTIONS IS REQUIRED, DCNP AND DCNPS MUSTDARA 550
C BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN. DARA 560
C DARA 570
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DARA 580
C DAPLL, DAPFS, DFRAT, DCNPS, DCNP DARA 590
C DCNP IS REQUIRED WITHIN DFRAT DARA 600
C DARA 610
C METHOD DARA 620
C THE ITERATIVE SCHEME USED FOR CALCULATION OF THE DARA 630
C APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS DARA 640
C WHICH ARE OBTAINED BY LINEARIZATION. DARA 650
C A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH DARA 660
C IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF DARA 670
C ZEROES WITHIN THE APPROXIMATION INTERVAL. DARA 680
C FOR REFERENCE SEE DARA 690
C D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN, DARA 700
C COMPUTING(1966), VOL.1, ED.3, PP.264-272. DARA 710
C D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION DARA 720
C OF NONLINEAR PARAMETERS, DARA 730
C JSIAM(1963), VOL.11, ED.2, PP.431-441. DARA 740
C DARA 750
C ..................................................................DARA 760
C DARA 770
SUBROUTINE DARAT(DATI,N,WORK,P,IP,IQ,IER) DARA 780
C DARA 790
C DARA 800
EXTERNAL DFRAT DARA 810
C DARA 820
C DIMENSIONED LOCAL VARIABLE DARA 830
DIMENSION IERV(3) DARA 840
C DARA 850
C DIMENSIONED DUMMY VARIABLES DARA 860
DIMENSION DATI(1),WORK(1),P(1) DARA 870
DOUBLE PRECISION DATI,WORK,P,T,OSUM,DIAG,RELAX,SUM,SSOE,SAVE DARA 880
C DARA 890
C INITIALIZE TESTVALUES DARA 900
LIMIT=20 DARA 910
ETA=1.E-29 DARA 920
EPS=1.E-14 DARA 930
C DARA 940
C CHECK FOR FORMAL ERRORS DARA 950
IF(N)4,4,1 DARA 960
1 IF(IP)4,4,2 DARA 970
2 IF(IQ)4,4,3 DARA 980
3 IPQ=IP+IQ DARA 990
IF(N-IPQ)4,5,5 DARA1000
C DARA1010
C ERROR RETURN IN CASE OF FORMAL ERRORS DARA1020
4 IER=-1 DARA1030
RETURN DARA1040
C DARA1050
C INITIALIZE ITERATION PROCESS DARA1060
5 KOUNT=0 DARA1070
IERV(2)=IP DARA1080
IERV(3)=IQ DARA1090
NDP=N+N+1 DARA1100
NNE=NDP+NDP DARA1110
IX=IPQ-1 DARA1120
IQP1=IQ+1 DARA1130
IRHS=NNE+IPQ*IX/2 DARA1140
IEND=IRHS+IX DARA1150
C DARA1160
C TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION DARA1170
IF(IER)8,6,8 DARA1180
C DARA1190
C INITIALIZE NUMERATOR AND DENOMINATOR DARA1200
6 DO 7 I=2,IPQ DARA1210
7 P(I)=0.D0 DARA1220
P(1)=1.D0 DARA1230
C DARA1240
C CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL DARA1250
C APPROXIMATION DARA1260
8 DO 9 J=1,N DARA1270
T=DATI(J) DARA1280
I=J+N DARA1290
CALL DCNPS(WORK(I),T,P(IQP1),IP) DARA1300
K=I+N DARA1310
9 CALL DCNPS(WORK(K),T,P,IQ) DARA1320
C DARA1330
C SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION) DARA1340
10 CALL DAPLL(DFRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV) DARA1350
C DARA1360
C CHECK FOR ZERO DENOMINATOR DARA1370
IF(IERV(1))4,11,4 DARA1380
11 INCR=0 DARA1390
RELAX=2.D0 DARA1400
C DARA1410
C RESTORE MATRIX IN WORKING STORAGE DARA1420
12 J=IEND DARA1430
DO 13 I=NNE,IEND DARA1440
J=J+1 DARA1450
13 WORK(I)=WORK(J) DARA1460
IF(KOUNT)14,14,15 DARA1470
C DARA1480
C SAVE SQUARE SUM OF ERRORS DARA1490
14 OSUM=WORK(IEND) DARA1500
DIAG=OSUM*EPS DARA1510
K=IQ DARA1520
C DARA1530
C ADD CONSTANT TO DIAGONAL DARA1540
IF(WORK(NNE))17,17,19 DARA1550
15 IF(INCR)19,19,16 DARA1560
16 K=IPQ DARA1570
17 J=NNE-1 DARA1580
DO 18 I=1,K DARA1590
WORK(J)=WORK(J)+DIAG DARA1600
18 J=J+I DARA1610
C DARA1620
C SOLVE NORMAL EQUATIONS DARA1630
19 CALL DAPFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER) DARA1640
C DARA1650
C CHECK FOR FAILURE OF EQUATION SOLVER DARA1660
IF(IRES)4,4,20 DARA1670
C DARA1680
C TEST FOR DEFECTIVE NORMALEQUATIONS DARA1690
20 IF(IRES-IX)21,24,24 DARA1700
21 IF(INCR)22,22,23 DARA1710
22 DIAG=DIAG*0.125D0 DARA1720
23 DIAG=DIAG+DIAG DARA1730
INCR=INCR+1 DARA1740
C DARA1750
C START WITH OVER RELAXATION DARA1760
RELAX=8.D0 DARA1770
IF(INCR-LIMIT)12,45,45 DARA1780
C DARA1790
C CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR DARA1800
24 L=NDP DARA1810
J=NNE+IRES*(IRES-1)/2-1 DARA1820
K=J+IQ DARA1830
WORK(J)=0.D0 DARA1840
IRQ=IQ DARA1850
IRP=IRES-IQ+1 DARA1860
IF(IRP)25,26,26 DARA1870
25 IRQ=IRES+1 DARA1880
26 DO 29 I=1,N DARA1890
T=DATI(I) DARA1900
WORK(I)=0.D0 DARA1910
CALL DCNPS(WORK(I),T,WORK(K),IRP) DARA1920
M=L+N DARA1930
CALL DCNPS(WORK(M),T,WORK(J),IRQ) DARA1940
IF(WORK(M)*WORK(L))27,29,29 DARA1950
27 SUM=WORK(L)/WORK(M) DARA1960
IF(RELAX+SUM)29,29,28 DARA1970
28 RELAX=-SUM DARA1980
29 L=L+1 DARA1990
C DARA2000
C MODIFY RELAXATION FACTOR IF NECESSARY DARA2010
SSOE=OSUM DARA2020
ITER=LIMIT DARA2030
30 SUM=0.D0 DARA2040
RELAX=RELAX*0.5D0 DARA2050
DO 32 I=1,N DARA2060
M=I+N DARA2070
K=M+N DARA2080
L=K+N DARA2090
SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L)) DARA2100
SAVE=SAVE*SAVE DARA2110
IF(DATI(NDP))32,32,31 DARA2120
31 SAVE=SAVE*DATI(K) DARA2130
32 SUM=SUM+SAVE DARA2140
IF(ITER)45,33,33 DARA2150
33 ITER=ITER-1 DARA2160
IF(SUM-OSUM)34,37,35 DARA2170
34 OSUM=SUM DARA2180
GOTO 30 DARA2190
C DARA2200
C TEST FOR IMPROVEMENT DARA2210
35 IF(OSUM-SSOE)36,30,30 DARA2220
36 RELAX=RELAX+RELAX DARA2230
37 T=0. DARA2240
SAVE=0.D0 DARA2250
K=IRES+1 DARA2260
DO 38 I=2,K DARA2270
J=J+1 DARA2280
T=T+DABS(P(I)) DARA2290
P(I)=P(I)+RELAX*WORK(J) DARA2300
38 SAVE=SAVE+DABS(P(I)) DARA2310
C DARA2320
C UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR DARA2330
DO 39 I=1,N DARA2340
J=I+N DARA2350
K=J+N DARA2360
L=K+N DARA2370
WORK(J)=WORK(J)+RELAX*WORK(I) DARA2380
39 WORK(K)=WORK(K)+RELAX*WORK(L) DARA2390
C DARA2400
C TEST FOR CONVERGENCE DARA2410
IF(INCR)40,40,42 DARA2420
40 IF(SSOE-OSUM-RELAX*OSUM*DBLE(EPS))46,46,41 DARA2430
41 IF(DABS(T-SAVE)-RELAX*SAVE*DBLE(EPS))46,46,42 DARA2440
42 IF(OSUM-SAVE*DBLE(ETA))46,46,43 DARA2450
43 KOUNT=KOUNT+1 DARA2460
IF(KOUNT-LIMIT)10,44,44 DARA2470
C DARA2480
C ERROR RETURN IN CASE OF POOR CONVERGENCE DARA2490
44 IER=2 DARA2500
RETURN DARA2510
45 IER=1 DARA2520
RETURN DARA2530
C DARA2540
C NORMAL RETURN DARA2550
46 IER=0 DARA2560
RETURN DARA2570
END DARA2580