oracular (3) PDL::Slatec.3pm.gz
NAME
PDL::Slatec - PDL interface to the slatec numerical programming library
SYNOPSIS
use PDL::Slatec;
($ndeg, $r, $ierr, $c) = 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. 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 ndarrays.
Most of the routines require an integer parameter called "skip"; if set to 1, then no checks on the
validity of the input data are made, otherwise these checks are made. The value of "skip" can be set to
1 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", or use "chbs" to convert a PCH
function into B-representation for use with the B-spline routines of SLATEC.
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 broadcasting
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, $c, $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 $c 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.
$c is a working array accessible to Slatec - you can feed it to several other Slatec routines to get nice
things out. It does not broadcast 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, $x);
detslatec
compute the determinant of an invertible matrix
$mat = identity(5); # 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 broadcasts as usual.
This routine was previously known as "det" which clashes now with det which is provided by
PDL::MatrixOps.
fft
Fast Fourier Transform
$v_in = pdl(1,0,1,0);
($azero,$x,$y) = 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".
rfft
reverse Fast Fourier Transform
$v_out = PDL::Slatec::rfft($azero,$x,$y);
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);[t]work(n);longlong job();longlong [o]info())
singular value decomposition of a matrix
svdc does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
poco
Signature: ([io]a(n,n);[o]rcond();[o]z(n);longlong [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 ndarrays if the flag is
set for any of the input ndarrays.
geco
Signature: (a(n,n);longlong [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 ndarrays if the flag is
set for any of the input ndarrays.
gefa
Signature: ([io]a(n,n);longlong [o]ipvt(n);longlong [o]info())
Factor a matrix using Gaussian elimination.
gefa does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
podi
Signature: ([io]a(n,n);[o]det(two=2);longlong 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 ndarrays if the flag is
set for any of the input ndarrays.
gedi
Signature: ([io]a(n,n);longlong ipvt(n);[o]det(two=2);[t]work(n);longlong 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 ndarrays if the flag is
set for any of the input ndarrays.
gesl
Signature: (a(lda,n);longlong ipvt(n);[io]b(n);longlong 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 bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
rs
Signature: (a(n,n);[o]w(n);longlong matz();[o]z(n,n);[t]fvone(n);[t]fvtwo(n);longlong [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 ndarrays if the flag is set
for any of the input ndarrays.
ezffti
Signature: (longlong n();[io]wsave(foo);longlong [o]ifac(ni=15))
Subroutine ezffti initializes the work array "wsave(3n or more)" and ifac() 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 bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
ezfftf
Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo);longlong ifac(ni=15))
ezfftf does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
ezfftb
Signature: ([o]r(n);azero();a(n);b(n);wsave(foo);longlong ifac(ni=15))
ezfftb does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
pcoef
Signature: (longlong 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 bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
polyvalue
Signature: (longlong l();x();[o]yfit();[o]yp(nder);a(foo); PDL_LongLong nder => nder)
Use the coefficients "c" 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 "x". Broadcasts correctly over
multiple "x" positions.
($yfit, $yp) = polyvalue($l, $nder, $x, $c);
polyvalue does not process bad values. It will set the bad-value flag of all output ndarrays if the flag
is set for any of the input ndarrays.
chim
Signature: (x(n);f(n);[o]d(n);longlong [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 bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chic
Signature: (longlong ic(two=2);vc(two=2);mflag();x(n);f(n);[o]d(n);[t]wk(nwk=CALC(2*$SIZE(n)));longlong [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 ndarrays, 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 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.
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)".
chic does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chsp
Signature: (longlong ic(two=2);vc(two=2);x(n);f(n);[o]d(n);[t]wk(nwk=CALC(2*$SIZE(n)));longlong [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
ndarrays. 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 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.
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.
chsp does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chfd
Signature: (x(n);f(n);d(n);int [io]skip();xe(ne);[o]fe(ne);[o]de(ne);longlong [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
"skip" 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 lower-level routine, which should never happen.
chfd does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chfe
Signature: (x(n);f(n);d(n);int [io]skip();xe(ne);[o]fe(ne);longlong [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 "skip" 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 bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chia
Signature: (x(n);f(n);d(n);int [io]skip();la();lb();[o]ans();longlong [o]ierr())
Integrate (x,f(x)) over arbitrary limits.
Evaluate the definite integral of a piecewise cubic Hermite function over an arbitrary interval, given by
"[$la,$lb]". $d should contain the derivative values, computed by "chim". See "chid" if the integration
limits are data points. Set "skip" to 0 to skip checks on the input data.
The values of $la and $lb 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 $la lies outside $x.
• 2 if $lb 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.
chia does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chid
Signature: (x(n);f(n);d(n);int [io]skip();longlong ia();longlong ib();[o]ans();longlong [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. $d should contain the
derivative values, computed by "chim". Set "skip" 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 bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chcm
Signature: (x(n);f(n);d(n);int [io]skip();longlong [o]ismon(n);longlong [o]ierr())
Check the given piecewise cubic Hermite function for monotonicity.
The outout ndarray $ismon indicates over which intervals the function is monotonic. Set "skip" 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 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.
chcm does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
chbs
Signature: (x(n);f(n);d(n);longlong knotyp();longlong [io]nknots();[io]t(tsize=CALC(2*$SIZE(n)+4));[o]bcoef(bsize=CALC(2*$SIZE(n)));longlong [o]ndim();longlong [o]kord();longlong [o]ierr())
Piecewise Cubic Hermite function to B-Spline converter.
The resulting B-spline representation of the data (i.e. "nknots", "t", "bcoeff", "ndim", and "kord") can
be evaluated by "bvalu".
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))" 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".
chbs does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
bvalu
Signature: (t(nplusk);a(n);longlong ideriv();x();[t]work(k3=CALC(3*($SIZE(nplusk)-$SIZE(n))));[o]ans())
Evaluate the B-representation of a B-spline at X for the function value or any of its derivatives.
Evaluates the B-representation "(T,A,N,K)" of a B-spline at "X" for the function value on "IDERIV = 0" or
any of its derivatives on "IDERIV = 1,2,...,K-1". Right limiting values (right derivatives) are returned
except at the right end point "X=T(N+1)" where left limiting values are computed. The spline is defined
on "T(K) .LE. X .LE. T(N+1)". BVALU returns a fatal error message when "X" is outside of this interval.
To compute left derivatives or left limiting values at a knot T(I), replace "N" by "I-1" and set
"X=T(I)", "I=K+1,N+1".
Description of Arguments
Input
T - knot vector of length N+K
A - B-spline coefficient vector of length N
N - number of B-spline coefficients
N = sum of knot multiplicities-K
K - order of the B-spline, K .GE. 1
IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1
IDERIV=0 returns the B-spline value
X - argument, T(K) .LE. X .LE. T(N+1)
INBV - an initialization parameter which must be set
to 1 the first time BVALU is called. (PDL sets to 1 for you)
Output
INBV - INBV contains information for efficient process-
ing after the initial call and INBV must not
be changed by the user. Distinct splines require
distinct INBV parameters. (unavailable in PDL)
WORK - work vector of length 3*K.
BVALU - value of the IDERIV-th derivative at X
Error Conditions
An improper input is a fatal error
bvalu does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is
set for any of the input ndarrays.
polfit
Signature: (x(n); y(n); w(n); longlong maxdeg(); longlong [o]ndeg(); [io]eps(); [o]r(n); longlong [o]ierr(); [o]a(foo=CALC(3*($SIZE(n) + $SIZE(bar)))); [o]coeffs(bar);[t]xtmp(n);[t]ytmp(n);[t]wtmp(n);[t]rtmp(n); IV max_maxdeg_plus1 => bar)
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 bad-value flag of all output ndarrays if the flag is set
for any of the input ndarrays.
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.