The routines for this sort of application (and their origins) are:
E01BAF finds the interpolating cubic spline interpolant f(x) for a set of points (x,y). The interpolant is evaluated with E02BBF, evaluated with derivatives by E02BCF, and integrated by E02BDF. PDA_DBINTK does this, the order of the splines can be changed as well. This routine needs to be given the knots, while E01BAF set them itself. Evaluation of the interpolant and its derivatives is done by PDA_DBVALU, integration by PDA_DBSQAD.
E02BAF finds the fitting cubic spline f(x) for a set of points and weights (x,y,w). The interior knots 5 ... n+3 must be given and are fixed. The function is evaluated with E02BBF, with derivatives by E02BCF, and integrated with E02BDF. PDA_DEFC does this. All knots must be given, not just the interior ones. Instead of weights PDA_DEFC takes standard deviations (x,y,sigma) and uses 1/sigma as weight. Evaluation of the function and its derivatives is done by PDA_DBVALU, integration by PDA_DBSQAD.
E02BBF and E02BCF evaluate an interpolating or fitting cubic spline and its derivatives. They follow a call to E01BAF or E02BAF. This function is taken over by PDA_DBVALU. E02BDF integrates an interpolating or fitting cubic spline. It follows a call to E01BAF or E02BAF. This function is taken over by PDA_DBSQAD.
E02BEF finds the fitting cubic spline f(x) for a set of points and weights (x,y,w). The knots are located automatically. The function is evaluated with E02BBF, with derivatives by E02BCF, and integrated with E02BDF. PDA_CURFIT (from the DIERCKX package) solves this problem. While the other routines are from SLATEC and for double precision, PDA_CURFIT is for single precision. Hence, PDA_SPLEV, PDA_SPLDER and PDA_SPLINT should be used to evaluate the spline, its n-th derivative, and its integral.
E02ADF finds the fitting Chebyshev series minimising r.m.s. The Chebyshev series is equivalent to an ordinary polynomial, but cannot be extrapolated. The polynomial is evaluated by E02AEF or E02AKF. In the latter routine the Chebyshev coefficients can be one column of an array. It will also take the real-world x argument instead of the normalised x-bar argument within the range 1 ... +1.
In this library, PDA_DPOLFT fits an ordinary polynomial as a sum of orthogonal polynomials. The representation returned is somewhat special. It can be converted to coefficients of a Taylor series with PDA_DPCOEF or directly evaluated with PDA_DP1VLU. PDA_DP1VLU will return in one call as many derivatives as requested.
The SLATEC equivalent of this routine is PDA_DBINTK with order K = 4. The NAG code would look like
While the size of the coefficient vector can be reduced to M for PDA_DBINTK, you now need two work spaces. The major difference is that PDA_DBINTK needs to be given the knots. So you have to calculate them in the same way as E01BAF would have done. The handling of the status is different.
PDA_DEFC has not yet been used anywhere to replace E02BAF. Thus the migration hints given here may contain errors or may be based on misunderstandings.
The SLATEC equivalent of this routine is PDA_DEFC with order K = 4. The NAG code would look like
Here N+7 is the number of knots. Since the order is 4 (cubic), the number of interior knots is then N+78 = N1. N is the number of intervals. The interior knots are T(5) ... T(N+3). The weights W are reciprocal errors of Y. Although C has a length of N+7 only N+3 coefficients are returned by E02BAF.
The interior knots T(5) ... T(N+3) are given arguments, but the remaining knots are set by E02BAF and thus returned arguments.
In PDA_DEFC the order is K = 4. There are N+K+3 knots T(). T(1) ... T(K1) and T(N+5) ... T(N+K+3) are end knots. The next inner knots T(K) and T(N+4) are presumably the first and last x value. Then T(K+1) ... T(N+3) would be truly interior knots just as in the NAG code.
PDA_DEFC does not generate knots by itself. Contrary to the NAG code above, the first K and last K knots must be calculated before the call.
The size of the work space is more complex to calculate. PDA_DEFC needs the standard deviation in Y instead of the weights SD = 1/W. PDA_DEFC returns a diagnostic J, which should have value J = 1 if no error occurred.
The equivalent of this routine in SLATEC is PDA_DBVALU with the requested derivative being zero and the order being K = 4 for a cubic spline. The NAG code would look like
Here N, T and C are the same as in E02BAF. If T and C originate from a call to E01BAF then for N+7 read M+4 with M the number of data points given to the interpolation.
PDA_DBVALU is a function rather than a subroutine. The dimension passed to PDA_DBVALU is not that of T, but that of C, i.e. N+3 (or M after interpolation). PDA_DBVALU returns the value of any derivative, the fifth argument is zero so that it returns the function value itself. INVB must be given 1 in the first call. For several evaluations of the same spline it should not be changed between calls. It is changed by PDA_DBVALU. So if PDA_DBVALU is called in a DO loop, the statement INVB = 1 is typically before and outside the loop.
E02BCF has not yet been replaced anywhere. Thus the migration hints given here may contain errors or may be based on misunderstandings.
The equivalent of this routine in SLATEC is PDA_DBVALU. Several calls are necessary, one for each derivative. For the function value itself set the number of derivative to zero. The order is K = 4 for a cubic spline. The NAG code would look like
Here N, T and C are the same as in E02BAF. If T and C originate from a call to E01BAF then for N+7 read M+4 with M the number of data points given to the interpolation. S(I) returns the I-th derivative. In the case that X coincides with a knot and the derivatives are not continuous at that knot, LEFT is used to decide which side of the knot to use.
PDA_DBVALU is a function rather than a subroutine. The dimension passed to PDA_DBVALU is not that of T, but that of C, i.e. N+3 (or M after interpolation). PDA_DBVALU returns the value of any derivative, as specified in the fifth argument.
There is no equivalent to the LEFT parameter in NAG. PDA_DBVALU returns right limiting values, except at the right end point.
INVB must be given 1 in the first call. For several evaluations of the same spline it should not be changed between calls. It is changed by PDA_DBVALU. So if PDA_DBVALU is called in a DO loop, the statement INVB = 1 is typically before and outside the loop. In the code below, IFAIL is reset inside the DO loop. Assuming that an error will quit the loop, the IFAIL = 0 statement could be before and outside the DO loop as well.
PDA_DBSQAD has not yet been used anywhere to replace E02BDF. Thus the migration hints given here may contain errors or may be based on misunderstandings.
The equivalent of this routine in SLATEC is PDA_DBSQAD. The NAG code would look like
Here N, T and C are the same as in E02BAF. If T and C originate from a call to E01BAF then for N+7 read M+4 with M the number of data points given to the interpolation. DEFINT returns the integral over the whole x range where the spline is defined. This is from T(4) to T(N+4), which are most probably the smallest and largest X used in the fit or interpolation.
The dimension passed to PDA_DBSQAD is not that of T, but that of C, i.e. N+3 (or M after interpolation). PDA_DBSQAD calculates the integral for any interval on which the spline is defined. For the same interval as in the NAG code, the two limiting knots are given to PDA_DBSQAD.
E02BEF is a more advanced routine than the other NAG spline routines. It places knots automatically while fitting a cubic spline. This library includes for this case a number of routines from DIERCKX. PDA_CURFIT performs the spline approximation with given or automatic knots. The fitted function is evaluated with PDA_SPLEV, its derivatives with PDA_SPLDER, its integral with PDA_SPLINT. The routines exist only for single precision arguments.
These routines are so far unused, so there are no migration hints.
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
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.
The equivalent of this routine in SLATEC is PDA_DP1VLU. The NAG code would look like
Here A1 is the coefficient matrix returned by E02ADF, A2 is the column extracted for the required degree. XX is the x value for which an evaluation is required, it must be scaled into the range 1 ... +1, using the original extreme x values passed to E02ADF. P returns the function value.
PDA_DP1VLU returns any number of derivatives in addition to the function value. The second argument specifies how many derivatives are required. No scaling of XX is necessary, and no processing of the coefficients A3 as returned by PDA_DPOLFT.
E02AKF is similar to E02AEF but provides a different interface. Namely, A1 (from E02ADF) and the row length can be given instead of the extraction A2. Also XX is given rather than the scaled XCAP. There are no detailed migration hints for this routine yet.
This is not a NAG routine, but a routine in Figaro’s GEN library. It converts the vector of Chebyshev coefficients into a vector of coefficients of an ordinary polynomial. This is particularly useful when these are to be written to output files. In future the routine will still be needed to read old files that contain Chebyshev coefficients.
This routine is mentioned here really, because it has an equivalent in SLATEC named PDA_DPCOEF. PDA_DPCOEF is more general, in that it returns coefficients for a Taylor series, i.e. a polynomial in (XX-X0) for given expansion point X0.
The old code would look like
This would be replaced by
Just in case you are tempted to evaluate the polynomial yourself rather than use PDA_DP1VLU, here is a piece of code to get P = f(XX):