From 1b2da5624be57caaadcc6b63bf20f07fdce93cdc Mon Sep 17 00:00:00 2001 From: Ed J Date: Tue, 17 Dec 2024 13:00:42 +0000 Subject: [PATCH] make chfe into PP kernel, remove pchip.c and defslatec --- MANIFEST | 1 - lib/PDL/Primitive-pchip.c | 218 -------------------------------------- lib/PDL/Primitive.pd | 176 ++++-------------------------- 3 files changed, 21 insertions(+), 374 deletions(-) delete mode 100644 lib/PDL/Primitive-pchip.c diff --git a/MANIFEST b/MANIFEST index 00d34b78e..e74ab6342 100644 --- a/MANIFEST +++ b/MANIFEST @@ -146,7 +146,6 @@ lib/PDL/PP/Dims.pm lib/PDL/PP/PDLCode.pm lib/PDL/PP/PdlParObj.pm lib/PDL/PP/Signature.pm -lib/PDL/Primitive-pchip.c lib/PDL/Primitive-xoshiro256plus.c lib/PDL/Primitive.pd lib/PDL/QuickStart.pod diff --git a/lib/PDL/Primitive-pchip.c b/lib/PDL/Primitive-pchip.c deleted file mode 100644 index 12a2d0481..000000000 --- a/lib/PDL/Primitive-pchip.c +++ /dev/null @@ -1,218 +0,0 @@ -#include -#include -#include - -/* f2c.h -- Standard Fortran to C header file */ -typedef long int integer; -typedef double doublereal; -typedef long int logical; - -#define TRUE_ (1) - -#define min(a,b) ((a) <= (b) ? (a) : (b)) -#define max(a,b) ((a) >= (b) ? (a) : (b)) -/* end f2c.h */ - -#define xermsg_(lib, func, errstr, nerr, ...) \ - fprintf(stderr, "%s::%s: %s (err=%ld)\n", lib, func, errstr, nerr) - -/* Programming notes: */ -/* 2. Most of the coding between the call to DCHFDV and the end of */ -/* the IR-loop could be eliminated if it were permissible to */ -/* assume that XE is ordered relative to X. */ -/* 3. DCHFDV does not assume that X1 is less than X2. thus, it would */ -/* be possible to write a version of DPCHFD that assumes a strict- */ -/* ly decreasing X-array by simply running the IR-loop backwards */ -/* (and reversing the order of appropriate tests). */ -/* 4. The present code has a minor bug, which I have decided is not */ -/* worth the effort that would be required to fix it. */ -/* If XE contains points in [X(N-1),X(N)], followed by points .LT. */ -/* X(N-1), followed by points .GT.X(N), the extrapolation points */ -/* will be counted (at least) twice in the total returned in IERR. */ -void dpchfd(integer n, doublereal *x, doublereal *f, - doublereal *d, logical *skip, integer ne, - doublereal *xe, doublereal *fe, doublereal *de, integer *ierr) -{ -/* Local variables */ - integer i, j, nj, ir, ierc, next[2], jfirst, located; -/* VALIDITY-CHECK ARGUMENTS. */ - if (ne < 1) { - *ierr = -4; - xermsg_("SLATEC", "DPCHFD", "NUMBER OF EVALUATION POINTS LESS THAN ONE", *ierr); - return; - } - if (!*skip) { - if (n < 2) { - *ierr = -1; - xermsg_("SLATEC", "DPCHFD", "NUMBER OF DATA POINTS LESS THAN TWO", *ierr); - return; - } - for (i = 1; i < n; ++i) { - if (x[i] <= x[i - 1]) { - *ierr = -3; - xermsg_("SLATEC", "DPCHFD", "X-ARRAY NOT STRICTLY INCREASING", *ierr); - return; - } - } - } -/* FUNCTION DEFINITION IS OK, GO ON. */ - *ierr = 0; - *skip = TRUE_; -/* LOOP OVER INTERVALS. ( INTERVAL INDEX IS IL = IR-1 . ) */ -/* ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) */ - jfirst = 0; - ir = 1; - while (1) { -/* SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. */ - if (jfirst >= ne) { - return; - } -/* LOCATE ALL POINTS IN INTERVAL. */ - located = 0; - for (j = jfirst; j < ne; ++j) { - if (xe[j] >= x[ir]) { - located = 1; - break; - } - } - if (!located || ir == n-1) { - j = ne; - } -/* HAVE LOCATED FIRST POINT BEYOND INTERVAL. */ - nj = j - jfirst; -/* SKIP EVALUATION IF NO POINTS IN INTERVAL. */ - if (nj == 0) { - ++ir; - if (ir < n) { - continue; - } - return; - } -/* EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 . */ -/* ---------------------------------------------------------------- */ - do { /* inline dchfdv */ -/* Local variables */ - doublereal x1 = x[ir - 1], x2 = x[ir]; - doublereal f1 = f[ir - 1], f2 = f[ir]; - doublereal d1 = d[ir - 1], d2 = d[ir]; - integer ne = nj; - doublereal *xe2 = &xe[jfirst], *fe2 = &fe[jfirst]; - doublereal *de2 = de ? &de[jfirst] : NULL; - doublereal h; - integer i; - doublereal x, c2, c3, c2t2, c3t3, xma, xmi, del1, del2, delta; - if (ne < 1) { - ierc = -1; - xermsg_("SLATEC", "DCHFDV", "NUMBER OF EVALUATION POINTS LESS THAN ONE", ierc); - break; - } - h = x2 - x1; - if (h == 0.) { - ierc = -2; - xermsg_("SLATEC", "DCHFDV", "INTERVAL ENDPOINTS EQUAL", ierc); - break; - } -/* INITIALIZE. */ - ierc = 0; - next[0] = 0; - next[1] = 0; - xmi = min(0.,h); - xma = max(0.,h); -/* COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). */ - delta = (f2 - f1) / h; - del1 = (d1 - delta) / h; - del2 = (d2 - delta) / h; -/* (DELTA IS NO LONGER NEEDED.) */ - c2 = -(del1 + del1 + del2); - c2t2 = c2 + c2; - c3 = (del1 + del2) / h; -/* (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) */ - c3t3 = c3 + c3 + c3; -/* EVALUATION LOOP. */ - for (i = 0; i < ne; ++i) { - x = xe2[i] - x1; - fe2[i] = f1 + x * (d1 + x * (c2 + x * c3)); - if (de2) de2[i] = d1 + x * (c2t2 + x * c3t3); -/* COUNT EXTRAPOLATION POINTS. */ - if (x < xmi) { - ++next[0]; - } - if (x > xma) { - ++next[1]; - } -/* (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) */ - } - } while (0); /* end inline dchfdv */ -/* ---------------------------------------------------------------- */ - if (ierc < 0) { - *ierr = -5; - xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr); - return; - } - if (next[1] != 0) { -/* IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE */ -/* RIGHT OF X(IR). */ - if (ir < n-1) { -/* WE SHOULD NEVER HAVE GOTTEN HERE. */ - *ierr = -5; - xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr); - return; - } -/* THESE ARE ACTUALLY EXTRAPOLATION POINTS. */ - *ierr += next[1]; - } - if (next[0] != 0) { -/* IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE */ -/* LEFT OF X(IR-1). */ - if (ir < 2) { -/* THESE ARE ACTUALLY EXTRAPOLATION POINTS. */ - *ierr += next[0]; - goto L49; - } -/* XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST */ -/* EVALUATION INTERVAL. */ -/* FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). */ - located = 0; - for (i = jfirst; i < j; ++i) { - if (xe[i] < x[ir - 1]) { - located = 1; - break; - } - } - if (!located) { -/* NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR */ -/* IN DCHFDV. */ - *ierr = -5; - xermsg_("SLATEC", "DPCHFD", "ERROR RETURN FROM DCHFDV -- FATAL", *ierr); - return; - } -/* RESET J. (THIS WILL BE THE NEW JFIRST.) */ - j = i; -/* NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. */ - for (i = 0; i < ir; ++i) { - if (xe[j] < x[i]) { - break; - } - } -/* NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1). */ -/* AT THIS POINT, EITHER XE(J) .LT. X(1) */ -/* OR X(I-1) .LE. XE(J) .LT. X(I) . */ -/* RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE */ -/* CYCLING. */ -/* Computing MAX */ - ir = max(0,i-1); - } - L49: - jfirst = j; -/* END OF IR-LOOP. */ - ++ir; - if (ir >= n) break; - } -} - -void dpchfe(integer n, doublereal *x, doublereal *f, - doublereal *d, logical *skip, integer ne, - doublereal *xe, doublereal *fe, integer *ierr) -{ - dpchfd(n, x, f, d, skip, ne, xe, fe, NULL, ierr); -} diff --git a/lib/PDL/Primitive.pd b/lib/PDL/Primitive.pd index a3417411d..fdc5e1fff 100644 --- a/lib/PDL/Primitive.pd +++ b/lib/PDL/Primitive.pd @@ -6,7 +6,7 @@ my $AF = [map $_->ppsym, grep !$_->integer, types]; $AF = [(grep $_ ne 'D', @$AF), 'D']; # so defaults to D if non-float given { no warnings 'once'; # pass info back to Makefile.PL -$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= " $::PDLBASE-xoshiro256plus\$(OBJ_EXT) $::PDLBASE-pchip\$(OBJ_EXT)"; +$PDL::Core::Dev::EXTRAS{$::PDLMOD}{OBJECT} .= " $::PDLBASE-xoshiro256plus\$(OBJ_EXT)"; } pp_addpm({At=>'Top'},<<'EOD'); @@ -4133,140 +4133,6 @@ Find all numbers less that 100 that are of the form 2*y and 3*x sub PDL::intersect { setops($_[0], 'AND', $_[1]) } EOD -sub intcopy { - my ($cvarname, $rvalue) = @_; # closures - (sub {"&".$cvarname->(@_)}, sub {"PDL_LongLong *"}, - sub {"PDL_LongLong ".$cvarname->(@_)." = ".$rvalue->(@_).";"}); -} -# defslatec( $pdlname, $funcnames, $argstr, $docstring, $funcref ) -# -# $pdlname is the name of the PDL function to be created -# $funcnames is a reference to a hash array, whose keys define -# the single (S), and double precision (D) names of the -# SLATEC routines to be linked to. -# -# $argstr is a list of arguments expected by the SLATEC routine -# - some of the allowed type names are: -# FuncRet -# - specifies that this is a function, not a subroutine, and -# that the output of the function should be stored in this -# variable -# IntVal -# - used for passing in a value to func, either calculated or a constant, -# but not a caller-supplied arg -# Logical -# - A Fortran Boolean, defined as same size as REAL/INTEGER so 'long' -# -# $docstring gives the text to be used as the function documentation -# -# $funcref gets placed within a '=for ref' pod statement at the -# start of the documentation - ie it is placed before the -# text within $docstring. This string gets printed out -# in the perldl shell after a '?? string' command -# -# [Par_forcetype(undef=none,ref=OtherPar), C_argin, C_paramtype, C_decl_value] -my %KIND2I = ( - Mat => ['', sub{"\$P($_[0][2])"}, sub {PDL::Type->new($_[0])->ctype." *"}], - FuncRet => ['', sub {()}, sub {()}], - Integer => ['longlong ', sub{"\$P($_[0][2])"}, sub {"PDL_LongLong *"}], - IntScalar => ['longlong ', sub{"\$$_[0][2]()"}, sub {"PDL_LongLong "}], - MatScalar => ['', sub{"\$$_[0][2]()"}, sub {PDL::Type->new($_[0])->ctype." "}], - Logical => ['int ', sub{"\$P($_[0][2])"}, sub {"PDL_Long *"}], - OtherParDim => [\'PDL_LongLong ', intcopy(sub {"dim_$_[0][2]"}, sub {"\$SIZE($_[0][2])"})], - IntVal => [undef, sub {$_[0][4]}, sub {"PDL_LongLong "}], -); -my %ignore_ppar = map +($_=>1), qw(IntVal OtherParDim); -my %prototype = ( F => "float", D => "double" ); -my %ftypes = (S => 'F', D => 'D'); -our $debug = 0; # print out calls to pp_def; use "local $debug=1" before func -sub defslatec { - my ($pname,$fnames,$argstr,$docstring,$funcref) = @_; - my @args = grep /\S/, split ';', $argstr; - my @args2 = map { - /^\s*([a-zA-Z]+)\s+ # "Type name" - ((?:\[[^]]*\])?)\s* # Options - ([a-zA-Z]+)\s* # Par name - (?: - ((?:\(.*\))?)\s*$ # Dims - | - =\s*((?:.*?)?)\s*$ # value for IntVal - ) - /x or die("Invalid slatec par $_"); - [$1,$2,$3,$4,$5]} @args; - # is this for a function (Type name eq "FuncRet") - # or a subroutine? - die "Only one FuncRet allowed in pars list.\n" - if (my @funcret = grep $_->[0] eq "FuncRet", @args2) > 1; - my $fpar = @funcret ? "\$$funcret[0][2]()" : undef; - my $pars = join '; ', map - +($KIND2I{$_->[0]}[0] // die "Invalid ppars ",join(',',@$_),"\n") . - join(' ', grep length, $_->[1], "$_->[2]$_->[3]"), - grep !$ignore_ppar{$_->[0]}, @args2; - my $otherpars = join ';', map - ${$KIND2I{$_->[0]}[0]} . join(' => ', @$_[2,2]), - grep ref $KIND2I{$_->[0]}[0], @args2; - my $argorder; - if ($otherpars) { - my @mand = map $_->[2], grep $_->[1] !~ /\bo\b/, @args2; - my @opt = map $_->[2], grep $_->[1] =~ /\bo\b/, @args2; - $argorder = [@mand, @opt]; - } - my @talts = map - [$ftypes{$_} // die("FTYPE $_ NOT THERE\n"),$fnames->{$_}], - sort keys %$fnames; - my $func = "\$T".join('',map {$_->[0]} @talts) . "(" . - join(',',map qq{$_->[1]}, @talts).")"; - if ( defined $fpar ) { $func = "$fpar = $func"; } - my $funcargs = join ',', map $KIND2I{$_->[0]}[1]->($_), @args2; - my @protos = map $KIND2I{$_->[0]}[2], @args2; - my @cheaders; - foreach my $t ( @talts ) { - my @decl_args; - for my $ind (0..$#protos) { - my $arg = $protos[$ind]; - next unless my @ret = $arg->($t->[0]); - push @decl_args, "$ret[0]$args2[$ind][2]"; - } - push @cheaders, (defined $fpar ? $prototype{$t->[0]} : "void") - . " $t->[1](" . - join(',', @decl_args) . ");"; - } - # add on the function reference, if supplied, to the start of - # the doc string - if ( defined $docstring ) { - $docstring = "=for ref\n\n$funcref\n\n$docstring" if defined $funcref; - } else { - $docstring = ''; - } - # IntVals - my @intcopy_code = map +($KIND2I{$_->[0]}[3]||sub{})->($_), @args2; - my $code = join("\n", @intcopy_code, "$func($funcargs);"); - print <<"ENDDBG" if $debug; -pp_def('$pname', - Pars => '$pars',@{[ !$otherpars ? '' : " - OtherPars => '$otherpars',"]}@{[ !$argorder ? '' : " - ArgOrder => [qw(@$argorder)],"]} - Code => pp_line_numbers(__LINE__, <<'EOF'), -$code -EOF - CHeader => '@cheaders', - GenericTypes => [qw(@{[join ", ", map $_->[0], @talts]})], - Doc => <<'EOF', -$docstring -EOF -); -ENDDBG - pp_def($pname, - Pars => $pars, - ($otherpars ? (OtherPars => $otherpars) : ()), - $argorder ? (ArgOrder => $argorder) : (), - Code => $code, - CHeader => join("\n", @cheaders), - GenericTypes => [map $_->[0], @talts], - Doc => $docstring, - ); -} - pp_add_macros( PCHDF => sub {my ($k, $x, $s, $out) = @_; ' /* PDL version: K, X, S are var names, 4th param output */ @@ -5687,19 +5553,25 @@ which should never happen. EOF ); -defslatec('pchip_chfe', - {D => 'dpchfe'}, - 'IntVal n=$SIZE(n); - Mat x (n); - Mat f (n); - Mat d (n); - Logical [io] skip (); - IntVal ne=$SIZE(ne); - Mat xe (ne); - Mat [o] fe (ne); - Integer [o] ierr (); - ', -'Given a piecewise cubic Hermite function - such as from +pp_def('pchip_chfe', + Pars => 'x(n); f(n); d(n); int [io] skip(); xe(ne); + [o] fe(ne); longlong [o] ierr()', + RedoDimsCode => <<'EOF', +if ($SIZE(n) < 2) $CROAK("NUMBER OF DATA POINTS LESS THAN TWO"); +if ($SIZE(ne) < 1) $CROAK("NUMBER OF EVALUATION POINTS LESS THAN ONE"); +EOF + Code => pp_line_numbers(__LINE__, <<'EOF'), +$CHFD(0); +EOF + GenericTypes => $F, + Doc => <<'EOF', +=for ref + +Evaluate a piecewise cubic Hermite function at an array of points. +May be used by itself for Hermite interpolation, or as an evaluator +for L or L. + +Given a piecewise cubic Hermite function - such as from L - evaluate the function (C<$fe>) at a set of points (C<$xe>). If derivative values are also required, use L. @@ -5767,13 +5639,7 @@ E0 if extrapolation was performed at C points (The FE-array has not been changed in any of these cases.) NOTE: The above errors are checked in the order listed, and following arguments have B been validated. - -=cut - -', -'Evaluate a piecewise cubic Hermite function at an array of points. -May be used by itself for Hermite interpolation, or as an evaluator -for L or L.' +EOF ); pp_def('pchip_chia',