1pp_addpm({At=>Top},<<'EOD'); 2 3=head1 NAME 4 5PDL::Slatec - PDL interface to the slatec numerical programming library 6 7=head1 SYNOPSIS 8 9 use PDL::Slatec; 10 11 ($ndeg, $r, $ierr, $a) = polyfit($x, $y, $w, $maxdeg, $eps); 12 13=head1 DESCRIPTION 14 15This module serves the dual purpose of providing an interface to 16parts of the slatec library and showing how to interface PDL 17to an external library. 18Using this library requires a fortran compiler; the source for the routines 19is provided for convenience. 20 21Currently available are routines to: 22manipulate matrices; calculate FFT's; 23fit data using polynomials; 24and interpolate/integrate data using piecewise cubic Hermite interpolation. 25 26=head2 Piecewise cubic Hermite interpolation (PCHIP) 27 28PCHIP is the slatec package of routines to perform piecewise cubic 29Hermite interpolation of data. 30It features software to produce a monotone and "visually pleasing" 31interpolant to monotone data. 32According to Fritsch & Carlson ("Monotone piecewise 33cubic interpolation", SIAM Journal on Numerical Analysis 3417, 2 (April 1980), pp. 238-246), 35such an interpolant may be more reasonable than a cubic spline if 36the data contains both "steep" and "flat" sections. 37Interpolation of cumulative probability distribution functions is 38another application. 39These routines are cryptically named (blame FORTRAN), 40beginning with 'ch', and accept either float or double piddles. 41 42Most of the routines require an integer parameter called C<check>; 43if set to 0, then no checks on the validity of the input data are 44made, otherwise these checks are made. 45The value of C<check> can be set to 0 if a routine 46such as L<chim|/chim> has already been successfully called. 47 48=over 4 49 50=item * 51 52If not known, estimate derivative values for the points 53using the L<chim|/chim>, L<chic|/chic>, or L<chsp|/chsp> routines 54(the following routines require both the function (C<f>) 55and derivative (C<d>) values at a set of points (C<x>)). 56 57=item * 58 59Evaluate, integrate, or differentiate the resulting PCH 60function using the routines: 61L<chfd|/chfd>; L<chfe|/chfe>; L<chia|/chia>; L<chid|/chid>. 62 63=item * 64 65If desired, you can check the monotonicity of your 66data using L<chcm|/chcm>. 67 68=back 69 70=cut 71 72EOD 73# ' un-confuse emacs 74 75# if define chbs, then add something like the following to point 3: 76# 77# or use L<chbs|/chbs> to convert a PCH function into B-representation 78# for use with the B-spline routines of slatec 79# (although no interface to them currently exist). 80# 81 82# add function definitions after finishing the first pp_addpm(), since this 83# adds a '=head1 FUNCTIONS' line at the end of the text 84 85pp_addpm(<<'END'); 86=head2 eigsys 87 88=for ref 89 90Eigenvalues and eigenvectors of a real positive definite symmetric matrix. 91 92=for usage 93 94 ($eigvals,$eigvecs) = eigsys($mat) 95 96Note: this function should be extended to calculate only eigenvalues if called 97in scalar context! 98 99=head2 matinv 100 101=for ref 102 103Inverse of a square matrix 104 105=for usage 106 107 ($inv) = matinv($mat) 108 109=head2 polyfit 110 111Convenience wrapper routine about the C<polfit> C<slatec> function. 112Separates supplied arguments and return values. 113 114=for ref 115 116Fit discrete data in a least squares sense by polynomials 117in one variable. Handles threading correctly--one can pass in a 2D PDL (as C<$y>) 118and it will pass back a 2D PDL, the rows of which are the polynomial regression 119results (in C<$r> corresponding to the rows of $y. 120 121=for usage 122 123 ($ndeg, $r, $ierr, $a, $coeffs, $rms) = polyfit($x, $y, $w, $maxdeg, [$eps]); 124 125 $coeffs = polyfit($x,$y,$w,$maxdeg,[$eps]); 126 127where on input: 128 129C<$x> and C<$y> are the values to fit to a polynomial. 130C<$w> are weighting factors 131C<$maxdeg> is the maximum degree of polynomial to use and 132C<$eps> is the required degree of fit. 133 134and the output switches on list/scalar context. 135 136In list context: 137 138C<$ndeg> is the degree of polynomial actually used 139C<$r> is the values of the fitted polynomial 140C<$ierr> is a return status code, and 141C<$a> is some working array or other (preserved for historical purposes) 142C<$coeffs> is the polynomial coefficients of the best fit polynomial. 143C<$rms> is the rms error of the fit. 144 145In scalar context, only $coeffs is returned. 146 147Historically, C<$eps> was modified in-place to be a return value of the 148rms error. This usage is deprecated, and C<$eps> is an optional parameter now. 149It is still modified if present. 150 151C<$a> is a working array accessible to Slatec - you can feed it to several 152other Slatec routines to get nice things out. It does not thread 153correctly and should probably be fixed by someone. If you are 154reading this, that someone might be you. 155 156=for bad 157 158This version of polyfit handles bad values correctly. Bad values in 159$y are ignored for the fit and give computed values on the fitted 160curve in the return. Bad values in $x or $w are ignored for the fit and 161result in bad elements in the output. 162 163=head2 polycoef 164 165Convenience wrapper routine around the C<pcoef> C<slatec> function. 166Separates supplied arguments and return values. 167 168=for ref 169 170Convert the C<polyfit>/C<polfit> coefficients to Taylor series form. 171 172=for usage 173 174 $tc = polycoef($l, $c, $a); 175 176=head2 polyvalue 177 178Convenience wrapper routine around the C<pvalue> C<slatec> function. 179Separates supplied arguments and return values. 180 181For multiple input x positions, a corresponding y position is calculated. 182 183The derivatives PDL is one dimensional (of size C<nder>) if a single x 184position is supplied, two dimensional if more than one x position is 185supplied. 186 187=for ref 188 189Use the coefficients generated by C<polyfit> (or C<polfit>) to evaluate 190the polynomial fit of degree C<l>, along with the first C<nder> of its 191derivatives, at a specified point. 192 193=for usage 194 195 ($yfit, $yp) = polyvalue($l, $nder, $x, $a); 196 197=head2 detslatec 198 199=for ref 200 201compute the determinant of an invertible matrix 202 203=for example 204 205 $mat = zeroes(5,5); $mat->diagonal(0,1) .= 1; # unity matrix 206 $det = detslatec $mat; 207 208Usage: 209 210=for usage 211 212 $determinant = detslatec $matrix; 213 214=for sig 215 216 Signature: detslatec(mat(n,m); [o] det()) 217 218C<detslatec> computes the determinant of an invertible matrix and barfs if 219the matrix argument provided is non-invertible. The matrix threads as usual. 220 221This routine was previously known as C<det> which clashes now with 222L<det|PDL::MatrixOps/det> which is provided by L<PDL::MatrixOps>. 223 224=head2 fft 225 226=for ref 227 228Fast Fourier Transform 229 230=for example 231 232 $v_in = pdl(1,0,1,0); 233 ($azero,$a,$b) = PDL::Slatec::fft($v_in); 234 235C<PDL::Slatec::fft> is a convenience wrapper for L<ezfftf|ezfftf>, and 236performs a Fast Fourier Transform on an input vector C<$v_in>. The 237return values are the same as for L<ezfftf|ezfftf>. 238 239=head2 rfft 240 241=for ref 242 243reverse Fast Fourier Transform 244 245=for example 246 247 $v_out = PDL::Slatec::rfft($azero,$a,$b); 248 print $v_in, $vout 249 [1 0 1 0] [1 0 1 0] 250 251C<PDL::Slatec::rfft> is a convenience wrapper for L<ezfftb|ezfftb>, 252and performs a reverse Fast Fourier Transform. The input is the same 253as the output of L<PDL::Slatec::fft|/PDL::Slatec::fft>, and the output 254of C<rfft> is a data vector, similar to what is input into 255L<PDL::Slatec::fft|/PDL::Slatec::fft>. 256 257=cut 258 259END 260 261use strict; 262 263# for MDim, ld[str] is interpreted as "leading dimension of ..." 264 265# Making the BAD BAD BAD assumption that PDL_Long == int 266# in fortran. BAD BAD BAD XXX (I'm going to regret this) 267 268my %ftypes = (S => 'F', D => 'D'); 269 270sub firstpar { 271 $_[0] =~ /^\(([^),]+)[),]/ or die "Can't find first par from $_[0]"; 272 $1 273} 274 275# whether or not to append undercores 276 277my $uscore = (-e "f77_underscore" ? "_" : ""); 278 279# used in defslatec() 280#my %ignore_ppar = ( Incfd => 1, CheckFlag => 1 ); 281my %ignore_ppar = ( Incfd => 1 ); 282my %prototype = ( F => "float", D => "double" ); 283 284# an alternative is to declare the function in the Code section 285# of pp_def(), using something like: 286# 287# my $codeproto = "\$T".(join '',map {$_->[0]} @talts)."(". 288# (join ',',map {$_->[1].$uscore} @talts).") ();"; 289# if ( defined $fpar ) { 290# $codeproto = "\$T".(join '',map {$_->[0]} @talts)."(float,double) $codeproto"; 291# } 292# $codeproto = "extern $codeproto"; 293# 294# and then add `$codeproto . "\n" .' to the beginning of the Code 295# section. 296# 297# this then gets rid of the need of the prototype file. 298# 299open( PROTOS, "> SlatecProtos.h" ); 300 301# defslatec( $pdlname, $funcnames, $argstr, $docstring, $funcref ) 302# 303# $pdlname is the name of the PDL function to be created 304# $funcnames is a reference to a hash array, whose keys define 305# the single (S), and double precision (D) names of the 306# SLATEC routines to be linked to. 307# 308# $argstr is a list of arguments expected by the SLATEC routine 309# - some of the allowed type names are: 310# FuncRet 311# - specifies that this is a function, not a subroutine, and 312# that the output of the function should be stored in this 313# variable 314# Incfd 315# - used in the PCHIP functions to specify the INCFD argument 316# that we force to be 1, so the user never has to specify it 317# (this allows the PCHIP routines to use 2D data, but as it's 318# done in FORTRAN array order, and PDL has a much richer way 319# of accessing parts of an array we force the data to be 1D). 320# CheckFlag 321# - the PCHIP routined may change the value from 0 to 1 if an 322# error occurs but the checks were successful. As this complicates 323# things we copy the user value to a temporary variable, 324# so that the sent in value is not changed. 325# FortranIndex 326# - pchid()/dpchid() require FORTRAN array indices, so 327# this type flags that we should add 1 onto the input values 328# before sending to the slatec function 329# 330# $docstring gives the text to be used as the function dicumentation 331# 332# $funcref gets placed within a '=for ref' pod statement at the 333# start of the documentation - ie it is placed before the 334# text within $docstring. This string gets printed out 335# in the perldl or pdl2 shell after a '?? string' command 336# 337sub defslatec { 338 339 my $debug = 0; # print out calls to pp_def 340 341 my($pname,$fnames,$argstr,$docstring,$funcref) = @_; 342 my @args = map {/^\s*$/ ? () : $_} split ';', $argstr; 343 my @args2 = map { 344 /^\s*([a-zA-Z]+)\s+ # "Type name" 345 ((?:\[[^]]*\])?)\s* # Options 346 ([a-zA-Z]+)\s* # Par name 347 ((?:\([^)]*\))?)\s*$ # Dims 348 /x or die("Invalid slatec par $_"); 349 [$1,$2,$3,$4]} @args; 350 351 # is this for a function (Type name eq "FuncRet") 352 # or a subroutine? 353 my $fpar; 354 foreach ( @args2 ) { 355 next unless $_->[0] eq "FuncRet"; 356 die "Only one FuncRet allowed in pars list.\n" if defined $fpar; 357 $fpar = "\$$_->[2]()"; 358 } 359 360 my @ppars = map { 361 if($_->[0] =~ /^M?Dim$/ or defined $ignore_ppar{$_->[0]} ) { 362 () 363 } else { 364 (($_->[0] eq "Mat" or $_->[0] eq "FuncRet") 365 and join '',@{$_}[1,2,3]) or 366 (($_->[0] eq "IntFlag" or $_->[0] eq "FortranIndex" or $_->[0] eq "CheckFlag") 367 and "int ".join '',@{$_}[1,2,3]) or 368 die "Invalid ppars ",(join ',',@$_),"\n"; 369 } 370 } @args2; 371 372 # uncomment the following line to see what perl thinks the input pars are 373 ##print "Pars: ",(join ';',@ppars),"\n"; 374 375 my @talts = map { 376 defined $ftypes{$_} or die "FTYPE $_ NOT THERE\n"; 377 [$ftypes{$_},$fnames->{$_}] 378 } sort keys %$fnames; 379 380 my $func = "\$T".(join '',map {$_->[0]} @talts) . "(" . 381 (join ',',map {$_->[1].$uscore} @talts).")"; 382 if ( defined $fpar ) { $func = "$fpar = $func"; } 383 384 my %lds = map { 385 ($_->[0] eq "Mat" and $_->[3] ne "()") ? 386 ("ld".$_->[2] => "&\$PRIV(__".firstpar($_->[3])."_size)") 387 : () 388 } @args2; 389 390 my @funcpars; 391 foreach ( @args2 ) { 392 next if $_->[0] eq "FuncRet"; 393 if ( $_->[0] eq "Mat" or $_->[0] eq "IntFlag" ) { 394 push @funcpars, "\$P($_->[2])"; 395 } elsif ( $_->[0] eq "Dim" ) { 396 push @funcpars, "&\$PRIV(__$_->[2]_size)"; 397 } elsif ( $_->[0] eq "MDim" ) { 398 push @funcpars, $lds{$_->[2]}; 399 } elsif ( $_->[0] eq "Incfd" or $_->[0] eq "CheckFlag" ) { 400 push @funcpars, "&_" . lc($_->[0]); 401 } elsif ( $_->[0] eq "FortranIndex" ) { 402 push @funcpars, "&_$_->[2]"; 403 } else { 404 die "Invalid args2"; 405 } 406 } 407 408 # _incfd = 1 makes sure PCHIP code treats piddle as 1D 409 # _checkflag - copy input data to a temporary variable, in case 410 # the PCHIP routine decides to change it 411 # 412 my @ifincode; 413 foreach ( @args2 ) { 414 if ( $_->[0] eq "Incfd" ) { 415 push @ifincode, "int _" . lc($_->[0]) . " = 1;"; 416 } elsif ( $_->[0] eq "CheckFlag" ) { 417 push @ifincode, "int _" . lc($_->[0]) . " = \$$_->[2]();"; 418 } elsif ( $_->[0] eq "FortranIndex" ) { 419 # convert from C to F77 index 420 push @ifincode, "int _$_->[2] = \$$_->[2]() + 1;" 421 } 422 } 423 424 foreach ( @talts ) { 425 my $codeproto = "extern "; 426 if ( defined $fpar ) { $codeproto .= "$prototype{$_->[0]} "; } 427 else { $codeproto .= "int "; } 428 $codeproto .= "$_->[1]$uscore ();"; 429 print PROTOS $codeproto . "\n"; 430 } 431 432 # add on the function reference, if supplied, to the start of 433 # the doc string 434 if ( defined $docstring ) { 435 $docstring = "\n=for ref\n\n$funcref\n\n$docstring" if defined $funcref; 436 } else { 437 $docstring = ''; 438 } 439 440 # If debug flag set, then print out pp_def call for each call to defslatec 441 if ($debug) { 442 my $pars = (join ';',@ppars); 443 my $code = (join '',@ifincode) . "\n " . $func . " (". (join ',',@funcpars) . ");\n"; 444 my $generictypes = "[" . join (", ", map {$_->[0]} @talts) . "],\n"; 445 print <<"ENDDBG"; 446pp_def($pname, 447 Pars => $pars, 448 OtherPars => '', 449 Code => $code, 450 GenericTypes => $generictypes, 451 Doc => $docstring 452); 453ENDDBG 454} 455 456 pp_def($pname, 457 Pars => (join ';',@ppars), 458 OtherPars => '', 459 Code => (join '',@ifincode) . "\n " . 460 $func . " (". (join ',',@funcpars) . ");\n", 461# . (join '',@ifoutcode), 462 GenericTypes => [map {$_->[0]} @talts], 463 Doc => $docstring 464# %$opts, 465 ); 466} # sub: defslatec() 467 468pp_addhdr(qq| 469#include "SlatecProtos.h" 470 471void MAIN__ () { 472 /* Cheat to define MAIN__ symbol */ 473 croak("This should never happen"); 474} 475 476void slatecbarf$uscore() { 477 croak("slatec called halt"); 478} 479 480|); 481 482pp_add_exported('',"eigsys matinv polyfit polycoef polyvalue"); 483 484pp_addpm(<<'END'); 485 486use PDL::Core; 487use PDL::Basic; 488use PDL::Primitive; 489use PDL::Ufunc; 490use strict; 491 492# Note: handles only real symmetric positive-definite. 493 494*eigsys = \&PDL::eigsys; 495 496sub PDL::eigsys { 497 my($h) = @_; 498 $h = float($h); 499 rs($h, 500 (my $eigval=PDL->null), 501 (long (pdl (1))),(my $eigmat=PDL->null), 502 (my $fvone = PDL->null),(my $fvtwo = PDL->null), 503 (my $errflag=PDL->null) 504 ); 505# print $covar,$eigval,$eigmat,$fvone,$fvtwo,$errflag; 506 if(sum($errflag) > 0) { 507 barf("Non-positive-definite matrix given to eigsys: $h\n"); 508 } 509 return ($eigval,$eigmat); 510} 511 512*matinv = \&PDL::matinv; 513 514sub PDL::matinv { 515 my($m) = @_; 516 my(@dims) = $m->dims; 517 518 # Keep from dumping core (FORTRAN does no error checking) 519 barf("matinv requires a 2-D square matrix") 520 unless( @dims >= 2 && $dims[0] == $dims[1] ); 521 522 $m = $m->copy(); # Make sure we don't overwrite :( 523 gefa($m,(my $ipvt=null),(my $info=null)); 524 if(sum($info) > 0) { 525 barf("Uninvertible matrix given to inv: $m\n"); 526 } 527 gedi($m,$ipvt,(pdl 0,0),(null),(long( pdl (1)))); 528 $m; 529} 530 531*detslatec = \&PDL::detslatec; 532sub PDL::detslatec { 533 my($m) = @_; 534 $m = $m->copy(); # Make sure we don't overwrite :( 535 gefa($m,(my $ipvt=null),(my $info=null)); 536 if(sum($info) > 0) { 537 barf("Uninvertible matrix given to inv: $m\n"); 538 } 539 gedi($m,$ipvt,(my $det=null),(null),(long( pdl (10)))); 540 return $det->slice('(0)')*10**$det->slice('(1)'); 541} 542 543 544sub prepfft { 545 my($n) = @_; 546 my $tmp = PDL->zeroes(float(),$n*3+15); 547 $n = pdl $n; 548 ezffti($n,$tmp); 549 return $tmp; 550} 551 552sub fft (;@) { 553 my($v) = @_; 554 my $ws = prepfft($v->getdim(0)); 555 ezfftf($v,(my $az = PDL->null), (my $a = PDL->null), 556 (my $b = PDL->null), $ws); 557 return ($az,$a,$b); 558} 559 560sub rfft { 561 my($az,$a,$b) = @_; 562 my $ws = prepfft($a->getdim(0)); 563 my $v = $a->copy(); 564 ezfftb($v,$az,$a,$b,$ws); 565 return $v; 566} 567 568# polynomial fitting routines 569# simple wrappers around the SLATEC implementations 570 571*polyfit = \&PDL::polyfit; 572sub PDL::polyfit { 573 barf 'Usage: polyfit($x, $y, $w, $maxdeg, [$eps]);' 574 unless (@_ == 5 || @_==4 ); 575 576 my ($x_in, $y_in, $w_in, $maxdeg_in, $eps_in) = @_; 577 578 # if $w_in does not match the data vectors ($x_in, $y_in), then we can resize 579 # it to match the size of $y_in : 580 $w_in = $w_in + $y_in->zeros; 581 582 # Create the output arrays 583 my $r = PDL->null; 584 585 # A array needs some work space 586 my $sz = ((3 * $x_in->getdim(0)) + (3*$maxdeg_in) + 3); # Buffer size called for by Slatec 587 my @otherdims = $_[0]->dims; shift @otherdims; # Thread dims 588 my $a = PDL::new_from_specification('PDL',$x_in->type,$sz,@otherdims); 589 my $coeffs = PDL::new_from_specification('PDL',$x_in->type, $maxdeg_in + 1, @otherdims); 590 591 my $ierr = PDL->null; 592 my $ndeg = PDL->null; 593 594 # Now call polfit 595 my $rms = pdl($eps_in); 596 polfit($x_in, $y_in, $w_in, $maxdeg_in, $ndeg, $rms, $r, $ierr, $a, $coeffs); 597 # Preserve historic compatibility by flowing rms error back into the argument 598 if( UNIVERSAL::isa($eps_in,'PDL') ){ 599 $eps_in .= $rms; 600 } 601 602 # Return the arrays 603 if(wantarray) { 604 return ($ndeg, $r, $ierr, $a, $coeffs, $rms ); 605 } else { 606 return $coeffs; 607 } 608} 609 610 611*polycoef = \&PDL::polycoef; 612sub PDL::polycoef { 613 barf 'Usage: polycoef($l, $c, $a);' 614 unless @_ == 3; 615 616 # Allocate memory for return PDL 617 # Simply l + 1 but cant see how to get PP to do this - TJ 618 # Not sure the return type since I do not know 619 # where PP will get the information from 620 my $tc = PDL->zeroes( abs($_[0]) + 1 ); 621 622 # Run the slatec routine 623 pcoef($_[0], $_[1], $tc, $_[2]); 624 625 # Return results 626 return $tc; 627 628} 629 630*polyvalue = \&PDL::polyvalue; 631sub PDL::polyvalue { 632 barf 'Usage: polyvalue($l, $nder, $x, $a);' 633 unless @_ == 4; 634 635 # Two output arrays 636 my $yfit = PDL->null; 637 638 # This one must be preallocated and must take into account 639 # the size of $x if greater than 1 640 my $yp; 641 if ($_[2]->getdim(0) == 1) { 642 $yp = $_[2]->zeroes($_[1]); 643 } else { 644 $yp = $_[2]->zeroes($_[1], $_[2]->getdim(0)); 645 } 646 647 # Run the slatec function 648 pvalue($_[0], $_[2], $yfit, $yp, $_[3]); 649 650 # Returns 651 return ($yfit, $yp); 652 653} 654 655END 656 657defslatec( 658 'svdc',{S => 'ssvdc'}, 659 'Mat x (n,p); 660 MDim ldx; 661 Dim n; 662 Dim p; 663 Mat [o] s (p); 664 Mat [o] e (p); 665 Mat [o] u (n,p); 666 MDim ldu; 667 Mat [o] v (p,p); 668 MDim ldv; 669 Mat [o] work (n); 670 IntFlag job (); 671 IntFlag [o] info (); 672 ', 673'singular value decomposition of a matrix' 674); 675 676defslatec( 677 'poco',{S => 'spoco', D => 'dpoco'}, 678 'Mat a (n,n); 679 MDim lda; 680 Dim n; 681 Mat rcond (); 682 Mat [o] z (n); 683 IntFlag [o] info (); 684 ', 685'Factor a real symmetric positive definite matrix 686and estimate the condition number of the matrix.' 687); 688 689defslatec( 690 'geco',{S => 'sgeco', D => 'dgeco'}, 691 'Mat a (n,n); 692 MDim lda; 693 Dim n; 694 IntFlag [o] ipvt (n); 695 Mat [o] rcond (); 696 Mat [o] z (n); 697 ', 698'Factor a matrix using Gaussian elimination and estimate 699the condition number of the matrix.' 700); 701 702defslatec( 703 'gefa',{S => 'sgefa', D => 'dgefa'}, 704 'Mat a (n,n); 705 MDim lda; 706 Dim n; 707 IntFlag [o] ipvt (n); 708 IntFlag [o] info (); 709 ', 710'Factor a matrix using Gaussian elimination.' 711); 712 713# XXX Ensure two == 2!! 714# 715# pofa and sqrdc aren't (yet?) implemented 716# 717defslatec( 718 'podi',{S => 'spodi', D => 'dpodi'}, 719 'Mat a (n,n); 720 MDim lda; 721 Dim n; 722 Mat [o] det (two=2); 723 IntFlag job (); 724 ', 725'Compute the determinant and inverse of a certain real 726symmetric positive definite matrix using the factors 727computed by L<poco|/poco>.' 728); 729 730defslatec( 731 'gedi',{S => 'sgedi', D => 'dgedi'}, 732 'Mat a (n,n); 733 MDim lda; 734 Dim n; 735 IntFlag [o] ipvt (n); 736 Mat [o] det (two=2); 737 Mat [o] work (n); 738 IntFlag job (); 739 ', 740'Compute the determinant and inverse of a matrix using the 741factors computed by L<geco|/geco> or L<gefa|/gefa>.' 742); 743 744 745defslatec( 746 'gesl',{S => 'sgesl', D => 'dgesl'}, 747 'Mat a (lda,n); 748 MDim lda; 749 Dim n; 750 IntFlag ipvt (n); 751 Mat b (n); 752 IntFlag job (); 753 ', 754'Solve the real system C<A*X=B> or C<TRANS(A)*X=B> using the 755factors computed by L<geco|/geco> or L<gefa|/gefa>.' 756); 757 758defslatec( 759 'rs', {S => 'rsfoo'}, 760 'MDim lda; 761 Dim n; 762 Mat a (n,n); 763 Mat [o] w (n); 764 IntFlag matz (); 765 Mat [o] z (n,n); 766 Mat [t] fvone (n); 767 Mat [t] fvtwo (n); 768 IntFlag [o] ierr (); 769 ', 770'This subroutine calls the recommended sequence of 771subroutines from the eigensystem subroutine package (EISPACK) 772to find the eigenvalues and eigenvectors (if desired) 773of a REAL SYMMETRIC matrix.' 774 775); 776 777# XXX wsave : at least 3n+15 778defslatec( 779 'ezffti', {S => 'ezffti'}, 780 'IntFlag n (); 781 Mat [o] wsave(foo); 782 ', 783'Subroutine ezffti initializes the work array C<wsave()> 784which is used in both L<ezfftf|/ezfftf> and 785L<ezfftb|/ezfftb>. 786The prime factorization 787of C<n> together with a tabulation of the trigonometric functions 788are computed and stored in C<wsave()>.' 789 790); 791 792# XXX Correct for azero, a and b 793defslatec( 794 'ezfftf', {S => 'ezfftf'}, 795 'Dim n; 796 Mat r(n); 797 Mat [o] azero(); 798 Mat [o] a(n); 799 Mat [o] b(n); 800 Mat wsave(foo); 801 ' 802); 803 804defslatec( 805 'ezfftb', {S => 'ezfftb'}, 806 'Dim n; 807 Mat [o] r(n); 808 Mat azero(); 809 Mat a(n); 810 Mat b(n); 811 Mat wsave(foo); 812 ' 813); 814 815################################################################## 816################################################################## 817 818defslatec( 819 'pcoef', {S => 'pcoef', D => 'dpcoef'}, 820 ' 821 IntFlag l (); 822 Mat c (); 823 Mat [o] tc (bar); 824 Mat a (foo); 825 ', 826'Convert the C<polfit> coefficients to Taylor series form. 827C<c> and C<a()> must be of the same type.' 828); 829 830defslatec( 831 'pvalue', {S => 'pvalue', D => 'dp1vlu'}, 832 ' 833 IntFlag l (); 834 Dim nder; 835 Mat x (); 836 Mat [o] yfit (); 837 Mat [o] yp (nder); 838 Mat a (foo); 839 ', 840'Use the coefficients generated by C<polfit> to evaluate the 841polynomial fit of degree C<l>, along with the first C<nder> of 842its derivatives, at a specified point. C<x> and C<a> must be of the 843same type.' 844); 845 846################################################################## 847################################################################## 848# 849# PCHIP library 850# 851defslatec( 852 'chim', {S => 'pchim', D => 'dpchim'}, 853 'Dim n; 854 Mat x (n); 855 Mat f (n); 856 Mat [o] d (n); 857 Incfd dummy; 858 IntFlag [o] ierr (); 859 ', 860'Calculate the derivatives at the given set of points (C<$x,$f>, 861where C<$x> is strictly increasing). 862The resulting set of points - C<$x,$f,$d>, referred to 863as the cubic Hermite representation - can then be used in 864other functions, such as L<chfe|/chfe>, L<chfd|/chfd>, 865and L<chia|/chia>. 866 867The boundary conditions are compatible with monotonicity, 868and if the data are only piecewise monotonic, the interpolant 869will have an extremum at the switch points; for more control 870over these issues use L<chic|/chic>. 871 872Error status returned by C<$ierr>: 873 874=over 4 875 876=item * 877 8780 if successful. 879 880=item * 881 882E<gt> 0 if there were C<ierr> switches in the direction of 883monotonicity (data still valid). 884 885=item * 886 887-1 if C<nelem($x) E<lt> 2>. 888 889=item * 890 891-3 if C<$x> is not strictly increasing. 892 893=back 894 895=cut 896 897', 898'Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.' 899); 900 901# switch has become mflag, since `switch' is a reserved word in 902# C. 903# 904# can not say (nwk=2*n) --- the rhs has to equal a number 905# -> could Basic/Gen/PP/Dims.pm be hacked to allow this? 906# 907# I didn't have much success with preceding wk by [t] 908# 909defslatec( 910 'chic', {S => 'pchic', D => 'dpchic'}, 911 'IntFlag ic (two=2); 912 Mat vc (two=2); 913 Mat mflag (); 914 Dim n; 915 Mat x (n); 916 Mat f (n); 917 Mat [o] d (n); 918 Incfd dummy; 919 Mat wk (nwk); 920 Dim nwk; 921 IntFlag [o] ierr (); 922 ', 923'Calculate the derivatives at the given points (C<$x,$f>, 924where C<$x> is strictly increasing). 925Control over the boundary conditions is given by the 926C<$ic> and C<$vc> piddles, and the value of C<$mflag> determines 927the treatment of points where monotoncity switches 928direction. A simpler, more restricted, interface is available 929using L<chim|/chim>. 930 931The first and second elements of C<$ic> determine the boundary 932conditions at the start and end of the data respectively. 933If the value is 0, then the default condition, as used by 934L<chim|/chim>, is adopted. 935If greater than zero, no adjustment for monotonicity is made, 936otherwise if less than zero the derivative will be adjusted. 937The allowed magnitudes for C<ic(0)> are: 938 939=over 4 940 941=item * 942 9431 if first derivative at C<x(0)> is given in C<vc(0)>. 944 945=item * 946 9472 if second derivative at C<x(0)> is given in C<vc(0)>. 948 949=item * 950 9513 to use the 3-point difference formula for C<d(0)>. 952(Reverts to the default b.c. if C<n E<lt> 3>) 953 954=item * 955 9564 to use the 4-point difference formula for C<d(0)>. 957(Reverts to the default b.c. if C<n E<lt> 4>) 958 959=item * 960 9615 to set C<d(0)> so that the second derivative is 962continuous at C<x(1)>. 963(Reverts to the default b.c. if C<n E<lt> 4>) 964 965=back 966 967The values for C<ic(1)> are the same as above, except that 968the first-derivative value is stored in C<vc(1)> for cases 1 and 2. 969The values of C<$vc> need only be set if options 1 or 2 are chosen 970for C<$ic>. 971 972Set C<$mflag = 0> if interpolant is required to be monotonic in 973each interval, regardless of the data. This causes C<$d> to be 974set to 0 at all switch points. Set C<$mflag> to be non-zero to 975use a formula based on the 3-point difference formula at switch 976points. If C<$mflag E<gt> 0>, then the interpolant at swich points 977is forced to not deviate from the data by more than C<$mflag*dfloc>, 978where C<dfloc> is the maximum of the change of C<$f> on this interval 979and its two immediate neighbours. 980If C<$mflag E<lt> 0>, no such control is to be imposed. 981 982The piddle C<$wk> is only needed for work space. However, I could 983not get it to work as a temporary variable, so you must supply 984it; it is a 1D piddle with C<2*n> elements. 985 986Error status returned by C<$ierr>: 987 988=over 4 989 990=item * 991 9920 if successful. 993 994=item * 995 9961 if C<ic(0) E<lt> 0> and C<d(0)> had to be adjusted for 997monotonicity. 998 999=item * 1000 10012 if C<ic(1) E<lt> 0> and C<d(n-1)> had to be adjusted 1002for monotonicity. 1003 1004=item * 1005 10063 if both 1 and 2 are true. 1007 1008=item * 1009 1010-1 if C<n E<lt> 2>. 1011 1012=item * 1013 1014-3 if C<$x> is not strictly increasing. 1015 1016=item * 1017 1018-4 if C<abs(ic(0)) E<gt> 5>. 1019 1020=item * 1021 1022-5 if C<abs(ic(1)) E<gt> 5>. 1023 1024=item * 1025 1026-6 if both -4 and -5 are true. 1027 1028=item * 1029 1030-7 if C<nwk E<lt> 2*(n-1)>. 1031 1032=back 1033 1034=cut 1035 1036', 1037'Calculate the derivatives of (x,f(x)) using cubic Hermite interpolation.' 1038); 1039 1040# as above, have made wk an actual piddle, rather than a [t] 1041defslatec( 1042 'chsp', {S => 'pchsp', D => 'dpchsp'}, 1043 'IntFlag ic (two=2); 1044 Mat vc (two=2); 1045 Dim n; 1046 Mat x (n); 1047 Mat f (n); 1048 Mat [o] d (n); 1049 Incfd dummy; 1050 Mat wk (nwk); 1051 Dim nwk; 1052 IntFlag [o] ierr (); 1053 ', 1054'Calculate the derivatives, using cubic spline interpolation, 1055at the given points (C<$x,$f>), with the specified 1056boundary conditions. 1057Control over the boundary conditions is given by the 1058C<$ic> and C<$vc> piddles. 1059The resulting values - C<$x,$f,$d> - can 1060be used in all the functions which expect a cubic 1061Hermite function. 1062 1063The first and second elements of C<$ic> determine the boundary 1064conditions at the start and end of the data respectively. 1065The allowed values for C<ic(0)> are: 1066 1067=over 4 1068 1069=item * 1070 10710 to set C<d(0)> so that the third derivative is 1072continuous at C<x(1)>. 1073 1074=item * 1075 10761 if first derivative at C<x(0)> is given in C<vc(0>). 1077 1078=item * 1079 10802 if second derivative at C<x(0>) is given in C<vc(0)>. 1081 1082=item * 1083 10843 to use the 3-point difference formula for C<d(0)>. 1085(Reverts to the default b.c. if C<n E<lt> 3>.) 1086 1087=item * 1088 10894 to use the 4-point difference formula for C<d(0)>. 1090(Reverts to the default b.c. if C<n E<lt> 4>.) 1091 1092=back 1093 1094The values for C<ic(1)> are the same as above, except that 1095the first-derivative value is stored in C<vc(1)> for cases 1 and 2. 1096The values of C<$vc> need only be set if options 1 or 2 are chosen 1097for C<$ic>. 1098 1099The piddle C<$wk> is only needed for work space. However, I could 1100not get it to work as a temporary variable, so you must supply 1101it; it is a 1D piddle with C<2*n> elements. 1102 1103Error status returned by C<$ierr>: 1104 1105=over 4 1106 1107=item * 1108 11090 if successful. 1110 1111=item * 1112 1113-1 if C<nelem($x) E<lt> 2>. 1114 1115=item * 1116 1117-3 if C<$x> is not strictly increasing. 1118 1119=item * 1120 1121-4 if C<ic(0) E<lt> 0> or C<ic(0) E<gt> 4>. 1122 1123=item * 1124 1125-5 if C<ic(1) E<lt> 0> or C<ic(1) E<gt> 4>. 1126 1127=item * 1128 1129-6 if both of the above are true. 1130 1131=item * 1132 1133-7 if C<nwk E<lt> 2*n>. 1134 1135=item * 1136 1137-8 in case of trouble solving the linear system 1138for the interior derivative values. 1139 1140=back 1141 1142=cut 1143 1144', 1145'Calculate the derivatives of (x,f(x)) using cubic spline interpolation.' 1146); 1147 1148defslatec( 1149 'chfd', {S => 'pchfd', D => 'dpchfd'}, 1150 'Dim n; 1151 Mat x (n); 1152 Mat f (n); 1153 Mat d (n); 1154 Incfd dummy; 1155 CheckFlag check (); 1156 Dim ne; 1157 Mat xe (ne); 1158 Mat [o] fe (ne); 1159 Mat [o] de (ne); 1160 IntFlag [o] ierr (); 1161 ', 1162'Given a piecewise cubic Hermite function - such as from 1163L<chim|/chim> - evaluate the function (C<$fe>) and 1164derivative (C<$de>) at a set of points (C<$xe>). 1165If function values alone are required, use L<chfe|/chfe>. 1166Set C<check> to 0 to skip checks on the input data. 1167 1168Error status returned by C<$ierr>: 1169 1170=over 4 1171 1172=item * 1173 11740 if successful. 1175 1176=item * 1177 1178E<gt>0 if extrapolation was performed at C<ierr> points 1179(data still valid). 1180 1181=item * 1182 1183-1 if C<nelem($x) E<lt> 2> 1184 1185=item * 1186 1187-3 if C<$x> is not strictly increasing. 1188 1189=item * 1190 1191-4 if C<nelem($xe) E<lt> 1>. 1192 1193=item * 1194 1195-5 if an error has occurred in a lower-level routine, 1196which should never happen. 1197 1198=back 1199 1200=cut 1201 1202', 1203'Interpolate function and derivative values.' 1204); 1205 1206defslatec( 1207 'chfe', {S => 'pchfe', D => 'dpchfe'}, 1208 'Dim n; 1209 Mat x (n); 1210 Mat f (n); 1211 Mat d (n); 1212 Incfd dummy; 1213 CheckFlag check (); 1214 Dim ne; 1215 Mat xe (ne); 1216 Mat [o] fe (ne); 1217 IntFlag [o] ierr (); 1218 ', 1219'Given a piecewise cubic Hermite function - such as from 1220L<chim|/chim> - evaluate the function (C<$fe>) at 1221a set of points (C<$xe>). 1222If derivative values are also required, use L<chfd|/chfd>. 1223Set C<check> to 0 to skip checks on the input data. 1224 1225Error status returned by C<$ierr>: 1226 1227=over 4 1228 1229=item * 1230 12310 if successful. 1232 1233=item * 1234 1235E<gt>0 if extrapolation was performed at C<ierr> points 1236(data still valid). 1237 1238=item * 1239 1240-1 if C<nelem($x) E<lt> 2> 1241 1242=item * 1243 1244-3 if C<$x> is not strictly increasing. 1245 1246=item * 1247 1248-4 if C<nelem($xe) E<lt> 1>. 1249 1250=back 1251 1252=cut 1253 1254', 1255'Interpolate function values.' 1256); 1257 1258defslatec( 1259 'chia', {S => 'pchia', D => 'dpchia'}, 1260 'Dim n; 1261 Mat x (n); 1262 Mat f (n); 1263 Mat d (n); 1264 Incfd dummy; 1265 CheckFlag check (); 1266 Mat a (); 1267 Mat b (); 1268 FuncRet [o] ans (); 1269 IntFlag [o] ierr (); 1270 ', 1271'Evaluate the definite integral of a a piecewise 1272cubic Hermite function over an arbitrary interval, 1273given by C<[$a,$b]>. 1274See L<chid|/chid> if the integration limits are 1275data points. 1276Set C<check> to 0 to skip checks on the input data. 1277 1278The values of C<$a> and C<$b> do not have 1279to lie within C<$x>, although the resulting integral 1280value will be highly suspect if they are not. 1281 1282Error status returned by C<$ierr>: 1283 1284=over 4 1285 1286=item * 1287 12880 if successful. 1289 1290=item * 1291 12921 if C<$a> lies outside C<$x>. 1293 1294=item * 1295 12962 if C<$b> lies outside C<$x>. 1297 1298=item * 1299 13003 if both 1 and 2 are true. 1301 1302=item * 1303 1304-1 if C<nelem($x) E<lt> 2> 1305 1306=item * 1307 1308-3 if C<$x> is not strictly increasing. 1309 1310=item * 1311 1312-4 if an error has occurred in a lower-level routine, 1313which should never happen. 1314 1315=back 1316 1317=cut 1318 1319', 1320'Integrate (x,f(x)) over arbitrary limits.' 1321); 1322 1323defslatec( 1324 'chid', {S => 'pchid', D => 'dpchid'}, 1325 'Dim n; 1326 Mat x (n); 1327 Mat f (n); 1328 Mat d (n); 1329 Incfd dummy; 1330 CheckFlag check (); 1331 FortranIndex ia (); 1332 FortranIndex ib (); 1333 FuncRet [o] ans (); 1334 IntFlag [o] ierr (); 1335 ', 1336'Evaluate the definite integral of a a piecewise 1337cubic Hermite function between C<x($ia)> and 1338C<x($ib)>. 1339 1340See L<chia|/chia> for integration between arbitrary 1341limits. 1342 1343Although using a fortran routine, the values of 1344C<$ia> and C<$ib> are zero offset. 1345Set C<check> to 0 to skip checks on the input data. 1346 1347Error status returned by C<$ierr>: 1348 1349=over 4 1350 1351=item * 1352 13530 if successful. 1354 1355=item * 1356 1357-1 if C<nelem($x) E<lt> 2>. 1358 1359=item * 1360 1361-3 if C<$x> is not strictly increasing. 1362 1363=item * 1364 1365-4 if C<$ia> or C<$ib> is out of range. 1366 1367=back 1368 1369=cut 1370 1371', 1372'Integrate (x,f(x)) between data points.' 1373); 1374 1375defslatec( 1376 'chcm', {S => 'pchcm', D => 'dpchcm'}, 1377 'Dim n; 1378 Mat x (n); 1379 Mat f (n); 1380 Mat d (n); 1381 Incfd dummy; 1382 CheckFlag check (); 1383 IntFlag [o] ismon (n); 1384 IntFlag [o] ierr (); 1385 ', 1386'The outout piddle C<$ismon> indicates over 1387which intervals the function is monotonic. 1388Set C<check> to 0 to skip checks on the input data. 1389 1390For the data interval C<[x(i),x(i+1)]>, the 1391values of C<ismon(i)> can be: 1392 1393=over 4 1394 1395=item * 1396 1397-3 if function is probably decreasing 1398 1399=item * 1400 1401-1 if function is strictly decreasing 1402 1403=item * 1404 14050 if function is constant 1406 1407=item * 1408 14091 if function is strictly increasing 1410 1411=item * 1412 14132 if function is non-monotonic 1414 1415=item * 1416 14173 if function is probably increasing 1418 1419=back 1420 1421If C<abs(ismon(i)) == 3>, the derivative values are 1422near the boundary of the monotonicity region. A small 1423increase produces non-monotonicity, whereas a decrease 1424produces strict monotonicity. 1425 1426The above applies to C<i = 0 .. nelem($x)-1>. The last element of 1427C<$ismon> indicates whether 1428the entire function is monotonic over $x. 1429 1430Error status returned by C<$ierr>: 1431 1432=over 4 1433 1434=item * 1435 14360 if successful. 1437 1438=item * 1439 1440-1 if C<n E<lt> 2>. 1441 1442=item * 1443 1444-3 if C<$x> is not strictly increasing. 1445 1446=back 1447 1448=cut 1449 1450', 1451'Check the given piecewise cubic Hermite function for monotonicity.' 1452); 1453 1454=pod 1455 1456ignore function for now although code is in slatec/ directory 1457 1458=cut 1459 1460# XXX tsize = 2*n+4 1461# bsize = 2*n 1462# 1463# ndim gets set to 2*n 1464# 1465# Changed by routine: 1466# nknots 1467# t 1468defslatec( 1469 'chbs', {S => 'pchbs', D => 'dpchbs'}, 1470 'Dim n; 1471 Mat x (n); 1472 Mat f (n); 1473 Mat d (n); 1474 Incfd dummy; 1475 IntFlag knotyp (); 1476 IntFlag nknots (); 1477 Mat t (tsize); 1478 Mat [o] bcoef (bsize); 1479 IntFlag [o] ndim (); 1480 IntFlag [o] kord (); 1481 IntFlag [o] ierr (); 1482 ', 1483'The resulting B-spline representation of the data 1484(i.e. C<nknots>, C<t>, C<bcoeff>, C<ndim>, and 1485C<kord>) can be evaluated by C<bvalu> (which is 1486currently not available). 1487 1488Array sizes: C<tsize = 2*n + 4>, C<bsize = 2*n>, 1489and C<ndim = 2*n>. 1490 1491C<knotyp> is a flag which controls the knot sequence. 1492The knot sequence C<t> is normally computed from C<$x> 1493by putting a double knot at each C<x> and setting the end knot pairs 1494according to the value of C<knotyp> (where C<m = ndim = 2*n>): 1495 1496=over 1497 1498=item * 1499 15000 - Quadruple knots at the first and last points. 1501 1502=item * 1503 15041 - Replicate lengths of extreme subintervals: 1505C<t( 0 ) = t( 1 ) = x(0) - (x(1)-x(0))> and 1506C<t(m+3) = t(m+2) = x(n-1) + (x(n-1)-x(n-2))> 1507 1508=item * 1509 15102 - Periodic placement of boundary knots: 1511C<t( 0 ) = t( 1 ) = x(0) - (x(n-1)-x(n-2))> and 1512C<t(m+3) = t(m+2) = x(n) + (x(1)-x(0))> 1513 1514=item * 1515 1516E<lt>0 - Assume the C<nknots> and C<t> were set in a previous call. 1517 1518=back 1519 1520C<nknots> is the number of knots and may be changed by the routine. 1521If C<knotyp E<gt>= 0>, C<nknots> will be set to C<ndim+4>, 1522otherwise it is an input variable, and an error will occur if its 1523value is not equal to C<ndim+4>. 1524 1525C<t> is the array of C<2*n+4> knots for the B-representation 1526and may be changed by the routine. 1527If C<knotyp E<gt>= 0>, C<t> will be changed so that the 1528interior double knots are equal to the x-values and the 1529boundary knots set as indicated above, 1530otherwise it is assumed that C<t> was set by a 1531previous call (no check is made to verify that the data 1532forms a legitimate knot sequence). 1533 1534Error status returned by C<$ierr>: 1535 1536=over 4 1537 1538=item * 1539 15400 if successful. 1541 1542=item * 1543 1544-4 if C<knotyp E<gt> 2>. 1545 1546=item * 1547 1548-5 if C<knotyp E<lt> 0> and C<nknots != 2*n + 4>. 1549 1550=back 1551 1552=cut 1553 1554', 1555'Piecewise Cubic Hermite function to B-Spline converter.' 1556); 1557 1558################################################################## 1559################################################################## 1560 1561# 1562# This version of polfit accepts bad values and also allows threading 1563# 1564 1565# 1566# indices: 1567# n runs across input points; 1568# foo runs across wacky Slatec buffer size; 1569# bar runs across polynomial coefficients. 1570# 1571pp_def('polfit', 1572 Pars => '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)', 1573 OtherPars => '', 1574 Code => ' 1575 int maxord; 1576 int ord; 1577 int k; 1578 $GENERIC() zero = 0; 1579 1580 $TFD(polfit'.$uscore.',dpolft'.$uscore.') (&$PRIV(__n_size),$P(x),$P(y),$P(w),$P(maxdeg),$P(ndeg),$P(eps),$P(r),$P(ierr),$P(a)); 1581 maxord = ($P(a))[0]+0.5; 1582 ord = ($P(a))[maxord * 3 + 2]; 1583 if(ord >= $maxdeg()) { 1584 ord = $maxdeg(); 1585 } 1586 $TFD(pcoef'.$uscore.',dpcoef'.$uscore.') ( &(ord), &(zero), $P(coeffs), $P(a)); 1587 for(k=ord+1; k<=$maxdeg(); k++) 1588 ($P(coeffs))[k] = 0; 1589', 1590 GenericTypes => ['F','D'], 1591 HandleBad => 1, 1592 NoBadifNaN => 1, 1593 BadCode => 'int ns = $SIZE(n); 1594 int i; 1595 int j = 0; 1596 if($SIZE(n)<$maxdeg()) { 1597 barf("polfit: Must have at least <n> points to fit <n> coefficients"); 1598 } 1599 1600 for (i=0;i<ns;i++) { /* get rid of bad values. Call polfit with [xyw]tmp instead of [xyz]. */ 1601 if ($ISGOOD(y(n=>i)) && $ISGOOD(x(n=>i)) && $ISGOOD(w(n=>i))) { 1602 $xtmp(n=>j) = $x(n=>i); 1603 $ytmp(n=>j) = $y(n=>i); 1604 $wtmp(n=>j) = $w(n=>i); 1605 j++; 1606 } 1607 } 1608 if (j <= $maxdeg()) { 1609 /* Not enough good datapoints -- set this whole row to BAD. */ 1610 for (i=0;i<ns;i++) { 1611 $SETBAD(r(n=>i)); 1612 } 1613 $ierr() = 2; 1614 } else { 1615 /* Found enough good datapoints for a fit -- do the fit */ 1616 int k; 1617 int ord; 1618 int maxord; 1619 $GENERIC() zero = 0; 1620 1621 /* Do the fit */ 1622 $TFD(polfit'.$uscore.',dpolft'.$uscore.') 1623 (&j,$P(xtmp),$P(ytmp),$P(wtmp),$P(maxdeg),$P(ndeg),$P(eps),$P(rtmp),$P(ierr),$P(a)); 1624 1625 maxord = ($P(a))[0]+0.5; 1626 ord = ($P(a))[maxord * 3 + 2]; 1627 if(ord >= $maxdeg()) { 1628 ord = $maxdeg(); 1629 } 1630 /* Extract the polynomial coefficients into coeffs -- used for bad values */ 1631 $TFD(pcoef'.$uscore.',dpcoef'.$uscore.') ( &(ord), &(zero), $P(coeffs), $P(a)); 1632 for(k=ord+1; k<=$maxdeg(); k++) 1633 ($P(coeffs))[k] = 0; 1634 j=0; 1635 for (i=0;i<ns;i++) { /* replace bad values */ 1636 if ($ISGOOD(y(n=>i))) { 1637 $r(n=>i) = $rtmp(n=>j); 1638 j++; 1639 } else if($ISGOOD(x(n=>i))) { 1640 /* Bad values are omitted from the call to polfit, so we must reconstitute them on return */ 1641 /* (just because a value is bad in y, does not mean the fit itself is bad there) */ 1642 /* */ 1643 /* The number in ord is not the number of coefficients in the polynomial, it is the highest */ 1644 /* order coefficient -- so 3 for a cubic, which has 4 coefficients. */ 1645 /* --CED */ 1646 int ii; 1647 $GENERIC() acc = 0; 1648 for( ii=ord; ii>0; ii-- ) { 1649 acc += $coeffs(bar=>ii); 1650 acc *= $x(n=>i); 1651 } 1652 1653 acc += $coeffs(bar=>0); 1654 $r(n=>i) = acc; 1655 } else { 1656 /* x and y are bad here... */ 1657 $SETBAD(r(n=>i)); 1658 } 1659 } 1660 }', 1661 1662 Doc => 'Fit discrete data in a least squares sense by polynomials 1663 in one variable. C<x()>, C<y()> and C<w()> must be of the same type. 1664 This version handles bad values appropriately', 1665); 1666 1667#these two need to be done manually because we don't use defslatec for them 1668print PROTOS "extern int polfit_ ();\nextern int dpolft_ ();\n"; 1669 1670close( PROTOS ); 1671 1672################################################################## 1673################################################################## 1674 1675pp_addpm(<<'EOD'); 1676 1677=head1 AUTHOR 1678 1679Copyright (C) 1997 Tuomas J. Lukka. 1680Copyright (C) 2000 Tim Jenness, Doug Burke. 1681All rights reserved. There is no warranty. You are allowed 1682to redistribute this software / documentation under certain 1683conditions. For details, see the file COPYING in the PDL 1684distribution. If this file is separated from the PDL distribution, 1685the copyright notice should be included in the file. 1686 1687=cut 1688 1689 1690EOD 1691 1692pp_done(); 1693 1694 1695