C
C     ..................................................................
C
C        SUBROUTINE DARAT
C
C        PURPOSE
C           CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE
C           FUNCTION IN THE LEAST SQUARES SENSE
C
C        USAGE
C           CALL DARAT(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                   DATI MUST BE OF DOUBLE PRECISION
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                   WORK MUST BE OF DOUBLE PRECISION
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                   P MUST BE OF DOUBLE PRECISION
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, DCNP AND DCNPS MUST
C           BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           DAPLL, DAPFS, DFRAT, DCNPS, DCNP
C           DCNP IS REQUIRED WITHIN DFRAT
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 DARAT(DATI,N,WORK,P,IP,IQ,IER)
C
C
      EXTERNAL DFRAT
C
C        DIMENSIONED LOCAL VARIABLE
      DIMENSION IERV(3)
C
C        DIMENSIONED DUMMY VARIABLES
      DIMENSION DATI(1),WORK(1),P(1)
      DOUBLE PRECISION DATI,WORK,P,T,OSUM,DIAG,RELAX,SUM,SSOE,SAVE
C
C        INITIALIZE TESTVALUES
      LIMIT=20
      ETA=1.E-29
      EPS=1.E-14
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.D0
      P(1)=1.D0
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 DCNPS(WORK(I),T,P(IQP1),IP)
      K=I+N
    9 CALL DCNPS(WORK(K),T,P,IQ)
C
C        SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION)
   10 CALL DAPLL(DFRAT,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.D0
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 DAPFS(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.125D0
   23 DIAG=DIAG+DIAG
      INCR=INCR+1
C
C        START WITH OVER RELAXATION
      RELAX=8.D0
      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.D0
      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.D0
      CALL DCNPS(WORK(I),T,WORK(K),IRP)
      M=L+N
      CALL DCNPS(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.D0
      RELAX=RELAX*0.5D0
      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.D0
      K=IRES+1
      DO 38 I=2,K
      J=J+1
      T=T+DABS(P(I))
      P(I)=P(I)+RELAX*WORK(J)
   38 SAVE=SAVE+DABS(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*OSUM*DBLE(EPS))46,46,41
   41 IF(DABS(T-SAVE)-RELAX*SAVE*DBLE(EPS))46,46,42
   42 IF(OSUM-SAVE*DBLE(ETA))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
