Replacing calls to E02ADF

Superficially, the equivalent SLATEC routine is PDA_DPOLFT. Since NAG has no routine to do a polynomial extrapolation, Figaro usurped this routine with a peculiar set of weights and very few degrees of freedom to extrapolate a polynomial. Although PDA_DPOLFT can be used even in that case, it is PDA_DPLINT that is intended for polynomial interpolation.

The NAG code would look like

      INTEGER I, M, K, NROWS, IFAIL
      DOUBLE PRECISION X(M), Y(M), W(M)
      DOUBLE PRECISION WORK1(3*M), WORK2( 2*(K+1) )
      DOUBLE PRECISION A1(NROWS,K+1), S(K+1)
      DOUBLE PRECISION A2(K+1)
      IFAIL = 1
      CALL E02ADF( M, K+1, NROWS, X, Y, W, WORK1, WORK2, A1, S, IFAIL )
      IF ( IFAIL .NE. 0 ) THEN
         An error has occurred
      END IF
      DO 1 I = 1, K+1
         A2(I) = A1(K+1,I)
    1 CONTINUE

Here W are the weights proportional to the reciprocal of the standard deviation of Y. A1 returns a matrix of Chebyshev coefficients, one set of coefficients for each polynomial degree from 0 to K. The DO loop extracts the coefficients for degree K into A2. Note that these coefficients form a column rather than a row in A1. S returns the r.m.s. for the fit of each degree from 0 to K.

The behaviour of PDA_DPOLFT is controlled by the given value of EPS, passing zero (0D0) makes it perform fits for all degrees from 0 to K. EPS is also a returned argument, it returns the r.m.s. for the highest degree fitted. What degree that was is returned in NDEG. An indication of the success is returned in IFAIL1.

The weights should be proportional to the reciprocal of the variance, i.e. the square of the weights used for NAG.

The returned description of the polynomials A3 is rather different from the Chebyshev coefficients returned by E02ADF. A3 would be passed on to PDA_DP1VLU to evaluate the polynomial or to PDA_DPCOEF to convert A3 to coefficients of a Taylor series. A3 contains sufficient information to evaluate the polynomial of any degree from 0 to K. The desired degree is specified to PDA_DP1VLU or PDA_DPCOEF.

R is a returned vector in which the fit of highest degree is evaluated at all given X. Often this may render any further evaluation calls obsolete.

      INTEGER I, M, K, NDEG, IFAIL1, IFAIL2
      DOUBLE PRECISION X(M), Y(M), W(M), W2(M), R(M)
      DOUBLE PRECISION A3( 3*M + 3*(K+1) )
      DOUBLE PRECISION EPS, S(K+1)
      DO 1 I = 1, M
         W2(I) = W(I) * W(I)
    1 CONTINUE
      IFAIL2 = 0
      EPS = 0D0
      CALL PDA_DPOLFT( M, X, Y, W2, K, NDEG, EPS, R, IFAIL1, A3, IFAIL2 )
      IF ( NDEG .NE. K .OR. IFAIL1 .NE. 1 .OR. IFAIL2 .NE. 0 ) THEN
         An error has occurred
      END IF
      DO 2 I = 1, K
         S(I) = 0D0
    2 CONTINUE
      S(K+1) = EPS