Slatec PCHIP routines [io] skip Pars, fix doc to match SLATEC code - …

…#468
PDLPorters · Apr 1, 2024 · 9469b35 · 9469b35
1 parent b0365c8
commit 9469b35
Show file tree

Hide file tree

Showing 4 changed files with 19 additions and 14 deletions.
diff --git a/Changes b/Changes
@@ -1,3 +1,5 @@
+- Slatec PCHIP routines now have idiomatic [io] `skip` Pars, with doc fixed to match SLATEC code (#468)
+
 2.086 2024-03-31
 - fix 64-bit problem in Minuit revealed by LTO (#468) - thanks @eli-schwartz for report
 - add pdl.ntrans_children to speed up destroying ndarrays (#421,#451)

diff --git a/Libtmp/Func.pm b/Libtmp/Func.pm
@@ -738,7 +738,13 @@ sub _interp_hermite {
     # get gradient
     my $g = $self->_get_value( 'g' );
 
-    my ( $yi, $ierr ) = chfe( $x, $y, $g, 0, $xi );
+    my @highest_dims;
+    for my $param ($x,$y,$g,$xi) {
+      my @this_dims = $param->dims;
+      @highest_dims = @this_dims if @this_dims > @highest_dims;
+    }
+    shift @highest_dims; # don't care about bottom one
+    my ( $yi, $ierr ) = chfe( $x, $y, $g, PDL->ones(PDL::long(),@highest_dims), $xi );
     $self->{flags}{routine} = "chfe";
     $self->_set_value( err => $ierr );
 

diff --git a/Libtmp/Slatec/slatec.pd b/Libtmp/Slatec/slatec.pd
@@ -295,11 +295,8 @@ sub intcopy {
 #       (this allows the PCHIP routines to use 2D data, but as it's
 #        done in FORTRAN array order, and PDL has a much richer way
 #        of accessing parts of an array we force the data to be 1D).
-#   IntCopy
-#     - the PCHIP routines may change the value from 0 to 1 if an
-#       error occurs but the checks were successful. As this complicates
-#       things we copy the user value to a temporary variable,
-#       so that the sent in value is not changed.
+#   Logical
+#     - A Fortran Boolean, defined as same size as REAL/INTEGER so 'long'
 #   FortranIndex
 #     - pchid()/dpchid() require FORTRAN array indices, so 
 #       this type flags that we should add 1 onto the input values
@@ -319,7 +316,7 @@ my %KIND2I = (
     # convert from C to F77 index so add 1
   FortranIndex => ['longlong ', intcopy(sub {$_[0][2]}, sub {"\$$_[0][2]() + 1"})],
   Integer => ['longlong ', sub{"\$P($_[0][2])"}, sub {"PDL_LongLong *"}],
-  IntCopy => ['longlong ', intcopy(sub {$_[0][2]}, sub {"\$$_[0][2]()"})],
+  Logical => ['int ', sub{"\$P($_[0][2])"}, sub {"PDL_Long *"}],
   Dim => [undef, intcopy(sub {"dim_$_[0][2]"}, sub {"\$SIZE($_[0][2])"})],
   MDim => [undef, intcopy(sub {$_[0][2]}, sub {$_[1]{$_[0][2]}})],
   Incfd => [undef, intcopy(sub {$_[0][2]}, sub {1})],
@@ -1085,7 +1082,7 @@ defslatec(
            Mat          f       (n);
            Mat          d       (n);
            Incfd        incfd;
-           IntCopy      skip    ();
+           Logical [io] skip    ();
            Dim          ne;
            Mat          xe      (ne);
            Mat     [o]  fe      (ne);
@@ -1143,7 +1140,7 @@ defslatec(
            Mat          f       (n);
            Mat          d       (n);
            Incfd        incfd;
-           IntCopy      skip    ();
+           Logical [io] skip    ();
            Dim          ne;
            Mat          xe      (ne);
            Mat     [o]  fe      (ne);
@@ -1195,7 +1192,7 @@ defslatec(
            Mat          f       (n);
            Mat          d       (n);
            Incfd        incfd;
-           IntCopy      skip    ();
+           Logical [io] skip    ();
            Mat          la    ();
            Mat          lb    ();
            FuncRet [o]  ans     ();
@@ -1260,7 +1257,7 @@ defslatec(
            Mat          f       (n);
            Mat          d       (n);
            Incfd        incfd;
-           IntCopy      skip    ();
+           Logical [io] skip    ();
            FortranIndex ia      ();
            FortranIndex ib      ();
            FuncRet [o]  ans     ();
@@ -1313,7 +1310,7 @@ defslatec(
            Mat          f       (n);
            Mat          d       (n);
            Incfd        incfd;
-           IntCopy      skip    ();
+           Logical [io] skip    ();
            Integer [o]  ismon   (n);
            Integer [o]  ierr    ();
           ',

diff --git a/Libtmp/Slatec/t/slatec.t b/Libtmp/Slatec/t/slatec.t
@@ -210,14 +210,14 @@ $f = $x * $x;
 ( $d, $err ) = chim($x, $f);
 
 $ans = pdl( 9.0**3, (8.0**3-1.0**3) ) / 3.0;
-( my $int, $err ) = chia($x, $f, $d, 1, pdl(0.0,1.0), pdl(9.0,8.0));
+( my $int, $err ) = chia($x, $f, $d, my $skip=zeroes(2), pdl(0.0,1.0), pdl(9.0,8.0));
 ok(all($err == 0));
 ok(all( abs($int-$ans) < 0.04 ) );
 
 my $hi = pdl( $x->at(9), $x->at(7) );
 my $lo = pdl( $x->at(0), $x->at(1) );
 $ans = ($hi**3 - $lo**3) / 3;
-( $int, $err ) = chid( $x, $f, $d, 1, pdl(0,1), pdl(9,7) );
+( $int, $err ) = chid( $x, $f, $d, $skip=zeroes(2), pdl(0,1), pdl(9,7) );
 ok(all($err == 0));
 ok(all( abs($int-$ans) < 0.06 ) );
 ## print STDERR "int=$int; ans=$ans; int-ans=".($int-$ans)."\n";