1Slatec(3) User Contributed Perl Documentation Slatec(3)
2
6 PDL::Slatec - PDL interface to the slatec numerical programming library
7
9 use PDL::Slatec;
10 ($ndeg, $r, $ierr, $a) = polyfit($x, $y, $w, $maxdeg, $eps);
12
14 This module serves the dual purpose of providing an interface to parts
15 of the slatec library and showing how to interface PDL to an external
16 library. Using this library requires a fortran compiler; the source
17 for the routines is provided for convenience.
18 Currently available are routines to: manipulate matrices; calculate
20 FFT's; fit data using polynomials; and interpolate/integrate data using
21 piecewise cubic Hermite interpolation.
22 Piecewise cubic Hermite interpolation (PCHIP)
24 PCHIP is the slatec package of routines to perform piecewise cubic
25 Hermite interpolation of data. It features software to produce a
26 monotone and "visually pleasing" interpolant to monotone data.
27 According to Fritsch & Carlson ("Monotone piecewise cubic
28 interpolation", SIAM Journal on Numerical Analysis 17, 2 (April 1980),
29 pp. 238-246), such an interpolant may be more reasonable than a cubic
30 spline if the data contains both "steep" and "flat" sections.
31 Interpolation of cumulative probability distribution functions is
32 another application. These routines are cryptically named (blame
33 FORTRAN), beginning with 'ch', and accept either float or double
34 piddles.
35 Most of the routines require an integer parameter called "check"; if
37 set to 0, then no checks on the validity of the input data are made,
38 otherwise these checks are made. The value of "check" can be set to 0
39 if a routine such as chim has already been successfully called.
40 · If not known, estimate derivative values for the points using the
42 chim, chic, or chsp routines (the following routines require both
43 the function ("f") and derivative ("d") values at a set of points
44 ("x")).
45 · Evaluate, integrate, or differentiate the resulting PCH function
47 using the routines: chfd; chfe; chia; chid.
48 · If desired, you can check the monotonicity of your data using chcm.
50
52 eigsys
53 Eigenvalues and eigenvectors of a real positive definite symmetric
54 matrix.
55 ($eigvals,$eigvecs) = eigsys($mat)
57 Note: this function should be extended to calculate only eigenvalues if
59 called in scalar context!
60 matinv
62 Inverse of a square matrix
63 ($inv) = matinv($mat)
65 polyfit
67 Convenience wrapper routine about the "polfit" "slatec" function.
68 Separates supplied arguments and return values.
69 Fit discrete data in a least squares sense by polynomials in one
71 variable. Handles threading correctly--one can pass in a 2D PDL (as
72 $y) and it will pass back a 2D PDL, the rows of which are the
73 polynomial regression results (in $r corresponding to the rows of $y.
74 ($ndeg, $r, $ierr, $a) = polyfit($x, $y, $w, $maxdeg, $eps);
76 where on input:
78 C<x> and C<y> are the values to fit to a polynomial.
80 C<w> are weighting factors
81 C<maxdeg> is the maximum degree of polynomial to use and
82 C<eps> is the required degree of fit.
83 and on output:
85 C<ndeg> is the degree of polynomial actually used
87 C<r> is the values of the fitted polynomial
88 C<ierr> is a return status code, and
89 C<a> is some working array or other
90 C<eps> is modified to contain the rms error of the fit.
91 This version of polyfit handles bad values correctly. It strips them
93 out of the $x variable and creates an appropriate $y variable
94 containing indices of the non-bad members of $x before calling the
95 Slatec routine "polfit".
96 polycoef
98 Convenience wrapper routine around the "pcoef" "slatec" function.
99 Separates supplied arguments and return values.
100 Convert the "polyfit"/"polfit" coefficients to Taylor series form.
102 $tc = polycoef($l, $c, $a);
104 polyvalue
106 Convenience wrapper routine around the "pvalue" "slatec" function.
107 Separates supplied arguments and return values.
108 For multiple input x positions, a corresponding y position is
110 calculated.
111 The derivatives PDL is one dimensional (of size "nder") if a single x
113 position is supplied, two dimensional if more than one x position is
114 supplied.
115 Use the coefficients generated by "polyfit" (or "polfit") to evaluate
117 the polynomial fit of degree "l", along with the first "nder" of its
118 derivatives, at a specified point.
119 ($yfit, $yp) = polyvalue($l, $nder, $x, $a);
121 detslatec
123 compute the determinant of an invertible matrix
124 $mat = zeroes(5,5); $mat->diagonal(0,1) .= 1; # unity matrix
126 $det = detslatec $mat;
127 Usage:
129 $determinant = detslatec $matrix;
131 Signature: detslatec(mat(n,m); [o] det())
133 "detslatec" computes the determinant of an invertible matrix and barfs
135 if the matrix argument provided is non-invertible. The matrix threads
136 as usual.
137 This routine was previously known as "det" which clashes now with det
139 which is provided by PDL::MatrixOps. For the moment PDL::Slatec will
140 also load PDL::MatrixOps thereby making sure that older scripts work.
141 PDL::Slatec::fft
143 Fast Fourier Transform
144 $v_in = pdl(1,0,1,0);
146 ($azero,$a,$b) = PDL::Slatec::fft($v_in);
147 "PDL::Slatec::fft" is a convenience wrapper for ezfftf, and performs a
149 Fast Fourier Transform on an input vector $v_in. The return values are
150 the same as for ezfftf.
151 PDL::Slatec::rfft
153 reverse Fast Fourier Transform
154 $v_out = PDL::Slatec::rfft($azero,$a,$b);
156 print $v_in, $vout
157 [1 0 1 0] [1 0 1 0]
158 "PDL::Slatec::rfft" is a convenience wrapper for ezfftb, and performs a
160 reverse Fast Fourier Transform. The input is the same as the output of
161 PDL::Slatec::fft, and the output of "rfft" is a data vector, similar to
162 what is input into PDL::Slatec::fft.
163 svdc
165 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())
166 singular value decomposition of a matrix
168 svdc does not process bad values. It will set the bad-value flag of
170 all output piddles if the flag is set for any of the input piddles.
171 poco
173 Signature: (a(n,n);rcond();[o]z(n);int [o]info())
174 Factor a real symmetric positive definite matrix and estimate the
176 condition number of the matrix.
177 poco does not process bad values. It will set the bad-value flag of
179 all output piddles if the flag is set for any of the input piddles.
180 geco
182 Signature: (a(n,n);int [o]ipvt(n);[o]rcond();[o]z(n))
183 Factor a matrix using Gaussian elimination and estimate the condition
185 number of the matrix.
186 geco does not process bad values. It will set the bad-value flag of
188 all output piddles if the flag is set for any of the input piddles.
189 gefa
191 Signature: (a(n,n);int [o]ipvt(n);int [o]info())
192 Factor a matrix using Gaussian elimination.
194 gefa does not process bad values. It will set the bad-value flag of
196 all output piddles if the flag is set for any of the input piddles.
197 podi
199 Signature: (a(n,n);[o]det(two=2);int job())
200 Compute the determinant and inverse of a certain real symmetric
202 positive definite matrix using the factors computed by poco.
203 podi does not process bad values. It will set the bad-value flag of
205 all output piddles if the flag is set for any of the input piddles.
206 gedi
208 Signature: (a(n,n);int [o]ipvt(n);[o]det(two=2);[o]work(n);int job())
209 Compute the determinant and inverse of a matrix using the factors
211 computed by geco or gefa.
212 gedi does not process bad values. It will set the bad-value flag of
214 all output piddles if the flag is set for any of the input piddles.
215 gesl
217 Signature: (a(lda,n);int ipvt(n);b(n);int job())
218 Solve the real system "A*X=B" or "TRANS(A)*X=B" using the factors
220 computed by geco or gefa.
221 gesl does not process bad values. It will set the bad-value flag of
223 all output piddles if the flag is set for any of the input piddles.
224 rs
226 Signature: (a(n,n);[o]w(n);int matz();[o]z(n,n);[t]fvone(n);[t]fvtwo(n);int [o]ierr())
227 This subroutine calls the recommended sequence of subroutines from the
229 eigensystem subroutine package (EISPACK) to find the eigenvalues and
230 eigenvectors (if desired) of a REAL SYMMETRIC matrix.
231 rs does not process bad values. It will set the bad-value flag of all
233 output piddles if the flag is set for any of the input piddles.
234 ezffti
236 Signature: (int n();[o]wsave(foo))
237 Subroutine ezffti initializes the work array "wsave()" which is used in
239 both ezfftf and ezfftb. The prime factorization of "n" together with a
240 tabulation of the trigonometric functions are computed and stored in
241 "wsave()".
242 ezffti does not process bad values. It will set the bad-value flag of
244 all output piddles if the flag is set for any of the input piddles.
245 ezfftf
247 Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo))
248 ezfftf does not process bad values. It will set the bad-value flag of
250 all output piddles if the flag is set for any of the input piddles.
251 ezfftb
253 Signature: ([o]r(n);azero();a(n);b(n);wsave(foo))
254 ezfftb does not process bad values. It will set the bad-value flag of
256 all output piddles if the flag is set for any of the input piddles.
257 pcoef
259 Signature: (int l();c();[o]tc(bar);a(foo))
260 Convert the "polfit" coefficients to Taylor series form. "c" and "a()"
262 must be of the same type.
263 pcoef does not process bad values. It will set the bad-value flag of
265 all output piddles if the flag is set for any of the input piddles.
266 pvalue
268 Signature: (int l();x();[o]yfit();[o]yp(nder);a(foo))
269 Use the coefficients generated by "polfit" to evaluate the polynomial
271 fit of degree "l", along with the first "nder" of its derivatives, at a
272 specified point. "x" and "a" must be of the same type.
273 pvalue does not process bad values. It will set the bad-value flag of
275 all output piddles if the flag is set for any of the input piddles.
276 chim
278 Signature: (x(n);f(n);[o]d(n);int [o]ierr())
279 Calculate the derivatives of (x,f(x)) using cubic Hermite
281 interpolation.
282 Calculate the derivatives at the given set of points ("$x,$f", where $x
284 is strictly increasing). The resulting set of points - "$x,$f,$d",
285 referred to as the cubic Hermite representation - can then be used in
286 other functions, such as chfe, chfd, and chia.
287 The boundary conditions are compatible with monotonicity, and if the
289 data are only piecewise monotonic, the interpolant will have an
290 extremum at the switch points; for more control over these issues use
291 chic.
292 Error status returned by $ierr:
294 · 0 if successful.
296 · > 0 if there were "ierr" switches in the direction of monotonicity
298 (data still valid).
299 · -1 if "nelem($x) < 2".
301 · -3 if $x is not strictly increasing.
303 chim does not process bad values. It will set the bad-value flag of
305 all output piddles if the flag is set for any of the input piddles.
306 chic
308 Signature: (int ic(two=2);vc(two=2);mflag();x(n);f(n);[o]d(n);wk(nwk);int [o]ierr())
309 Calculate the derivatives of (x,f(x)) using cubic Hermite
311 interpolation.
312 Calculate the derivatives at the given points ("$x,$f", where $x is
314 strictly increasing). Control over the boundary conditions is given by
315 the $ic and $vc piddles, and the value of $mflag determines the
316 treatment of points where monotoncity switches direction. A simpler,
317 more restricted, interface is available using chim.
318 The first and second elements of $ic determine the boundary conditions
320 at the start and end of the data respectively. If the value is 0, then
321 the default condition, as used by chim, is adopted. If greater than
322 zero, no adjustment for monotonicity is made, otherwise if less than
323 zero the derivative will be adjusted. The allowed magnitudes for ic(0)
324 are:
325 · 1 if first derivative at x(0) is given in vc(0).
327 · 2 if second derivative at x(0) is given in vc(0).
329 · 3 to use the 3-point difference formula for d(0). (Reverts to the
331 default b.c. if "n < 3")
332 · 4 to use the 4-point difference formula for d(0). (Reverts to the
334 default b.c. if "n < 4")
335 · 5 to set d(0) so that the second derivative is continuous at x(1).
337 (Reverts to the default b.c. if "n < 4")
338 The values for ic(1) are the same as above, except that the first-
340 derivative value is stored in vc(1) for cases 1 and 2. The values of
341 $vc need only be set if options 1 or 2 are chosen for $ic.
342 Set "$mflag = 0" if interpolant is required to be monotonic in each
344 interval, regardless of the data. This causes $d to be set to 0 at all
345 switch points. Set $mflag to be non-zero to use a formula based on the
346 3-point difference formula at switch points. If "$mflag > 0", then the
347 interpolant at swich points is forced to not deviate from the data by
348 more than "$mflag*dfloc", where "dfloc" is the maximum of the change of
349 $f on this interval and its two immediate neighbours. If "$mflag < 0",
350 no such control is to be imposed.
351 The piddle $wk is only needed for work space. However, I could not get
353 it to work as a temporary variable, so you must supply it; it is a 1D
354 piddle with "2*n" elements.
355 Error status returned by $ierr:
357 · 0 if successful.
359 · 1 if "ic(0) < 0" and d(0) had to be adjusted for monotonicity.
361 · 2 if "ic(1) < 0" and "d(n-1)" had to be adjusted for monotonicity.
363 · 3 if both 1 and 2 are true.
365 · -1 if "n < 2".
367 · -3 if $x is not strictly increasing.
369 · -4 if "abs(ic(0)) > 5".
371 · -5 if "abs(ic(1)) > 5".
373 · -6 if both -4 and -5 are true.
375 · -7 if "nwk < 2*(n-1)".
377 chic does not process bad values. It will set the bad-value flag of
379 all output piddles if the flag is set for any of the input piddles.
380 chsp
382 Signature: (int ic(two=2);vc(two=2);x(n);f(n);[o]d(n);wk(nwk);int [o]ierr())
383 Calculate the derivatives of (x,f(x)) using cubic spline interpolation.
385 Calculate the derivatives, using cubic spline interpolation, at the
387 given points ("$x,$f"), with the specified boundary conditions.
388 Control over the boundary conditions is given by the $ic and $vc
389 piddles. The resulting values - "$x,$f,$d" - can be used in all the
390 functions which expect a cubic Hermite function.
391 The first and second elements of $ic determine the boundary conditions
393 at the start and end of the data respectively. The allowed values for
394 ic(0) are:
395 · 0 to set d(0) so that the third derivative is continuous at x(1).
397 · 1 if first derivative at x(0) is given in "vc(0").
399 · 2 if second derivative at "x(0") is given in vc(0).
401 · 3 to use the 3-point difference formula for d(0). (Reverts to the
403 default b.c. if "n < 3".)
404 · 4 to use the 4-point difference formula for d(0). (Reverts to the
406 default b.c. if "n < 4".)
407 The values for ic(1) are the same as above, except that the first-
409 derivative value is stored in vc(1) for cases 1 and 2. The values of
410 $vc need only be set if options 1 or 2 are chosen for $ic.
411 The piddle $wk is only needed for work space. However, I could not get
413 it to work as a temporary variable, so you must supply it; it is a 1D
414 piddle with "2*n" elements.
415 Error status returned by $ierr:
417 · 0 if successful.
419 · -1 if "nelem($x) < 2".
421 · -3 if $x is not strictly increasing.
423 · -4 if "ic(0) < 0" or "ic(0) > 4".
425 · -5 if "ic(1) < 0" or "ic(1) > 4".
427 · -6 if both of the above are true.
429 · -7 if "nwk < 2*n".
431 · -8 in case of trouble solving the linear system for the interior
433 derivative values.
434 chsp does not process bad values. It will set the bad-value flag of
436 all output piddles if the flag is set for any of the input piddles.
437 chfd
439 Signature: (x(n);f(n);d(n);int check();xe(ne);[o]fe(ne);[o]de(ne);int [o]ierr())
440 Interpolate function and derivative values.
442 Given a piecewise cubic Hermite function - such as from chim - evaluate
444 the function ($fe) and derivative ($de) at a set of points ($xe). If
445 function values alone are required, use chfe. Set "check" to 0 to skip
446 checks on the input data.
447 Error status returned by $ierr:
449 · 0 if successful.
451 · >0 if extrapolation was performed at "ierr" points (data still
453 valid).
454 · -1 if "nelem($x) < 2"
456 · -3 if $x is not strictly increasing.
458 · -4 if "nelem($xe) < 1".
460 · -5 if an error has occurred in a lower-level routine, which should
462 never happen.
463 chfd does not process bad values. It will set the bad-value flag of
465 all output piddles if the flag is set for any of the input piddles.
466 chfe
468 Signature: (x(n);f(n);d(n);int check();xe(ne);[o]fe(ne);int [o]ierr())
469 Interpolate function values.
471 Given a piecewise cubic Hermite function - such as from chim - evaluate
473 the function ($fe) at a set of points ($xe). If derivative values are
474 also required, use chfd. Set "check" to 0 to skip checks on the input
475 data.
476 Error status returned by $ierr:
478 · 0 if successful.
480 · >0 if extrapolation was performed at "ierr" points (data still
482 valid).
483 · -1 if "nelem($x) < 2"
485 · -3 if $x is not strictly increasing.
487 · -4 if "nelem($xe) < 1".
489 chfe does not process bad values. It will set the bad-value flag of
491 all output piddles if the flag is set for any of the input piddles.
492 chia
494 Signature: (x(n);f(n);d(n);int check();a();b();[o]ans();int [o]ierr())
495 Integrate (x,f(x)) over arbitrary limits.
497 Evaluate the definite integral of a a piecewise cubic Hermite function
499 over an arbitrary interval, given by "[$a,$b]". See chid if the
500 integration limits are data points. Set "check" to 0 to skip checks on
501 the input data.
502 The values of $a and $b do not have to lie within $x, although the
504 resulting integral value will be highly suspect if they are not.
505 Error status returned by $ierr:
507 · 0 if successful.
509 · 1 if $a lies outside $x.
511 · 2 if $b lies outside $x.
513 · 3 if both 1 and 2 are true.
515 · -1 if "nelem($x) < 2"
517 · -3 if $x is not strictly increasing.
519 · -4 if an error has occurred in a lower-level routine, which should
521 never happen.
522 chia does not process bad values. It will set the bad-value flag of
524 all output piddles if the flag is set for any of the input piddles.
525 chid
527 Signature: (x(n);f(n);d(n);int check();int ia();int ib();[o]ans();int [o]ierr())
528 Integrate (x,f(x)) between data points.
530 Evaluate the definite integral of a a piecewise cubic Hermite function
532 between "x($ia)" and "x($ib)".
533 See chia for integration between arbitrary limits.
535 Although using a fortran routine, the values of $ia and $ib are zero
537 offset. Set "check" to 0 to skip checks on the input data.
538 Error status returned by $ierr:
540 · 0 if successful.
542 · -1 if "nelem($x) < 2".
544 · -3 if $x is not strictly increasing.
546 · -4 if $ia or $ib is out of range.
548 chid does not process bad values. It will set the bad-value flag of
550 all output piddles if the flag is set for any of the input piddles.
551 chcm
553 Signature: (x(n);f(n);d(n);int check();int [o]ismon(n);int [o]ierr())
554 Check the given piecewise cubic Hermite function for monotonicity.
556 The outout piddle $ismon indicates over which intervals the function is
558 monotonic. Set "check" to 0 to skip checks on the input data.
559 For the data interval "[x(i),x(i+1)]", the values of ismon(i) can be:
561 · -3 if function is probably decreasing
563 · -1 if function is strictly decreasing
565 · 0 if function is constant
567 · 1 if function is strictly increasing
569 · 2 if function is non-monotonic
571 · 3 if function is probably increasing
573 If "abs(ismon(i)) == 3", the derivative values are near the boundary of
575 the monotonicity region. A small increase produces non-monotonicity,
576 whereas a decrease produces strict monotonicity.
577 The above applies to "i = 0 .. nelem($x)-1". The last element of $ismon
579 indicates whether the entire function is monotonic over $x.
580 Error status returned by $ierr:
582 · 0 if successful.
584 · -1 if "n < 2".
586 · -3 if $x is not strictly increasing.
588 chcm does not process bad values. It will set the bad-value flag of
590 all output piddles if the flag is set for any of the input piddles.
591 polfit
593 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))
594 Fit discrete data in a least squares sense by polynomials
596 in one variable. "x()", "y()" and "w()" must be of the same
597 type. This version handles bad values appropriately
598 polfit does handle bad values. It will set the bad-value flag of all
600 output piddles if the flag is set for any of the input piddles.
601
603 Copyright (C) 1997 Tuomas J. Lukka. Copyright (C) 2000 Tim Jenness,
604 Doug Burke. All rights reserved. There is no warranty. You are allowed
605 to redistribute this software / documentation under certain conditions.
606 For details, see the file COPYING in the PDL distribution. If this file
607 is separated from the PDL distribution, the copyright notice should be
608 included in the file.
609perl v5.12.3 2011-03-31 Slatec(3)