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, $coeffs, $rms) = polyfit($x, $y, $w, $maxdeg, [$eps]);
76 $coeffs = polyfit($x,$y,$w,$maxdeg,[$eps]);
78 where on input:
80 $x and $y are the values to fit to a polynomial. $w are weighting
82 factors $maxdeg is the maximum degree of polynomial to use and $eps is
83 the required degree of fit.
84 and the output switches on list/scalar context.
86 In list context:
88 $ndeg is the degree of polynomial actually used $r is the values of the
90 fitted polynomial $ierr is a return status code, and $a is some working
91 array or other (preserved for historical purposes) $coeffs is the
92 polynomial coefficients of the best fit polynomial. $rms is the rms
93 error of the fit.
94 In scalar context, only $coeffs is returned.
96 Historically, $eps was modified in-place to be a return value of the
98 rms error. This usage is deprecated, and $eps is an optional parameter
99 now. It is still modified if present.
100 $a is a working array accessible to Slatec - you can feed it to several
102 other Slatec routines to get nice things out. It does not thread
103 correctly and should probably be fixed by someone. If you are reading
104 this, that someone might be you.
105 This version of polyfit handles bad values correctly. Bad values in $y
107 are ignored for the fit and give computed values on the fitted curve in
108 the return. Bad values in $x or $w are ignored for the fit and result
109 in bad elements in the output.
110 polycoef
112 Convenience wrapper routine around the "pcoef" "slatec" function.
113 Separates supplied arguments and return values.
114 Convert the "polyfit"/"polfit" coefficients to Taylor series form.
116 $tc = polycoef($l, $c, $a);
118 polyvalue
120 Convenience wrapper routine around the "pvalue" "slatec" function.
121 Separates supplied arguments and return values.
122 For multiple input x positions, a corresponding y position is
124 calculated.
125 The derivatives PDL is one dimensional (of size "nder") if a single x
127 position is supplied, two dimensional if more than one x position is
128 supplied.
129 Use the coefficients generated by "polyfit" (or "polfit") to evaluate
131 the polynomial fit of degree "l", along with the first "nder" of its
132 derivatives, at a specified point.
133 ($yfit, $yp) = polyvalue($l, $nder, $x, $a);
135 detslatec
137 compute the determinant of an invertible matrix
138 $mat = zeroes(5,5); $mat->diagonal(0,1) .= 1; # unity matrix
140 $det = detslatec $mat;
141 Usage:
143 $determinant = detslatec $matrix;
145 Signature: detslatec(mat(n,m); [o] det())
147 "detslatec" computes the determinant of an invertible matrix and barfs
149 if the matrix argument provided is non-invertible. The matrix threads
150 as usual.
151 This routine was previously known as "det" which clashes now with det
153 which is provided by PDL::MatrixOps.
154 fft
156 Fast Fourier Transform
157 $v_in = pdl(1,0,1,0);
159 ($azero,$a,$b) = PDL::Slatec::fft($v_in);
160 "PDL::Slatec::fft" is a convenience wrapper for ezfftf, and performs a
162 Fast Fourier Transform on an input vector $v_in. The return values are
163 the same as for ezfftf.
164 rfft
166 reverse Fast Fourier Transform
167 $v_out = PDL::Slatec::rfft($azero,$a,$b);
169 print $v_in, $vout
170 [1 0 1 0] [1 0 1 0]
171 "PDL::Slatec::rfft" is a convenience wrapper for ezfftb, and performs a
173 reverse Fast Fourier Transform. The input is the same as the output of
174 PDL::Slatec::fft, and the output of "rfft" is a data vector, similar to
175 what is input into PDL::Slatec::fft.
176 svdc
178 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())
179 singular value decomposition of a matrix
181 svdc does not process bad values. It will set the bad-value flag of
183 all output piddles if the flag is set for any of the input piddles.
184 poco
186 Signature: (a(n,n);rcond();[o]z(n);int [o]info())
187 Factor a real symmetric positive definite matrix and estimate the
189 condition number of the matrix.
190 poco does not process bad values. It will set the bad-value flag of
192 all output piddles if the flag is set for any of the input piddles.
193 geco
195 Signature: (a(n,n);int [o]ipvt(n);[o]rcond();[o]z(n))
196 Factor a matrix using Gaussian elimination and estimate the condition
198 number of the matrix.
199 geco does not process bad values. It will set the bad-value flag of
201 all output piddles if the flag is set for any of the input piddles.
202 gefa
204 Signature: (a(n,n);int [o]ipvt(n);int [o]info())
205 Factor a matrix using Gaussian elimination.
207 gefa does not process bad values. It will set the bad-value flag of
209 all output piddles if the flag is set for any of the input piddles.
210 podi
212 Signature: (a(n,n);[o]det(two=2);int job())
213 Compute the determinant and inverse of a certain real symmetric
215 positive definite matrix using the factors computed by poco.
216 podi does not process bad values. It will set the bad-value flag of
218 all output piddles if the flag is set for any of the input piddles.
219 gedi
221 Signature: (a(n,n);int [o]ipvt(n);[o]det(two=2);[o]work(n);int job())
222 Compute the determinant and inverse of a matrix using the factors
224 computed by geco or gefa.
225 gedi does not process bad values. It will set the bad-value flag of
227 all output piddles if the flag is set for any of the input piddles.
228 gesl
230 Signature: (a(lda,n);int ipvt(n);b(n);int job())
231 Solve the real system "A*X=B" or "TRANS(A)*X=B" using the factors
233 computed by geco or gefa.
234 gesl does not process bad values. It will set the bad-value flag of
236 all output piddles if the flag is set for any of the input piddles.
237 rs
239 Signature: (a(n,n);[o]w(n);int matz();[o]z(n,n);[t]fvone(n);[t]fvtwo(n);int [o]ierr())
240 This subroutine calls the recommended sequence of subroutines from the
242 eigensystem subroutine package (EISPACK) to find the eigenvalues and
243 eigenvectors (if desired) of a REAL SYMMETRIC matrix.
244 rs does not process bad values. It will set the bad-value flag of all
246 output piddles if the flag is set for any of the input piddles.
247 ezffti
249 Signature: (int n();[o]wsave(foo))
250 Subroutine ezffti initializes the work array "wsave()" which is used in
252 both ezfftf and ezfftb. The prime factorization of "n" together with a
253 tabulation of the trigonometric functions are computed and stored in
254 "wsave()".
255 ezffti does not process bad values. It will set the bad-value flag of
257 all output piddles if the flag is set for any of the input piddles.
258 ezfftf
260 Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo))
261 ezfftf does not process bad values. It will set the bad-value flag of
263 all output piddles if the flag is set for any of the input piddles.
264 ezfftb
266 Signature: ([o]r(n);azero();a(n);b(n);wsave(foo))
267 ezfftb does not process bad values. It will set the bad-value flag of
269 all output piddles if the flag is set for any of the input piddles.
270 pcoef
272 Signature: (int l();c();[o]tc(bar);a(foo))
273 Convert the "polfit" coefficients to Taylor series form. "c" and "a()"
275 must be of the same type.
276 pcoef does not process bad values. It will set the bad-value flag of
278 all output piddles if the flag is set for any of the input piddles.
279 pvalue
281 Signature: (int l();x();[o]yfit();[o]yp(nder);a(foo))
282 Use the coefficients generated by "polfit" to evaluate the polynomial
284 fit of degree "l", along with the first "nder" of its derivatives, at a
285 specified point. "x" and "a" must be of the same type.
286 pvalue does not process bad values. It will set the bad-value flag of
288 all output piddles if the flag is set for any of the input piddles.
289 chim
291 Signature: (x(n);f(n);[o]d(n);int [o]ierr())
292 Calculate the derivatives of (x,f(x)) using cubic Hermite
294 interpolation.
295 Calculate the derivatives at the given set of points ("$x,$f", where $x
297 is strictly increasing). The resulting set of points - "$x,$f,$d",
298 referred to as the cubic Hermite representation - can then be used in
299 other functions, such as chfe, chfd, and chia.
300 The boundary conditions are compatible with monotonicity, and if the
302 data are only piecewise monotonic, the interpolant will have an
303 extremum at the switch points; for more control over these issues use
304 chic.
305 Error status returned by $ierr:
307 · 0 if successful.
309 · > 0 if there were "ierr" switches in the direction of monotonicity
311 (data still valid).
312 · -1 if "nelem($x) < 2".
314 · -3 if $x is not strictly increasing.
316 chim does not process bad values. It will set the bad-value flag of
318 all output piddles if the flag is set for any of the input piddles.
319 chic
321 Signature: (int ic(two=2);vc(two=2);mflag();x(n);f(n);[o]d(n);wk(nwk);int [o]ierr())
322 Calculate the derivatives of (x,f(x)) using cubic Hermite
324 interpolation.
325 Calculate the derivatives at the given points ("$x,$f", where $x is
327 strictly increasing). Control over the boundary conditions is given by
328 the $ic and $vc piddles, and the value of $mflag determines the
329 treatment of points where monotoncity switches direction. A simpler,
330 more restricted, interface is available using chim.
331 The first and second elements of $ic determine the boundary conditions
333 at the start and end of the data respectively. If the value is 0, then
334 the default condition, as used by chim, is adopted. If greater than
335 zero, no adjustment for monotonicity is made, otherwise if less than
336 zero the derivative will be adjusted. The allowed magnitudes for ic(0)
337 are:
338 · 1 if first derivative at x(0) is given in vc(0).
340 · 2 if second derivative at x(0) is given in vc(0).
342 · 3 to use the 3-point difference formula for d(0). (Reverts to the
344 default b.c. if "n < 3")
345 · 4 to use the 4-point difference formula for d(0). (Reverts to the
347 default b.c. if "n < 4")
348 · 5 to set d(0) so that the second derivative is continuous at x(1).
350 (Reverts to the default b.c. if "n < 4")
351 The values for ic(1) are the same as above, except that the first-
353 derivative value is stored in vc(1) for cases 1 and 2. The values of
354 $vc need only be set if options 1 or 2 are chosen for $ic.
355 Set "$mflag = 0" if interpolant is required to be monotonic in each
357 interval, regardless of the data. This causes $d to be set to 0 at all
358 switch points. Set $mflag to be non-zero to use a formula based on the
359 3-point difference formula at switch points. If "$mflag > 0", then the
360 interpolant at swich points is forced to not deviate from the data by
361 more than "$mflag*dfloc", where "dfloc" is the maximum of the change of
362 $f on this interval and its two immediate neighbours. If "$mflag < 0",
363 no such control is to be imposed.
364 The piddle $wk is only needed for work space. However, I could not get
366 it to work as a temporary variable, so you must supply it; it is a 1D
367 piddle with "2*n" elements.
368 Error status returned by $ierr:
370 · 0 if successful.
372 · 1 if "ic(0) < 0" and d(0) had to be adjusted for monotonicity.
374 · 2 if "ic(1) < 0" and "d(n-1)" had to be adjusted for monotonicity.
376 · 3 if both 1 and 2 are true.
378 · -1 if "n < 2".
380 · -3 if $x is not strictly increasing.
382 · -4 if "abs(ic(0)) > 5".
384 · -5 if "abs(ic(1)) > 5".
386 · -6 if both -4 and -5 are true.
388 · -7 if "nwk < 2*(n-1)".
390 chic does not process bad values. It will set the bad-value flag of
392 all output piddles if the flag is set for any of the input piddles.
393 chsp
395 Signature: (int ic(two=2);vc(two=2);x(n);f(n);[o]d(n);wk(nwk);int [o]ierr())
396 Calculate the derivatives of (x,f(x)) using cubic spline interpolation.
398 Calculate the derivatives, using cubic spline interpolation, at the
400 given points ("$x,$f"), with the specified boundary conditions.
401 Control over the boundary conditions is given by the $ic and $vc
402 piddles. The resulting values - "$x,$f,$d" - can be used in all the
403 functions which expect a cubic Hermite function.
404 The first and second elements of $ic determine the boundary conditions
406 at the start and end of the data respectively. The allowed values for
407 ic(0) are:
408 · 0 to set d(0) so that the third derivative is continuous at x(1).
410 · 1 if first derivative at x(0) is given in "vc(0").
412 · 2 if second derivative at "x(0") is given in vc(0).
414 · 3 to use the 3-point difference formula for d(0). (Reverts to the
416 default b.c. if "n < 3".)
417 · 4 to use the 4-point difference formula for d(0). (Reverts to the
419 default b.c. if "n < 4".)
420 The values for ic(1) are the same as above, except that the first-
422 derivative value is stored in vc(1) for cases 1 and 2. The values of
423 $vc need only be set if options 1 or 2 are chosen for $ic.
424 The piddle $wk is only needed for work space. However, I could not get
426 it to work as a temporary variable, so you must supply it; it is a 1D
427 piddle with "2*n" elements.
428 Error status returned by $ierr:
430 · 0 if successful.
432 · -1 if "nelem($x) < 2".
434 · -3 if $x is not strictly increasing.
436 · -4 if "ic(0) < 0" or "ic(0) > 4".
438 · -5 if "ic(1) < 0" or "ic(1) > 4".
440 · -6 if both of the above are true.
442 · -7 if "nwk < 2*n".
444 · -8 in case of trouble solving the linear system for the interior
446 derivative values.
447 chsp does not process bad values. It will set the bad-value flag of
449 all output piddles if the flag is set for any of the input piddles.
450 chfd
452 Signature: (x(n);f(n);d(n);int check();xe(ne);[o]fe(ne);[o]de(ne);int [o]ierr())
453 Interpolate function and derivative values.
455 Given a piecewise cubic Hermite function - such as from chim - evaluate
457 the function ($fe) and derivative ($de) at a set of points ($xe). If
458 function values alone are required, use chfe. Set "check" to 0 to skip
459 checks on the input data.
460 Error status returned by $ierr:
462 · 0 if successful.
464 · >0 if extrapolation was performed at "ierr" points (data still
466 valid).
467 · -1 if "nelem($x) < 2"
469 · -3 if $x is not strictly increasing.
471 · -4 if "nelem($xe) < 1".
473 · -5 if an error has occurred in a lower-level routine, which should
475 never happen.
476 chfd does not process bad values. It will set the bad-value flag of
478 all output piddles if the flag is set for any of the input piddles.
479 chfe
481 Signature: (x(n);f(n);d(n);int check();xe(ne);[o]fe(ne);int [o]ierr())
482 Interpolate function values.
484 Given a piecewise cubic Hermite function - such as from chim - evaluate
486 the function ($fe) at a set of points ($xe). If derivative values are
487 also required, use chfd. Set "check" to 0 to skip checks on the input
488 data.
489 Error status returned by $ierr:
491 · 0 if successful.
493 · >0 if extrapolation was performed at "ierr" points (data still
495 valid).
496 · -1 if "nelem($x) < 2"
498 · -3 if $x is not strictly increasing.
500 · -4 if "nelem($xe) < 1".
502 chfe does not process bad values. It will set the bad-value flag of
504 all output piddles if the flag is set for any of the input piddles.
505 chia
507 Signature: (x(n);f(n);d(n);int check();a();b();[o]ans();int [o]ierr())
508 Integrate (x,f(x)) over arbitrary limits.
510 Evaluate the definite integral of a a piecewise cubic Hermite function
512 over an arbitrary interval, given by "[$a,$b]". See chid if the
513 integration limits are data points. Set "check" to 0 to skip checks on
514 the input data.
515 The values of $a and $b do not have to lie within $x, although the
517 resulting integral value will be highly suspect if they are not.
518 Error status returned by $ierr:
520 · 0 if successful.
522 · 1 if $a lies outside $x.
524 · 2 if $b lies outside $x.
526 · 3 if both 1 and 2 are true.
528 · -1 if "nelem($x) < 2"
530 · -3 if $x is not strictly increasing.
532 · -4 if an error has occurred in a lower-level routine, which should
534 never happen.
535 chia does not process bad values. It will set the bad-value flag of
537 all output piddles if the flag is set for any of the input piddles.
538 chid
540 Signature: (x(n);f(n);d(n);int check();int ia();int ib();[o]ans();int [o]ierr())
541 Integrate (x,f(x)) between data points.
543 Evaluate the definite integral of a a piecewise cubic Hermite function
545 between "x($ia)" and "x($ib)".
546 See chia for integration between arbitrary limits.
548 Although using a fortran routine, the values of $ia and $ib are zero
550 offset. Set "check" to 0 to skip checks on the input data.
551 Error status returned by $ierr:
553 · 0 if successful.
555 · -1 if "nelem($x) < 2".
557 · -3 if $x is not strictly increasing.
559 · -4 if $ia or $ib is out of range.
561 chid does not process bad values. It will set the bad-value flag of
563 all output piddles if the flag is set for any of the input piddles.
564 chcm
566 Signature: (x(n);f(n);d(n);int check();int [o]ismon(n);int [o]ierr())
567 Check the given piecewise cubic Hermite function for monotonicity.
569 The outout piddle $ismon indicates over which intervals the function is
571 monotonic. Set "check" to 0 to skip checks on the input data.
572 For the data interval "[x(i),x(i+1)]", the values of ismon(i) can be:
574 · -3 if function is probably decreasing
576 · -1 if function is strictly decreasing
578 · 0 if function is constant
580 · 1 if function is strictly increasing
582 · 2 if function is non-monotonic
584 · 3 if function is probably increasing
586 If "abs(ismon(i)) == 3", the derivative values are near the boundary of
588 the monotonicity region. A small increase produces non-monotonicity,
589 whereas a decrease produces strict monotonicity.
590 The above applies to "i = 0 .. nelem($x)-1". The last element of $ismon
592 indicates whether the entire function is monotonic over $x.
593 Error status returned by $ierr:
595 · 0 if successful.
597 · -1 if "n < 2".
599 · -3 if $x is not strictly increasing.
601 chcm does not process bad values. It will set the bad-value flag of
603 all output piddles if the flag is set for any of the input piddles.
604 chbs
606 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())
607 Piecewise Cubic Hermite function to B-Spline converter.
609 The resulting B-spline representation of the data (i.e. "nknots", "t",
611 "bcoeff", "ndim", and "kord") can be evaluated by "bvalu" (which is
612 currently not available).
613 Array sizes: "tsize = 2*n + 4", "bsize = 2*n", and "ndim = 2*n".
615 "knotyp" is a flag which controls the knot sequence. The knot sequence
617 "t" is normally computed from $x by putting a double knot at each "x"
618 and setting the end knot pairs according to the value of "knotyp"
619 (where "m = ndim = 2*n"):
620 · 0 - Quadruple knots at the first and last points.
622 · 1 - Replicate lengths of extreme subintervals: "t( 0 ) = t( 1 ) =
624 x(0) - (x(1)-x(0))" and "t(m+3) = t(m+2) = x(n-1) +
625 (x(n-1)-x(n-2))"
626 · 2 - Periodic placement of boundary knots: "t( 0 ) = t( 1 ) = x(0)
628 - (x(n-1)-x(n-2))" and "t(m+3) = t(m+2) = x(n) + (x(1)-x(0))"
629 · <0 - Assume the "nknots" and "t" were set in a previous call.
631 "nknots" is the number of knots and may be changed by the routine. If
633 "knotyp >= 0", "nknots" will be set to "ndim+4", otherwise it is an
634 input variable, and an error will occur if its value is not equal to
635 "ndim+4".
636 "t" is the array of "2*n+4" knots for the B-representation and may be
638 changed by the routine. If "knotyp >= 0", "t" will be changed so that
639 the interior double knots are equal to the x-values and the boundary
640 knots set as indicated above, otherwise it is assumed that "t" was set
641 by a previous call (no check is made to verify that the data forms a
642 legitimate knot sequence).
643 Error status returned by $ierr:
645 · 0 if successful.
647 · -4 if "knotyp > 2".
649 · -5 if "knotyp < 0" and "nknots != 2*n + 4".
651 chbs does not process bad values. It will set the bad-value flag of
653 all output piddles if the flag is set for any of the input piddles.
654 polfit
656 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))
657 Fit discrete data in a least squares sense by polynomials
659 in one variable. "x()", "y()" and "w()" must be of the same
660 type. This version handles bad values appropriately
661 polfit processes bad values. It will set the bad-value flag of all
663 output piddles if the flag is set for any of the input piddles.
664
666 Copyright (C) 1997 Tuomas J. Lukka. Copyright (C) 2000 Tim Jenness,
667 Doug Burke. All rights reserved. There is no warranty. You are allowed
668 to redistribute this software / documentation under certain conditions.
669 For details, see the file COPYING in the PDL distribution. If this file
670 is separated from the PDL distribution, the copyright notice should be
671 included in the file.
672perl v5.30.0 2019-09-05 Slatec(3)