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, $c) = 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 ndarrays.
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
43 both the function ("f") and derivative ("d") values at a set of
44 points ("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
50 "chcm".
51
53 eigsys
54 Eigenvalues and eigenvectors of a real positive definite symmetric
55 matrix.
56 ($eigvals,$eigvecs) = eigsys($mat)
58 Note: this function should be extended to calculate only eigenvalues if
60 called in scalar context!
61 matinv
63 Inverse of a square matrix
64 ($inv) = matinv($mat)
66 polyfit
68 Convenience wrapper routine about the "polfit" "slatec" function.
69 Separates supplied arguments and return values.
70 Fit discrete data in a least squares sense by polynomials in one
72 variable. Handles broadcasting correctly--one can pass in a 2D PDL (as
73 $y) and it will pass back a 2D PDL, the rows of which are the
74 polynomial regression results (in $r corresponding to the rows of $y.
75 ($ndeg, $r, $ierr, $c, $coeffs, $rms) = polyfit($x, $y, $w, $maxdeg, [$eps]);
77 $coeffs = polyfit($x,$y,$w,$maxdeg,[$eps]);
79 where on input:
81 $x and $y are the values to fit to a polynomial. $w are weighting
83 factors $maxdeg is the maximum degree of polynomial to use and $eps is
84 the required degree of fit.
85 and the output switches on list/scalar context.
87 In list context:
89 $ndeg is the degree of polynomial actually used $r is the values of the
91 fitted polynomial $ierr is a return status code, and $c is some working
92 array or other (preserved for historical purposes) $coeffs is the
93 polynomial coefficients of the best fit polynomial. $rms is the rms
94 error of the fit.
95 In scalar context, only $coeffs is returned.
97 Historically, $eps was modified in-place to be a return value of the
99 rms error. This usage is deprecated, and $eps is an optional parameter
100 now. It is still modified if present.
101 $c is a working array accessible to Slatec - you can feed it to several
103 other Slatec routines to get nice things out. It does not broadcast
104 correctly and should probably be fixed by someone. If you are reading
105 this, that someone might be you.
106 This version of polyfit handles bad values correctly. Bad values in $y
108 are ignored for the fit and give computed values on the fitted curve in
109 the return. Bad values in $x or $w are ignored for the fit and result
110 in bad elements in the output.
111 polycoef
113 Convenience wrapper routine around the "pcoef" "slatec" function.
114 Separates supplied arguments and return values.
115 Convert the "polyfit"/"polfit" coefficients to Taylor series form.
117 $tc = polycoef($l, $c, $x);
119 polyvalue
121 Convenience wrapper routine around the "pvalue" "slatec" function.
122 Separates supplied arguments and return values.
123 For multiple input x positions, a corresponding y position is
125 calculated.
126 The derivatives PDL is one dimensional (of size "nder") if a single x
128 position is supplied, two dimensional if more than one x position is
129 supplied.
130 Use the coefficients "c" generated by "polyfit" (or "polfit") to
132 evaluate the polynomial fit of degree "l", along with the first "nder"
133 of its derivatives, at a specified point "x".
134 ($yfit, $yp) = polyvalue($l, $nder, $x, $c);
136 detslatec
138 compute the determinant of an invertible matrix
139 $mat = zeroes(5,5); $mat->diagonal(0,1) .= 1; # unity matrix
141 $det = detslatec $mat;
142 Usage:
144 $determinant = detslatec $matrix;
146 Signature: detslatec(mat(n,m); [o] det())
148 "detslatec" computes the determinant of an invertible matrix and barfs
150 if the matrix argument provided is non-invertible. The matrix
151 broadcasts as usual.
152 This routine was previously known as "det" which clashes now with det
154 which is provided by PDL::MatrixOps.
155 fft
157 Fast Fourier Transform
158 $v_in = pdl(1,0,1,0);
160 ($azero,$x,$y) = PDL::Slatec::fft($v_in);
161 "PDL::Slatec::fft" is a convenience wrapper for "ezfftf", and performs
163 a Fast Fourier Transform on an input vector $v_in. The return values
164 are the same as for "ezfftf".
165 rfft
167 reverse Fast Fourier Transform
168 $v_out = PDL::Slatec::rfft($azero,$x,$y);
170 print $v_in, $vout
171 [1 0 1 0] [1 0 1 0]
172 "PDL::Slatec::rfft" is a convenience wrapper for "ezfftb", and performs
174 a reverse Fast Fourier Transform. The input is the same as the output
175 of "PDL::Slatec::fft", and the output of "rfft" is a data vector,
176 similar to what is input into "PDL::Slatec::fft".
177 svdc
179 Signature: (x(n,p);[o]s(p);[o]e(p);[o]u(n,p);[o]v(p,p);[o]work(n);longlong job();longlong [o]info())
180 singular value decomposition of a matrix
182 svdc does not process bad values. It will set the bad-value flag of
184 all output ndarrays if the flag is set for any of the input ndarrays.
185 poco
187 Signature: (a(n,n);rcond();[o]z(n);longlong [o]info())
188 Factor a real symmetric positive definite matrix and estimate the
190 condition number of the matrix.
191 poco does not process bad values. It will set the bad-value flag of
193 all output ndarrays if the flag is set for any of the input ndarrays.
194 geco
196 Signature: (a(n,n);longlong [o]ipvt(n);[o]rcond();[o]z(n))
197 Factor a matrix using Gaussian elimination and estimate the condition
199 number of the matrix.
200 geco does not process bad values. It will set the bad-value flag of
202 all output ndarrays if the flag is set for any of the input ndarrays.
203 gefa
205 Signature: (a(n,n);longlong [o]ipvt(n);longlong [o]info())
206 Factor a matrix using Gaussian elimination.
208 gefa does not process bad values. It will set the bad-value flag of
210 all output ndarrays if the flag is set for any of the input ndarrays.
211 podi
213 Signature: (a(n,n);[o]det(two=2);longlong job())
214 Compute the determinant and inverse of a certain real symmetric
216 positive definite matrix using the factors computed by "poco".
217 podi does not process bad values. It will set the bad-value flag of
219 all output ndarrays if the flag is set for any of the input ndarrays.
220 gedi
222 Signature: (a(n,n);longlong [o]ipvt(n);[o]det(two=2);[o]work(n);longlong job())
223 Compute the determinant and inverse of a matrix using the factors
225 computed by "geco" or "gefa".
226 gedi does not process bad values. It will set the bad-value flag of
228 all output ndarrays if the flag is set for any of the input ndarrays.
229 gesl
231 Signature: (a(lda,n);longlong ipvt(n);b(n);longlong job())
232 Solve the real system "A*X=B" or "TRANS(A)*X=B" using the factors
234 computed by "geco" or "gefa".
235 gesl does not process bad values. It will set the bad-value flag of
237 all output ndarrays if the flag is set for any of the input ndarrays.
238 rs
240 Signature: (a(n,n);[o]w(n);longlong matz();[o]z(n,n);[t]fvone(n);[t]fvtwo(n);longlong [o]ierr())
241 This subroutine calls the recommended sequence of subroutines from the
243 eigensystem subroutine package (EISPACK) to find the eigenvalues and
244 eigenvectors (if desired) of a REAL SYMMETRIC matrix.
245 rs does not process bad values. It will set the bad-value flag of all
247 output ndarrays if the flag is set for any of the input ndarrays.
248 ezffti
250 Signature: (longlong n();[o]wsave(foo))
251 Subroutine ezffti initializes the work array "wsave()" which is used in
253 both "ezfftf" and "ezfftb". The prime factorization of "n" together
254 with a tabulation of the trigonometric functions are computed and
255 stored in "wsave()".
256 ezffti does not process bad values. It will set the bad-value flag of
258 all output ndarrays if the flag is set for any of the input ndarrays.
259 ezfftf
261 Signature: (r(n);[o]azero();[o]a(n);[o]b(n);wsave(foo))
262 ezfftf does not process bad values. It will set the bad-value flag of
264 all output ndarrays if the flag is set for any of the input ndarrays.
265 ezfftb
267 Signature: ([o]r(n);azero();a(n);b(n);wsave(foo))
268 ezfftb does not process bad values. It will set the bad-value flag of
270 all output ndarrays if the flag is set for any of the input ndarrays.
271 pcoef
273 Signature: (longlong l();c();[o]tc(bar);a(foo))
274 Convert the "polfit" coefficients to Taylor series form. "c" and "a()"
276 must be of the same type.
277 pcoef does not process bad values. It will set the bad-value flag of
279 all output ndarrays if the flag is set for any of the input ndarrays.
280 pvalue
282 Signature: (longlong l();x();[o]yfit();[o]yp(nder);a(foo))
283 Use the coefficients generated by "polfit" to evaluate the polynomial
285 fit of degree "l", along with the first "nder" of its derivatives, at a
286 specified point. "x" and "a" must be of the same type.
287 pvalue does not process bad values. It will set the bad-value flag of
289 all output ndarrays if the flag is set for any of the input ndarrays.
290 chim
292 Signature: (x(n);f(n);[o]d(n);longlong [o]ierr())
293 Calculate the derivatives of (x,f(x)) using cubic Hermite
295 interpolation.
296 Calculate the derivatives at the given set of points ("$x,$f", where $x
298 is strictly increasing). The resulting set of points - "$x,$f,$d",
299 referred to as the cubic Hermite representation - can then be used in
300 other functions, such as "chfe", "chfd", and "chia".
301 The boundary conditions are compatible with monotonicity, and if the
303 data are only piecewise monotonic, the interpolant will have an
304 extremum at the switch points; for more control over these issues use
305 "chic".
306 Error status returned by $ierr:
308 • 0 if successful.
310 • > 0 if there were "ierr" switches in the direction of monotonicity
312 (data still valid).
313 • -1 if "nelem($x) < 2".
315 • -3 if $x is not strictly increasing.
317 chim does not process bad values. It will set the bad-value flag of
319 all output ndarrays if the flag is set for any of the input ndarrays.
320 chic
322 Signature: (longlong ic(two=2);vc(two=2);mflag();x(n);f(n);[o]d(n);wk(nwk);longlong [o]ierr())
323 Calculate the derivatives of (x,f(x)) using cubic Hermite
325 interpolation.
326 Calculate the derivatives at the given points ("$x,$f", where $x is
328 strictly increasing). Control over the boundary conditions is given by
329 the $ic and $vc ndarrays, and the value of $mflag determines the
330 treatment of points where monotoncity switches direction. A simpler,
331 more restricted, interface is available using "chim".
332 The first and second elements of $ic determine the boundary conditions
334 at the start and end of the data respectively. If the value is 0, then
335 the default condition, as used by "chim", is adopted. If greater than
336 zero, no adjustment for monotonicity is made, otherwise if less than
337 zero the derivative will be adjusted. The allowed magnitudes for ic(0)
338 are:
339 • 1 if first derivative at x(0) is given in vc(0).
341 • 2 if second derivative at x(0) is given in vc(0).
343 • 3 to use the 3-point difference formula for d(0). (Reverts to the
345 default b.c. if "n < 3")
346 • 4 to use the 4-point difference formula for d(0). (Reverts to the
348 default b.c. if "n < 4")
349 • 5 to set d(0) so that the second derivative is continuous at x(1).
351 (Reverts to the default b.c. if "n < 4")
352 The values for ic(1) are the same as above, except that the first-
354 derivative value is stored in vc(1) for cases 1 and 2. The values of
355 $vc need only be set if options 1 or 2 are chosen for $ic.
356 Set "$mflag = 0" if interpolant is required to be monotonic in each
358 interval, regardless of the data. This causes $d to be set to 0 at all
359 switch points. Set $mflag to be non-zero to use a formula based on the
360 3-point difference formula at switch points. If "$mflag > 0", then the
361 interpolant at swich points is forced to not deviate from the data by
362 more than "$mflag*dfloc", where "dfloc" is the maximum of the change of
363 $f on this interval and its two immediate neighbours. If "$mflag < 0",
364 no such control is to be imposed.
365 The ndarray $wk is only needed for work space. However, I could not get
367 it to work as a temporary variable, so you must supply it; it is a 1D
368 ndarray with "2*n" elements.
369 Error status returned by $ierr:
371 • 0 if successful.
373 • 1 if "ic(0) < 0" and d(0) had to be adjusted for monotonicity.
375 • 2 if "ic(1) < 0" and "d(n-1)" had to be adjusted for monotonicity.
377 • 3 if both 1 and 2 are true.
379 • -1 if "n < 2".
381 • -3 if $x is not strictly increasing.
383 • -4 if "abs(ic(0)) > 5".
385 • -5 if "abs(ic(1)) > 5".
387 • -6 if both -4 and -5 are true.
389 • -7 if "nwk < 2*(n-1)".
391 chic does not process bad values. It will set the bad-value flag of
393 all output ndarrays if the flag is set for any of the input ndarrays.
394 chsp
396 Signature: (longlong ic(two=2);vc(two=2);x(n);f(n);[o]d(n);wk(nwk);longlong [o]ierr())
397 Calculate the derivatives of (x,f(x)) using cubic spline interpolation.
399 Calculate the derivatives, using cubic spline interpolation, at the
401 given points ("$x,$f"), with the specified boundary conditions.
402 Control over the boundary conditions is given by the $ic and $vc
403 ndarrays. The resulting values - "$x,$f,$d" - can be used in all the
404 functions which expect a cubic Hermite function.
405 The first and second elements of $ic determine the boundary conditions
407 at the start and end of the data respectively. The allowed values for
408 ic(0) are:
409 • 0 to set d(0) so that the third derivative is continuous at x(1).
411 • 1 if first derivative at x(0) is given in "vc(0").
413 • 2 if second derivative at "x(0") is given in vc(0).
415 • 3 to use the 3-point difference formula for d(0). (Reverts to the
417 default b.c. if "n < 3".)
418 • 4 to use the 4-point difference formula for d(0). (Reverts to the
420 default b.c. if "n < 4".)
421 The values for ic(1) are the same as above, except that the first-
423 derivative value is stored in vc(1) for cases 1 and 2. The values of
424 $vc need only be set if options 1 or 2 are chosen for $ic.
425 The ndarray $wk is only needed for work space. However, I could not get
427 it to work as a temporary variable, so you must supply it; it is a 1D
428 ndarray with "2*n" elements.
429 Error status returned by $ierr:
431 • 0 if successful.
433 • -1 if "nelem($x) < 2".
435 • -3 if $x is not strictly increasing.
437 • -4 if "ic(0) < 0" or "ic(0) > 4".
439 • -5 if "ic(1) < 0" or "ic(1) > 4".
441 • -6 if both of the above are true.
443 • -7 if "nwk < 2*n".
445 • -8 in case of trouble solving the linear system for the interior
447 derivative values.
448 chsp does not process bad values. It will set the bad-value flag of
450 all output ndarrays if the flag is set for any of the input ndarrays.
451 chfd
453 Signature: (x(n);f(n);d(n);longlong check();xe(ne);[o]fe(ne);[o]de(ne);longlong [o]ierr())
454 Interpolate function and derivative values.
456 Given a piecewise cubic Hermite function - such as from "chim" -
458 evaluate the function ($fe) and derivative ($de) at a set of points
459 ($xe). If function values alone are required, use "chfe". Set "check"
460 to 0 to skip checks on the input data.
461 Error status returned by $ierr:
463 • 0 if successful.
465 • >0 if extrapolation was performed at "ierr" points (data still
467 valid).
468 • -1 if "nelem($x) < 2"
470 • -3 if $x is not strictly increasing.
472 • -4 if "nelem($xe) < 1".
474 • -5 if an error has occurred in a lower-level routine, which should
476 never happen.
477 chfd does not process bad values. It will set the bad-value flag of
479 all output ndarrays if the flag is set for any of the input ndarrays.
480 chfe
482 Signature: (x(n);f(n);d(n);longlong check();xe(ne);[o]fe(ne);longlong [o]ierr())
483 Interpolate function values.
485 Given a piecewise cubic Hermite function - such as from "chim" -
487 evaluate the function ($fe) at a set of points ($xe). If derivative
488 values are also required, use "chfd". Set "check" to 0 to skip checks
489 on the input data.
490 Error status returned by $ierr:
492 • 0 if successful.
494 • >0 if extrapolation was performed at "ierr" points (data still
496 valid).
497 • -1 if "nelem($x) < 2"
499 • -3 if $x is not strictly increasing.
501 • -4 if "nelem($xe) < 1".
503 chfe does not process bad values. It will set the bad-value flag of
505 all output ndarrays if the flag is set for any of the input ndarrays.
506 chia
508 Signature: (x(n);f(n);d(n);longlong check();la();lb();[o]ans();longlong [o]ierr())
509 Integrate (x,f(x)) over arbitrary limits.
511 Evaluate the definite integral of a piecewise cubic Hermite function
513 over an arbitrary interval, given by "[$la,$lb]". $d should contain the
514 derivative values, computed by "chim". See "chid" if the integration
515 limits are data points. Set "check" to 0 to skip checks on the input
516 data.
517 The values of $la and $lb do not have to lie within $x, although the
519 resulting integral value will be highly suspect if they are not.
520 Error status returned by $ierr:
522 • 0 if successful.
524 • 1 if $la lies outside $x.
526 • 2 if $lb lies outside $x.
528 • 3 if both 1 and 2 are true.
530 • -1 if "nelem($x) < 2"
532 • -3 if $x is not strictly increasing.
534 • -4 if an error has occurred in a lower-level routine, which should
536 never happen.
537 chia does not process bad values. It will set the bad-value flag of
539 all output ndarrays if the flag is set for any of the input ndarrays.
540 chid
542 Signature: (x(n);f(n);d(n);longlong check();longlong ia();longlong ib();[o]ans();longlong [o]ierr())
543 Integrate (x,f(x)) between data points.
545 Evaluate the definite integral of a a piecewise cubic Hermite function
547 between "x($ia)" and "x($ib)".
548 See "chia" for integration between arbitrary limits.
550 Although using a fortran routine, the values of $ia and $ib are zero
552 offset. $d should contain the derivative values, computed by "chim".
553 Set "check" to 0 to skip checks on the input data.
554 Error status returned by $ierr:
556 • 0 if successful.
558 • -1 if "nelem($x) < 2".
560 • -3 if $x is not strictly increasing.
562 • -4 if $ia or $ib is out of range.
564 chid does not process bad values. It will set the bad-value flag of
566 all output ndarrays if the flag is set for any of the input ndarrays.
567 chcm
569 Signature: (x(n);f(n);d(n);longlong check();longlong [o]ismon(n);longlong [o]ierr())
570 Check the given piecewise cubic Hermite function for monotonicity.
572 The outout ndarray $ismon indicates over which intervals the function
574 is monotonic. Set "check" to 0 to skip checks on the input data.
575 For the data interval "[x(i),x(i+1)]", the values of ismon(i) can be:
577 • -3 if function is probably decreasing
579 • -1 if function is strictly decreasing
581 • 0 if function is constant
583 • 1 if function is strictly increasing
585 • 2 if function is non-monotonic
587 • 3 if function is probably increasing
589 If "abs(ismon(i)) == 3", the derivative values are near the boundary of
591 the monotonicity region. A small increase produces non-monotonicity,
592 whereas a decrease produces strict monotonicity.
593 The above applies to "i = 0 .. nelem($x)-1". The last element of $ismon
595 indicates whether the entire function is monotonic over $x.
596 Error status returned by $ierr:
598 • 0 if successful.
600 • -1 if "n < 2".
602 • -3 if $x is not strictly increasing.
604 chcm does not process bad values. It will set the bad-value flag of
606 all output ndarrays if the flag is set for any of the input ndarrays.
607 chbs
609 Signature: (x(n);f(n);d(n);longlong knotyp();longlong nknots();t(tsize);[o]bcoef(bsize);longlong [o]ndim();longlong [o]kord();longlong [o]ierr())
610 Piecewise Cubic Hermite function to B-Spline converter.
612 The resulting B-spline representation of the data (i.e. "nknots", "t",
614 "bcoeff", "ndim", and "kord") can be evaluated by "bvalu" (which is
615 currently not available).
616 Array sizes: "tsize = 2*n + 4", "bsize = 2*n", and "ndim = 2*n".
618 "knotyp" is a flag which controls the knot sequence. The knot sequence
620 "t" is normally computed from $x by putting a double knot at each "x"
621 and setting the end knot pairs according to the value of "knotyp"
622 (where "m = ndim = 2*n"):
623 • 0 - Quadruple knots at the first and last points.
625 • 1 - Replicate lengths of extreme subintervals: "t( 0 ) = t( 1 ) =
627 x(0) - (x(1)-x(0))" and "t(m+3) = t(m+2) = x(n-1) +
628 (x(n-1)-x(n-2))"
629 • 2 - Periodic placement of boundary knots: "t( 0 ) = t( 1 ) = x(0)
631 - (x(n-1)-x(n-2))" and "t(m+3) = t(m+2) = x(n) + (x(1)-x(0))"
632 • <0 - Assume the "nknots" and "t" were set in a previous call.
634 "nknots" is the number of knots and may be changed by the routine. If
636 "knotyp >= 0", "nknots" will be set to "ndim+4", otherwise it is an
637 input variable, and an error will occur if its value is not equal to
638 "ndim+4".
639 "t" is the array of "2*n+4" knots for the B-representation and may be
641 changed by the routine. If "knotyp >= 0", "t" will be changed so that
642 the interior double knots are equal to the x-values and the boundary
643 knots set as indicated above, otherwise it is assumed that "t" was set
644 by a previous call (no check is made to verify that the data forms a
645 legitimate knot sequence).
646 Error status returned by $ierr:
648 • 0 if successful.
650 • -4 if "knotyp > 2".
652 • -5 if "knotyp < 0" and "nknots != 2*n + 4".
654 chbs does not process bad values. It will set the bad-value flag of
656 all output ndarrays if the flag is set for any of the input ndarrays.
657 polfit
659 Signature: (x(n); y(n); w(n); longlong maxdeg(); longlong [o]ndeg(); [o]eps(); [o]r(n); longlong [o]ierr(); [o]a(foo); [o]coeffs(bar);[t]xtmp(n);[t]ytmp(n);[t]wtmp(n);[t]rtmp(n))
660 Fit discrete data in a least squares sense by polynomials
662 in one variable. "x()", "y()" and "w()" must be of the same
663 type. This version handles bad values appropriately
664 polfit processes bad values. It will set the bad-value flag of all
666 output ndarrays if the flag is set for any of the input ndarrays.
667
669 Copyright (C) 1997 Tuomas J. Lukka. Copyright (C) 2000 Tim Jenness,
670 Doug Burke. All rights reserved. There is no warranty. You are allowed
671 to redistribute this software / documentation under certain conditions.
672 For details, see the file COPYING in the PDL distribution. If this file
673 is separated from the PDL distribution, the copyright notice should be
674 included in the file.
675perl v5.34.0 2022-02-28 Slatec(3)