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 CWcheck; if set to 0, then no checks on the validity of the input data are made, otherwise these checks are made. The value of CWcheck 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 (CWf) and derivative (CWd) values at a set of points (CWx)).

•

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 CWpolfit CWslatec function. Separates supplied arguments and return values.

Fit discrete data in a least squares sense by polynomials in one variable. Handles threading correctlyone can pass in a 2D PDL (as CW$y) and it will pass back a 2D PDL, the rows of which are the polynomial regression results (in CW$r corresponding to the rows of CW$y.

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

 where on input:

 C<x> and C<y> are the values to fit to a polynomial.
 C<w> are weighting factors
 C<maxdeg> is the maximum degree of polynomial to use and 
 C<eps> is the required degree of fit.

 and on output:

 C<ndeg> is the degree of polynomial actually used
 C<r> is the values of the fitted polynomial 
 C<ierr> is a return status code, and
 C<a> is some working array or other
 C<eps> is modified to contain the rms error of the fit.

This version of polyfit handles bad values correctly. It strips them out of the CW$x variable and creates an appropriate CW$y variable containing indices of the non-bad members of CW$x before calling the Slatec routine CWpolfit.

Convenience wrapper routine around the CWpcoef CWslatec function. Separates supplied arguments and return values.

Convert the CWpolyfit/CWpolfit coefficients to Taylor series form.

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

Convenience wrapper routine around the CWpvalue CWslatec 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 CWnder) if a single x position is supplied, two dimensional if more than one x position is supplied.

Use the coefficients generated by CWpolyfit (or CWpolfit) to evaluate the polynomial fit of degree CWl, along with the first CWnder 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())

CWdetslatec 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 CWdet which clashes now with <det:PDL::MatrixOps/det> which is provided by PDL::MatrixOps. For the moment PDL::Slatec will also load PDL::MatrixOps thereby making sure that older scripts work.

  Signature: (x(n,p);[o]s(p);[o]e(p);[o]u(n,p);[o]v(p,p);[o]work(n);int job();int [o]info())

singular value decomposition of a matrix

svdc does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(n,n);rcond();[o]z(n);int [o]info())

Factor a real symmetric positive definite matrix and estimate the condition number of the matrix.

poco does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(n,n);int [o]ipvt(n);[o]rcond();[o]z(n))

Factor a matrix using Gaussian elimination and estimate the condition number of the matrix.

geco does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(n,n);int [o]ipvt(n);int [o]info())

Factor a matrix using Gaussian elimination.

gefa does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(n,n);[o]det(two=2);int job())

Compute the determinant and inverse of a certain real symmetric positive definite matrix using the factors computed by poco.

podi does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(n,n);int [o]ipvt(n);[o]det(two=2);[o]work(n);int job())

Compute the determinant and inverse of a matrix using the factors computed by geco or gefa.

gedi does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(lda,n);int ipvt(n);b(n);int job())

Solve the real system CWA*X=B or CWTRANS(A)*X=B using the factors computed by geco or gefa.

gesl does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (a(n,n);[o]w(n);int matz();[o]z(n,n);[t]fvone(n);[t]fvtwo(n);int [o]ierr())

This subroutine calls the recommended sequence of subroutines from the eigensystem subroutine package (EISPACK) to find the eigenvalues and eigenvectors (if desired) of a REAL SYMMETRIC matrix.

rs does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (int n();[o]wsave(foo))

Subroutine ezffti initializes the work array CWwsave() which is used in both ezfftf and ezfftb. The prime factorization of CWn together with a tabulation of the trigonometric functions are computed and stored in CWwsave().

ezffti does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo))

ezfftf does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: ([o]r(n);azero();a(n);b(n);wsave(foo))

ezfftb does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (int l();c();[o]tc(bar);a(foo))

Convert the CWpolfit coefficients to Taylor series form. CWc and CWa() must be of the same type.

pcoef does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (int l();x();[o]yfit();[o]yp(n);a(foo))

Use the coefficients generated by CWpolfit to evaluate the polynomial fit of degree CWl, along with the first CWnder of its derivatives, at a specified point. CWx and CWa must be of the same type.

pvalue does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);f(n);[o]d(n);int [o]ierr())

Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.

Calculate the derivatives at the given set of points (CW$x,$f, where CW$x is strictly increasing). The resulting set of points - CW$x,$f,$d, referred to as the cubic Hermite representation - can then be used in other functions, such as chfe, chfd, and chia.

The boundary conditions are compatible with monotonicity, and if the data are only piecewise monotonic, the interpolant will have an extremum at the switch points; for more control over these issues use chic.

Error status returned by CW$ierr:

•

0 if successful.

•

> 0 if there were CWierr switches in the direction of monotonicity (data still valid).

•

-1 if CWnelem($x) < 2.

•

-3 if CW$x is not strictly increasing.

chim does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (int ic(two=2);vc(two=2);mflag();x(n);f(n);[o]d(n);wk(n);int [o]ierr())

Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.

Calculate the derivatives at the given points (CW$x,$f, where CW$x is strictly increasing). Control over the boundary conditions is given by the CW$ic and CW$vc piddles, and the value of CW$mflag determines the treatment of points where monotoncity switches direction. A simpler, more restricted, interface is available using chim.

The first and second elements of CW$ic determine the boundary conditions at the start and end of the data respectively. If the value is 0, then the default condition, as used by chim, is adopted. If greater than zero, no adjustment for monotonicity is made, otherwise if less than zero the derivative will be adjusted. The allowed magnitudes for CWic(0) are:

•

1 if first derivative at CWx(0) is given in CWvc(0).

•

2 if second derivative at CWx(0) is given in CWvc(0).

•

3 to use the 3-point difference formula for CWd(0). (Reverts to the default b.c. if CWn < 3)

•

4 to use the 4-point difference formula for CWd(0). (Reverts to the default b.c. if CWn < 4)

•

5 to set CWd(0) so that the second derivative is continuous at CWx(1). (Reverts to the default b.c. if CWn < 4)

The values for CWic(1) are the same as above, except that the first-derivative value is stored in CWvc(1) for cases 1 and 2. The values of CW$vc need only be set if options 1 or 2 are chosen for CW$ic.

Set CW$mflag = 0 if interpolant is required to be monotonic in each interval, regardless of the data. This causes CW$d to be set to 0 at all switch points. Set CW$mflag to be non-zero to use a formula based on the 3-point difference formula at switch points. If CW$mflag > 0, then the interpolant at swich points is forced to not deviate from the data by more than CW$mflag*dfloc, where CWdfloc is the maximum of the change of CW$f on this interval and its two immediate neighbours. If CW$mflag < 0, no such control is to be imposed.

The piddle CW$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 CW2*n elements.

Error status returned by CW$ierr:

•

0 if successful.

•

1 if CWic(0) < 0 and CWd(0) had to be adjusted for monotonicity.

•

2 if CWic(1) < 0 and CWd(n-1) had to be adjusted for monotonicity.

•

3 if both 1 and 2 are true.

•

-1 if CWn < 2.

•

-3 if CW$x is not strictly increasing.

•

-4 if CWabs(ic(0)) > 5.

•

-5 if CWabs(ic(1)) > 5.

•

-6 if both -4 and -5 are true.

•

-7 if CWnwk < 2*(n-1).

chic does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (int ic(two=2);vc(two=2);x(n);f(n);[o]d(n);wk(n);int [o]ierr())

Calculate the derivatives of (x,f(x)) using cubic spline interpolation.

Calculate the derivatives, using cubic spline interpolation, at the given points (CW$x,$f), with the specified boundary conditions. Control over the boundary conditions is given by the CW$ic and CW$vc piddles. The resulting values - CW$x,$f,$d - can be used in all the functions which expect a cubic Hermite function.

The first and second elements of CW$ic determine the boundary conditions at the start and end of the data respectively. The allowed values for CWic(0) are:

•

0 to set CWd(0) so that the third derivative is continuous at CWx(1).

•

1 if first derivative at CWx(0) is given in CWvc(0).

•

2 if second derivative at CWx(0) is given in CWvc(0).

•

3 to use the 3-point difference formula for CWd(0). (Reverts to the default b.c. if CWn < 3.)

•

4 to use the 4-point difference formula for CWd(0). (Reverts to the default b.c. if CWn < 4.)

The values for CWic(1) are the same as above, except that the first-derivative value is stored in CWvc(1) for cases 1 and 2. The values of CW$vc need only be set if options 1 or 2 are chosen for CW$ic.

The piddle CW$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 CW2*n elements.

Error status returned by CW$ierr:

•

0 if successful.

•

-1 if CWnelem($x) < 2.

•

-3 if CW$x is not strictly increasing.

•

-4 if CWic(0) < 0 or CWic(0) > 4.

•

-5 if CWic(1) < 0 or CWic(1) > 4.

•

-6 if both of the above are true.

•

-7 if CWnwk < 2*n.

•

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

chsp does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);f(n);d(n);int check();xe(n);[o]fe(n);[o]de(n);int [o]ierr())

Interpolate function and derivative values.

Given a piecewise cubic Hermite function - such as from chim - evaluate the function (CW$fe) and derivative (CW$de) at a set of points (CW$xe). If function values alone are required, use chfe. Set CWcheck to 0 to skip checks on the input data.

Error status returned by CW$ierr:

•

0 if successful.

•

>0 if extrapolation was performed at CWierr points (data still valid).

•

-1 if CWnelem($x) < 2

•

-3 if CW$x is not strictly increasing.

•

-4 if CWnelem($xe) < 1.

•

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

chfd does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);f(n);d(n);int check();xe(n);[o]fe(n);int [o]ierr())

Interpolate function values.

Given a piecewise cubic Hermite function - such as from chim - evaluate the function (CW$fe) at a set of points (CW$xe). If derivative values are also required, use chfd. Set CWcheck to 0 to skip checks on the input data.

Error status returned by CW$ierr:

•

0 if successful.

•

>0 if extrapolation was performed at CWierr points (data still valid).

•

-1 if CWnelem($x) < 2

•

-3 if CW$x is not strictly increasing.

•

-4 if CWnelem($xe) < 1.

chfe does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);f(n);d(n);int check();a();b();[o]ans();int [o]ierr())

Integrate (x,f(x)) over arbitrary limits.

Evaluate the definite integral of a a piecewise cubic Hermite function over an arbitrary interval, given by CW[$a,$b]. See chid if the integration limits are data points. Set CWcheck to 0 to skip checks on the input data.

The values of CW$a and CW$b do not have to lie within CW$x, although the resulting integral value will be highly suspect if they are not.

Error status returned by CW$ierr:

•

0 if successful.

•

1 if CW$a lies outside CW$x.

•

2 if CW$b lies outside CW$x.

•

3 if both 1 and 2 are true.

•

-1 if CWnelem($x) < 2

•

-3 if CW$x is not strictly increasing.

•

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

chia does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);f(n);d(n);int check();int ia();int ib();[o]ans();int [o]ierr())

Integrate (x,f(x)) between data points.

Evaluate the definite integral of a a piecewise cubic Hermite function between CWx($ia) and CWx($ib).

See chia for integration between arbitrary limits.

Although using a fortran routine, the values of CW$ia and CW$ib are zero offset. Set CWcheck to 0 to skip checks on the input data.

Error status returned by CW$ierr:

•

0 if successful.

•

-1 if CWnelem($x) < 2.

•

-3 if CW$x is not strictly increasing.

•

-4 if CW$ia or CW$ib is out of range.

chid does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);f(n);d(n);int check();int [o]ismon(n);int [o]ierr())

Check the given piecewise cubic Hermite function for monotonicity.

The outout piddle CW$ismon indicates over which intervals the function is monotonic. Set CWcheck to 0 to skip checks on the input data.

For the data interval CW[x(i),x(i+1)], the values of CWismon(i) can be:

•

-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 CWabs(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 CWi = 0 .. nelem($x)-1. The last element of CW$ismon indicates whether the entire function is monotonic over CW$x.

Error status returned by CW$ierr:

•

0 if successful.

•

-1 if CWn < 2.

•

-3 if CW$x is not strictly increasing.

chcm does not process bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

  Signature: (x(n);y(n);w(n);int maxdeg();int [o]ndeg();[o]eps();[o]r(n);int [o]ierr();[o]a(foo);[t]xtmp(n);[t]ytmp(n);[t]wtmp(n);[t]rtmp(n))

Fit discrete data in a least squares sense by polynomials in one variable. CWx(), CWy() and CWw() must be of the same type. This version handles bad values appropriately

polfit does handle bad values. It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.

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.