chic
 NAME
 SYNOPSIS
 DESCRIPTION
 FUNCTIONS
 eigsys
 matinv
 polyfit
 polycoef
 polyvalue
 detslatec
 PDL::Slatec::fft
 PDL::Slatec::rfft
 svdc
 poco
 geco
 gefa
 podi
 gedi
 gesl
 rs
 ezffti
 ezfftf
 ezfftb
 pcoef
 pvalue
 chim
 chic
 chsp
 chfd
 chfe
 chia
 chid
 chcm
 chbs
 polfit
 AUTHOR
NAME
PDL::Slatec  PDL interface to the slatec numerical programming library
SYNOPSIS
use PDL::Slatec;
($ndeg, $r, $ierr, $a) = polyfit($x, $y, $w, $maxdeg, $eps);
DESCRIPTION
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.
Piecewise cubic Hermite interpolation (PCHIP)
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. 238246), 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.
=cut
FUNCTIONS
eigsys
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!
matinv
Inverse of a square matrix
($inv) = matinv($mat)
polyfit
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 correctlyone 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 inplace 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.
C<$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.
polycoef
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);
polyvalue
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);
detslatec
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 noninvertible. 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.
PDL::Slatec::fft
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.
PDL::Slatec::rfft
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.
svdc
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
poco
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
geco
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
gefa
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
podi
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
gedi
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
gesl
Signature: (a(lda,n);int ipvt(n);b(n);int job())
Solve the real system A*X=B
or TRANS(A)*X=B
using the
factors computed by geco or gefa.
gesl does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
rs
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 badvalue flag of all output piddles if the flag is set for any of the input piddles.
ezffti
Signature: (int n();[o]wsave(foo))
Subroutine ezffti initializes the work array wsave()
which is used in both ezfftf and
ezfftb.
The prime factorization
of n
together with a tabulation of the trigonometric functions
are computed and stored in wsave()
.
ezffti does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
ezfftf
Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo))
ezfftf does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
ezfftb
Signature: ([o]r(n);azero();a(n);b(n);wsave(foo))
ezfftb does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
pcoef
Signature: (int l();c();[o]tc(bar);a(foo))
Convert the polfit
coefficients to Taylor series form.
c
and a()
must be of the same type.
pcoef does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
pvalue
Signature: (int l();x();[o]yfit();[o]yp(nder);a(foo))
Use the coefficients generated by polfit
to evaluate the
polynomial fit of degree l
, along with the first nder
of
its derivatives, at a specified point. x
and a
must be of the
same type.
pvalue does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chim
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 ($x,$f
,
where $x
is strictly increasing).
The resulting set of points  $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 $ierr
:

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.
chim does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chic
Signature: (int ic(two=2);vc(two=2);mflag();x(n);f(n);[o]d(n);wk(nwk);int [o]ierr())
Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.
Calculate the derivatives at the given points ($x,$f
,
where $x
is strictly increasing).
Control over the boundary conditions is given by the
$ic
and $vc
piddles, and the value of $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 $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 ic(0)
are:

1 if first derivative at
x(0)
is given invc(0)
. 
2 if second derivative at
x(0)
is given invc(0)
. 
3 to use the 3point difference formula for
d(0)
. (Reverts to the default b.c. ifn < 3
) 
4 to use the 4point difference formula for
d(0)
. (Reverts to the default b.c. ifn < 4
) 
5 to set
d(0)
so that the second derivative is continuous atx(1)
. (Reverts to the default b.c. ifn < 4
)
The values for ic(1)
are the same as above, except that
the firstderivative 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 nonzero to
use a formula based on the 3point 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
andd(0)
had to be adjusted for monotonicity. 
2 if
ic(1) < 0
andd(n1)
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*(n1)
.
chic does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chsp
Signature: (int ic(two=2);vc(two=2);x(n);f(n);[o]d(n);wk(nwk);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 ($x,$f
), with the specified
boundary conditions.
Control over the boundary conditions is given by the
$ic
and $vc
piddles.
The resulting values  $x,$f,$d
 can
be used in all the functions which expect a cubic
Hermite function.
The first and second elements of $ic
determine the boundary
conditions at the start and end of the data respectively.
The allowed values for ic(0)
are:

0 to set
d(0)
so that the third derivative is continuous atx(1)
. 
1 if first derivative at
x(0)
is given invc(0
). 
2 if second derivative at
x(0
) is given invc(0)
. 
3 to use the 3point difference formula for
d(0)
. (Reverts to the default b.c. ifn < 3
.) 
4 to use the 4point difference formula for
d(0)
. (Reverts to the default b.c. ifn < 4
.)
The values for ic(1)
are the same as above, except that
the firstderivative 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
oric(0) > 4
. 
5 if
ic(1) < 0
oric(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.
chsp does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chfd
Signature: (x(n);f(n);d(n);int check();xe(ne);[o]fe(ne);[o]de(ne);int [o]ierr())
Interpolate function and derivative values.
Given a piecewise cubic Hermite function  such as from
chim  evaluate the function ($fe
) and
derivative ($de
) at a set of points ($xe
).
If function values alone are required, use chfe.
Set check
to 0 to skip checks on the input data.
Error status returned by $ierr
:

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 lowerlevel routine, which should never happen.
chfd does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chfe
Signature: (x(n);f(n);d(n);int check();xe(ne);[o]fe(ne);int [o]ierr())
Interpolate function values.
Given a piecewise cubic Hermite function  such as from
chim  evaluate the function ($fe
) at
a set of points ($xe
).
If derivative values are also required, use chfd.
Set check
to 0 to skip checks on the input data.
Error status returned by $ierr
:

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
.
chfe does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chia
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 [$a,$b]
.
See chid if the integration limits are
data points.
Set check
to 0 to skip checks on the input data.
The values of $a
and $b
do not have
to lie within $x
, although the resulting integral
value will be highly suspect if they are not.
Error status returned by $ierr
:

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 lowerlevel routine, which should never happen.
chia does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chid
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 x($ia)
and
x($ib)
.
See chia for integration between arbitrary limits.
Although using a fortran routine, the values of
$ia
and $ib
are zero offset.
Set check
to 0 to skip checks on the input data.
Error status returned by $ierr
:

0 if successful.

1 if
nelem($x) < 2
. 
3 if
$x
is not strictly increasing. 
4 if
$ia
or$ib
is out of range.
chid does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chcm
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 $ismon
indicates over
which intervals the function is monotonic.
Set check
to 0 to skip checks on the input data.
For the data interval [x(i),x(i+1)]
, the
values of ismon(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 nonmonotonic

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 nonmonotonicity, 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.
chcm does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
chbs
Signature: (x(n);f(n);d(n);int knotyp();int nknots();t(tsize);[o]bcoef(bsize);int [o]ndim();int [o]kord();int [o]ierr())
Piecewise Cubic Hermite function to BSpline converter.
The resulting Bspline representation of the data
(i.e. nknots
, t
, bcoeff
, ndim
, and
kord
) can be evaluated by bvalu
(which is
currently not available).
Array sizes: tsize = 2*n + 4
, bsize = 2*n
,
and ndim = 2*n
.
knotyp
is a flag which controls the knot sequence.
The knot sequence t
is normally computed from $x
by putting a double knot at each x
and setting the end knot pairs
according to the value of knotyp
(where m = ndim = 2*n
):

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))
andt(m+3) = t(m+2) = x(n1) + (x(n1)x(n2))

2  Periodic placement of boundary knots:
t( 0 ) = t( 1 ) = x(0)  (x(n1)x(n2))
andt(m+3) = t(m+2) = x(n) + (x(1)x(0))

<0  Assume the
nknots
andt
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 Brepresentation
and may be changed by the routine.
If knotyp >= 0
, t
will be changed so that the
interior double knots are equal to the xvalues 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
andnknots != 2*n + 4
.
chbs does not process bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
polfit
Signature: (x(n); y(n); w(n); int maxdeg(); int [o]ndeg(); [o]eps(); [o]r(n); int [o]ierr(); [o]a(foo); [o]coeffs(bar);[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. x()
, y()
and w()
must be of the same type.
This version handles bad values appropriately
polfit processes bad values. It will set the badvalue flag of all output piddles if the flag is set for any of the input piddles.
AUTHOR
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.