Directory of image this file is from
This file as a plain text file
C
C ..................................................................
C
C SUBROUTINE ARAT
C
C PURPOSE
C CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE
C FUNCTION IN THE LEAST SQUARES SENSE
C
C USAGE
C CALL ARAT(DATI,N,WORK,P,IP,IQ,IER)
C
C DESCRIPTION OF PARAMETERS
C DATI - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS
C THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS,
C THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND
C THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY.
C IF NO WEIGHTS ARE TO BE USED THEN THE THIRD
C COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT
C WHICH MUST CONTAIN A NONPOSITIVE VALUE
C N - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION
C WORK - WORKING STORAGE WHICH IS OF DIMENSION
C (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST.
C ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED
C IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF
C THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO
C WORK(3*N)
C P - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND
C NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ
C LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP
C LOCATIONS.
C COEFFICIENTS ARE ORDERED FROM LOW TO HIGH.
C IP - DIMENSION OF THE NUMERATOR (INPUT VALUE)
C IQ - DIMENSION OF THE DENOMINATOR (INPUT VALUE)
C IER - RESULTANT ERROR PARAMETER
C IER =-1 MEANS FORMAL ERRORS
C IER = 0 MEANS NO ERRORS
C IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION
C IER IS ALSO USED AS INPUT VALUE
C A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN
C INITIAL APPROXIMATION STORED IN P
C
C REMARKS
C THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR
C OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P
C STARTING WITH LOW POWERS (DENOMINATOR FIRST).
C IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE.
C SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL
C FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL
C (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEAR
C TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS.
C IF A FIT IN OTHER FUNCTIONS IS REQUIRED, CNP AND CNPS MUST
C BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C APLL, APFS, FRAT, CNPS, CNP
C CNP IS REQUIRED WITHIN FRAT
C
C METHOD
C THE ITERATIVE SCHEME USED FOR CALCULATION OF THE
C APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS
C WHICH ARE OBTAINED BY LINEARIZATION.
C A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH
C IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF
C ZEROES WITHIN THE APPROXIMATION INTERVAL.
C FOR REFERENCE SEE
C D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN,
C COMPUTING(1966), VOL.1, ED.3, PP.264-272.
C D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION
C OF NONLINEAR PARAMETERS,
C JSIAM(1963), VOL.11, ED.2, PP.431-441.
C
C ..................................................................
C
SUBROUTINE ARAT(DATI,N,WORK,P,IP,IQ,IER)
C
C
EXTERNAL FRAT
C
C DIMENSIONED LOCAL VARIABLE
DIMENSION IERV(3)
C
C DIMENSIONED DUMMY VARIABLES
DIMENSION DATI(1),WORK(1),P(1)
C
C INITIALIZE TESTVALUES
LIMIT=20
ETA =1.E-11
EPS=1.E-5
C
C CHECK FOR FORMAL ERRORS
IF(N)4,4,1
1 IF(IP)4,4,2
2 IF(IQ)4,4,3
3 IPQ=IP+IQ
IF(N-IPQ)4,5,5
C
C ERROR RETURN IN CASE OF FORMAL ERRORS
4 IER=-1
RETURN
C
C INITIALIZE ITERATION PROCESS
5 KOUNT=0
IERV(2)=IP
IERV(3)=IQ
NDP=N+N+1
NNE=NDP+NDP
IX=IPQ-1
IQP1=IQ+1
IRHS=NNE+IPQ*IX/2
IEND=IRHS+IX
C
C TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION
IF(IER)8,6,8
C
C INITIALIZE NUMERATOR AND DENOMINATOR
6 DO 7 I=2,IPQ
7 P(I)=0.
P(1)=1.
C
C CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL
C APPROXIMATION
8 DO 9 J=1,N
T=DATI(J)
I=J+N
CALL CNPS(WORK(I),T,P(IQP1),IP)
K=I+N
9 CALL CNPS(WORK(K),T,P,IQ)
C
C SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION)
10 CALL APLL(FRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV)
C
C CHECK FOR ZERO DENOMINATOR
IF(IERV(1))4,11,4
11 INCR=0
RELAX=2.
C
C RESTORE MATRIX IN WORKING STORAGE
12 J=IEND
DO 13 I=NNE,IEND
J=J+1
13 WORK(I)=WORK(J)
IF(KOUNT)14,14,15
C
C SAVE SQUARE SUM OF ERRORS
14 OSUM=WORK(IEND)
DIAG=OSUM*EPS
K=IQ
C
C ADD CONSTANT TO DIAGONAL
IF(WORK(NNE))17,17,19
15 IF(INCR)19,19,16
16 K=IPQ
17 J=NNE-1
DO 18 I=1,K
WORK(J)=WORK(J)+DIAG
18 J=J+I
C
C SOLVE NORMAL EQUATIONS
19 CALL APFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER)
C
C CHECK FOR FAILURE OF EQUATION SOLVER
IF(IRES)4,4,20
C
C TEST FOR DEFECTIVE NORMALEQUATIONS
20 IF(IRES-IX)21,24,24
21 IF(INCR)22,22,23
22 DIAG=DIAG*0.125
23 DIAG=DIAG+DIAG
INCR=INCR+1
C
C START WITH OVER RELAXATION
RELAX=8.
IF(INCR-LIMIT)12,45,45
C
C CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR
24 L=NDP
J=NNE+IRES*(IRES-1)/2-1
K=J+IQ
WORK(J)=0.
IRQ=IQ
IRP=IRES-IQ+1
IF(IRP)25,26,26
25 IRQ=IRES+1
26 DO 29 I=1,N
T=DATI(I)
WORK(I)=0.
CALL CNPS(WORK(I),T,WORK(K),IRP)
M=L+N
CALL CNPS(WORK(M),T,WORK(J),IRQ)
IF(WORK(M)*WORK(L))27,29,29
27 SUM=WORK(L)/WORK(M)
IF(RELAX+SUM)29,29,28
28 RELAX=-SUM
29 L=L+1
C
C MODIFY RELAXATION FACTOR IF NECESSARY
SSOE=OSUM
ITER=LIMIT
30 SUM=0.
RELAX=RELAX*0.5
DO 32 I=1,N
M=I+N
K=M+N
L=K+N
SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L))
SAVE=SAVE*SAVE
IF(DATI(NDP))32,32,31
31 SAVE=SAVE*DATI(K)
32 SUM=SUM+SAVE
IF(ITER)45,33,33
33 ITER=ITER-1
IF(SUM-OSUM)34,37,35
34 OSUM=SUM
GOTO 30
C
C TEST FOR IMPROVEMENT
35 IF(OSUM-SSOE)36,30,30
36 RELAX=RELAX+RELAX
37 T=0.
SAVE=0.
K=IRES+1
DO 38 I=2,K
J=J+1
T=T+ABS(P(I))
P(I)=P(I)+RELAX*WORK(J)
38 SAVE=SAVE+ABS(P(I))
C
C UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR
DO 39 I=1,N
J=I+N
K=J+N
L=K+N
WORK(J)=WORK(J)+RELAX*WORK(I)
39 WORK(K)=WORK(K)+RELAX*WORK(L)
C
C TEST FOR CONVERGENCE
IF(INCR)40,40,42
40 IF(SSOE-OSUM-RELAX*EPS*OSUM)46,46,41
41 IF(ABS(T-SAVE)-RELAX*EPS*SAVE)46,46,42
42 IF(OSUM-ETA*SAVE)46,46,43
43 KOUNT=KOUNT+1
IF(KOUNT-LIMIT)10,44,44
C
C ERROR RETURN IN CASE OF POOR CONVERGENCE
44 IER=2
RETURN
45 IER=1
RETURN
C
C NORMAL RETURN
46 IER=0
RETURN
END