Slatec(3) User Contributed Perl Documentation Slatec(3)NAMEPDL::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. 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.
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 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) = 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 $x variable and creates an appropriate $y variable
containing indices of the non-bad members of $x before calling the
Slatec routine "polfit".
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 non-invertible. The matrix threads
as usual.
This routine was previously known as "det" 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.
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
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.
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.
gefa
Signature: (a(n,n);int [o]ipvt(n);int [o]info())
Factor a matrix using Gaussian elimination.
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.
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.
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.
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.
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()".
ezfftf
Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo))
ezfftb
Signature: ([o]r(n);azero();a(n);b(n);wsave(foo))
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.
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.
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.
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 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)".
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 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.
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 lower-level routine, which should
never happen.
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".
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 lower-level routine, which should
never happen.
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.
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 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.
polfit
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. "x()", "y()" and "w()" must be of the same
type. This version handles bad values appropriately
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.
perl v5.10.0 2008-08-29 Slatec(3)
