NAME
PDL::Slatec - PDL interface to some LINPACK and EISPACK routines - DEPRECATED
SYNOPSIS
use PDL::Slatec;
($ndeg, $r, $ierr, $c) = polyfit($x, $y, $w, $maxdeg, $eps);DESCRIPTION
This module is now deprecated in favour of PDL::LinearAlgebra.
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; and fit data using polynomials.
Piecewise cubic Hermite interpolation (PCHIP)
These routines are now in PDL::Primitive as of PDL 2.096.
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!
This is the EISPACK routine rs.
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".
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);
EOF
); ################################################################## ##################################################################
# This version of polfit accepts bad values and also allows broadcasting
# indices: # n runs across input points; # foo runs across wacky Slatec buffer size; # bar runs across polynomial coefficients. pp_def('polfit', Pars => '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)', OtherPars => 'IV max_maxdeg_plus1 => bar', GenericTypes => ['F','D'], HandleBad => 1, NoBadifNaN => 1, CHeader => <<'EOF', void PDL_FORTRAN(dpcoef)(PDL_LongLong *l,PDL_Double *c,PDL_Double *tc,PDL_Double *a); void PDL_FORTRAN(pcoef)(PDL_LongLong *l,PDL_Float *c,PDL_Float *tc,PDL_Float *a); EOF Code => 'PDL_LongLong howmany = PDL_IF_BAD(0,$SIZE(n)), i; if($SIZE(n)<$maxdeg()) $CROAK("Must have at least <n> points to fit <n> coefficients"); PDL_IF_BAD( loop(n) %{ /* get rid of bad values. Call polfit with [xyw]tmp instead of [xyz]. */ if ($ISGOOD(y()) && $ISGOOD(x()) && $ISGOOD(w())) { $xtmp(n=>howmany) = $x(); $ytmp(n=>howmany) = $y(); $wtmp(n=>howmany) = $w(); howmany++; } %} if (howmany <= $maxdeg()) { /* Not enough good datapoints -- set this whole row to BAD. */ loop(n) %{ $SETBAD(r()); %} $ierr() = 2; } else,) { /* Found enough good datapoints for a fit -- do the fit */ $GENERIC() zero = 0; /* Do the fit */ PDL_FORTRAN($TFD(polfit,dpolft))( &howmany, PDL_IF_BAD($P(xtmp),$P(x)), PDL_IF_BAD($P(ytmp),$P(y)), PDL_IF_BAD($P(wtmp),$P(w)), $P(maxdeg),$P(ndeg),$P(eps), PDL_IF_BAD($P(rtmp),$P(r)), $P(ierr),$P(a) ); PDL_LongLong maxord = $a(foo=>0)+0.5; PDL_LongLong ord = PDLMIN($maxdeg(),$a(foo=>(maxord*3)+2)); /* Extract the polynomial coefficients into coeffs -- used for bad values */ PDL_FORTRAN($TFD(,d)pcoef)(&(ord), &(zero), $P(coeffs), $P(a)); for(i=ord+1; i<=$maxdeg(); i++) $coeffs(bar=>i) = 0; PDL_IF_BAD( howmany=0; loop(n) %{ /* replace bad values */ if ($ISGOOD(y())) { $r() = $rtmp(n=>howmany); howmany++; } else if($ISGOOD(x())) { /* Bad values are omitted from the call to polfit, so we must reconstitute them on return */ /* (just because a value is bad in y, does not mean the fit itself is bad there) */ /* */ /* The number in ord is not the number of coefficients in the polynomial, it is the highest */ /* order coefficient -- so 3 for a cubic, which has 4 coefficients. */ PDL_LongLong ii; $GENERIC() acc = 0; for( ii=ord; ii>0; ii-- ) { acc += $coeffs(bar=>ii); acc *= $x(); } acc += $coeffs(bar=>0); $r() = acc; } else { /* x and y are bad here... */ $SETBAD(r()); } %},) }', Doc => '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', );
################################################################## ##################################################################
pp_addpm(<<'EOD'); =head1 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.