Chris Marshall >
PDL-2.007_03 >
PDL::Slatec

PDL::Slatec - PDL interface to the slatec numerical programming library

use PDL::Slatec; ($ndeg, $r, $ierr, $a) = polyfit($x, $y, $w, $maxdeg, $eps);

This module serves the dual purpose of providing an interface to parts of the slatec library and showing how to interface PDL to an external library. Using this library requires a fortran compiler; the source for the routines is provided for convenience.

Currently available are routines to: manipulate matrices; calculate FFT's; fit data using polynomials; and interpolate/integrate data using piecewise cubic Hermite interpolation.

PCHIP is the slatec package of routines to perform piecewise cubic Hermite interpolation of data. It features software to produce a monotone and "visually pleasing" interpolant to monotone data. According to Fritsch & Carlson ("Monotone piecewise cubic interpolation", SIAM Journal on Numerical Analysis 17, 2 (April 1980), pp. 238-246), such an interpolant may be more reasonable than a cubic spline if the data contains both "steep" and "flat" sections. Interpolation of cumulative probability distribution functions is another application. These routines are cryptically named (blame FORTRAN), beginning with 'ch', and accept either float or double piddles.

Most of the routines require an integer parameter called `check`

; if set to 0, then no checks on the validity of the input data are made, otherwise these checks are made. The value of `check`

can be set to 0 if a routine such as chim has already been successfully called.

- If not known, estimate derivative values for the points using the chim, chic, or chsp routines (the following routines require both the function (
`f`

) and derivative (`d`

) values at a set of points (`x`

)). - Evaluate, integrate, or differentiate the resulting PCH function using the routines: chfd; chfe; chia; chid.
- If desired, you can check the monotonicity of your data using chcm.

Eigenvalues and eigenvectors of a real positive definite symmetric matrix.

($eigvals,$eigvecs) = eigsys($mat)

Note: this function should be extended to calculate only eigenvalues if called in scalar context!

Inverse of a square matrix

($inv) = matinv($mat)

Convenience wrapper routine about the `polfit`

`slatec`

function. Separates supplied arguments and return values.

Fit discrete data in a least squares sense by polynomials in one variable. Handles threading correctly--one can pass in a 2D PDL (as `$y`

) and it will pass back a 2D PDL, the rows of which are the polynomial regression results (in `$r`

corresponding to the rows of $y.

($ndeg, $r, $ierr, $a, $coeffs, $rms) = polyfit($x, $y, $w, $maxdeg, [$eps]); $coeffs = polyfit($x,$y,$w,$maxdeg,[$eps]);

where on input:

`$x`

and `$y`

are the values to fit to a polynomial. `$w`

are weighting factors `$maxdeg`

is the maximum degree of polynomial to use and `$eps`

is the required degree of fit.

and the output switches on list/scalar context.

In list context:

`$ndeg`

is the degree of polynomial actually used `$r`

is the values of the fitted polynomial `$ierr`

is a return status code, and `$a`

is some working array or other (preserved for historical purposes) `$coeffs`

is the polynomial coefficients of the best fit polynomial. `$rms`

is the rms error of the fit.

In scalar context, only $coeffs is returned.

Historically, `$eps`

was modified in-place to be a return value of the rms error. This usage is deprecated, and `$eps`

is an optional parameter now. It is still modified if present.

`$a`

is a working array accessible to Slatec - you can feed it to several other Slatec routines to get nice things out. It does not thread correctly and should probably be fixed by someone. If you are reading this, that someone might be you.

This version of polyfit handles bad values correctly. Bad values in $y are ignored for the fit and give computed values on the fitted curve in the return. Bad values in $x or $w are ignored for the fit and result in bad elements in the output.

Convenience wrapper routine around the `pcoef`

`slatec`

function. Separates supplied arguments and return values.

Convert the `polyfit`

/`polfit`

coefficients to Taylor series form.

$tc = polycoef($l, $c, $a);

Convenience wrapper routine around the `pvalue`

`slatec`

function. Separates supplied arguments and return values.

For multiple input x positions, a corresponding y position is calculated.

The derivatives PDL is one dimensional (of size `nder`

) if a single x position is supplied, two dimensional if more than one x position is supplied.

Use the coefficients generated by `polyfit`

(or `polfit`

) to evaluate the polynomial fit of degree `l`

, along with the first `nder`

of its derivatives, at a specified point.

($yfit, $yp) = polyvalue($l, $nder, $x, $a);

compute the determinant of an invertible matrix

$mat = zeroes(5,5); $mat->diagonal(0,1) .= 1; # unity matrix $det = detslatec $mat;

Usage:

$determinant = detslatec $matrix;

Signature: detslatec(mat(n,m); [o] det())

`detslatec`

computes the determinant of an invertible matrix and barfs if the matrix argument provided is non-invertible. The matrix threads as usual.

This routine was previously known as `det`

which clashes now with det which is provided by PDL::MatrixOps. For the moment PDL::Slatec will also load PDL::MatrixOps thereby making sure that older scripts work.

Fast Fourier Transform

$v_in = pdl(1,0,1,0); ($azero,$a,$b) = PDL::Slatec::fft($v_in);

`PDL::Slatec::fft`

is a convenience wrapper for ezfftf, and performs a Fast Fourier Transform on an input vector `$v_in`

. The return values are the same as for ezfftf.

reverse Fast Fourier Transform

$v_out = PDL::Slatec::rfft($azero,$a,$b); print $v_in, $vout [1 0 1 0] [1 0 1 0]

`PDL::Slatec::rfft`

is a convenience wrapper for ezfftb, and performs a reverse Fast Fourier Transform. The input is the same as the output of PDL::Slatec::fft, and the output of `rfft`

is a data vector, similar to what is input into PDL::Slatec::fft.

- 0 if successful.
- > 0 if there were
`ierr`

switches in the direction of monotonicity (data still valid). - -1 if
`nelem($x) < 2`

. - -3 if
`$x`

is not strictly increasing.

- 1 if first derivative at
`x(0)`

is given in`vc(0)`

. - 2 if second derivative at
`x(0)`

is given in`vc(0)`

. - 3 to use the 3-point difference formula for
`d(0)`

. (Reverts to the default b.c. if`n < 3`

) - 4 to use the 4-point difference formula for
`d(0)`

. (Reverts to the default b.c. if`n < 4`

) - 5 to set
`d(0)`

so that the second derivative is continuous at`x(1)`

. (Reverts to the default b.c. if`n < 4`

)

The values for `ic(1)`

are the same as above, except that the first-derivative value is stored in `vc(1)`

for cases 1 and 2. The values of `$vc`

need only be set if options 1 or 2 are chosen for `$ic`

.

Set `$mflag = 0`

if interpolant is required to be monotonic in each interval, regardless of the data. This causes `$d`

to be set to 0 at all switch points. Set `$mflag`

to be non-zero to use a formula based on the 3-point difference formula at switch points. If `$mflag > 0`

, then the interpolant at swich points is forced to not deviate from the data by more than `$mflag*dfloc`

, where `dfloc`

is the maximum of the change of `$f`

on this interval and its two immediate neighbours. If `$mflag < 0`

, no such control is to be imposed.

The piddle `$wk`

is only needed for work space. However, I could not get it to work as a temporary variable, so you must supply it; it is a 1D piddle with `2*n`

elements.

Error status returned by `$ierr`

:

- 0 if successful.
- 1 if
`ic(0) < 0`

and`d(0)`

had to be adjusted for monotonicity. - 2 if
`ic(1) < 0`

and`d(n-1)`

had to be adjusted for monotonicity. - 3 if both 1 and 2 are true.
- -1 if
`n < 2`

. - -3 if
`$x`

is not strictly increasing. - -4 if
`abs(ic(0)) > 5`

. - -5 if
`abs(ic(1)) > 5`

. - -6 if both -4 and -5 are true.
- -7 if
`nwk < 2*(n-1)`

.

- 0 to set
`d(0)`

so that the third derivative is continuous at`x(1)`

. - 1 if first derivative at
`x(0)`

is given in`vc(0`

). - 2 if second derivative at
`x(0`

) is given in`vc(0)`

. - 3 to use the 3-point difference formula for
`d(0)`

. (Reverts to the default b.c. if`n < 3`

.) - 4 to use the 4-point difference formula for
`d(0)`

. (Reverts to the default b.c. if`n < 4`

.)

The values for `ic(1)`

are the same as above, except that the first-derivative value is stored in `vc(1)`

for cases 1 and 2. The values of `$vc`

need only be set if options 1 or 2 are chosen for `$ic`

.

The piddle `$wk`

is only needed for work space. However, I could not get it to work as a temporary variable, so you must supply it; it is a 1D piddle with `2*n`

elements.

Error status returned by `$ierr`

:

- 0 if successful.
- -1 if
`nelem($x) < 2`

. - -3 if
`$x`

is not strictly increasing. - -4 if
`ic(0) < 0`

or`ic(0) > 4`

. - -5 if
`ic(1) < 0`

or`ic(1) > 4`

. - -6 if both of the above are true.
- -7 if
`nwk < 2*n`

. - -8 in case of trouble solving the linear system for the interior derivative values.

- 0 if successful.
- >0 if extrapolation was performed at
`ierr`

points (data still valid). - -1 if
`nelem($x) < 2`

- -3 if
`$x`

is not strictly increasing. - -4 if
`nelem($xe) < 1`

. - -5 if an error has occurred in a lower-level routine, which should never happen.

- 0 if successful.
- >0 if extrapolation was performed at
`ierr`

points (data still valid). - -1 if
`nelem($x) < 2`

- -3 if
`$x`

is not strictly increasing. - -4 if
`nelem($xe) < 1`

.

- 0 if successful.
- 1 if
`$a`

lies outside`$x`

. - 2 if
`$b`

lies outside`$x`

. - 3 if both 1 and 2 are true.
- -1 if
`nelem($x) < 2`

- -3 if
`$x`

is not strictly increasing. - -4 if an error has occurred in a lower-level routine, which should never happen.

- 0 if successful.
- -1 if
`nelem($x) < 2`

. - -3 if
`$x`

is not strictly increasing. - -4 if
`$ia`

or`$ib`

is out of range.

- -3 if function is probably decreasing
- -1 if function is strictly decreasing
- 0 if function is constant
- 1 if function is strictly increasing
- 2 if function is non-monotonic
- 3 if function is probably increasing

If `abs(ismon(i)) == 3`

, the derivative values are near the boundary of the monotonicity region. A small increase produces non-monotonicity, whereas a decrease produces strict monotonicity.

The above applies to `i = 0 .. nelem($x)-1`

. The last element of `$ismon`

indicates whether the entire function is monotonic over $x.

Error status returned by `$ierr`

:

- 0 if successful.
- -1 if
`n < 2`

. - -3 if
`$x`

is not strictly increasing.

ignore function for now although code is in slatec/ directory

- 0 - Quadruple knots at the first and last points.
- 1 - Replicate lengths of extreme subintervals:
`t( 0 ) = t( 1 ) = x(0) - (x(1)-x(0))`

and`t(m+3) = t(m+2) = x(n-1) + (x(n-1)-x(n-2))`

- 2 - Periodic placement of boundary knots:
`t( 0 ) = t( 1 ) = x(0) - (x(n-1)-x(n-2))`

and`t(m+3) = t(m+2) = x(n) + (x(1)-x(0))`

- <0 - Assume the
`nknots`

and`t`

were set in a previous call.

`nknots`

is the number of knots and may be changed by the routine. If `knotyp >= 0`

, `nknots`

will be set to `ndim+4`

, otherwise it is an input variable, and an error will occur if its value is not equal to `ndim+4`

.

`t`

is the array of `2*n+4`

knots for the B-representation and may be changed by the routine. If `knotyp >= 0`

, `t`

will be changed so that the interior double knots are equal to the x-values and the boundary knots set as indicated above, otherwise it is assumed that `t`

was set by a previous call (no check is made to verify that the data forms a legitimate knot sequence).

Error status returned by `$ierr`

:

- 0 if successful.
- -4 if
`knotyp > 2`

. - -5 if
`knotyp < 0`

and`nknots != 2*n + 4`

.

Copyright (C) 1997 Tuomas J. Lukka. Copyright (C) 2000 Tim Jenness, Doug Burke. All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file.

syntax highlighting: