File Coverage

lib/PDL/Primitive-pp-pchip_chfe.c

Criterion	Covered	Total	%
statement	53	75	70.6
branch	47	246	19.1
condition			n/a
subroutine			n/a
pod			n/a
total	100	321	31.1

#line 453 "lib/PDL/PP.pm"
/*
 * THIS FILE WAS GENERATED BY PDL::PP from lib/PDL/Primitive.pd! Do not modify!
 */

#define PDL_FREE_CODE(trans, destroy, comp_free_code, ntpriv_free_code) \
    if (destroy) { \
	comp_free_code \
    } \
    if ((trans)->dims_redone) { \
	ntpriv_free_code \
    }

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"
#include "pdl.h"
#include "pdlcore.h"
#define PDL PDL_Primitive
extern Core* PDL; /* Structure hold core C functions */
#line 23 "lib/PDL/Primitive-pp-pchip_chfe.c"
extern int pdl_srand_threads;
extern uint64_t *pdl_rand_state;
void pdl_srand(uint64_t **s, uint64_t seed, int n);
double pdl_drand(uint64_t *s);
#define PDL_MAYBE_SRAND \
  if (pdl_srand_threads < 0) \
    pdl_srand(&pdl_rand_state, PDL->pdl_seed(), PDL->online_cpus());
#define PDL_RAND_SET_OFFSET(v, thr, pdl) \
  if (v < 0) { \
    if (thr.mag_nthr >= 0) { \
      int thr_no = PDL->magic_get_thread(pdl); \
      if (thr_no < 0) return PDL->make_error_simple(PDL_EFATAL, "Invalid pdl_magic_get_thread!"); \
      v = thr_no == 0 ? thr_no : thr_no % PDL->online_cpus(); \
    } else { \
      v = 0; \
    } \
  }

#line 1857 "lib/PDL/PP.pm"
pdl_error pdl_pchip_chfe_redodims(pdl_trans *__privtrans) {
  pdl_error PDL_err = {0, NULL, 0};
#line 45 "lib/PDL/Primitive-pp-pchip_chfe.c"
  #ifndef PDL_DECLARE_PARAMS_pchip_chfe_0
  #define PDL_DECLARE_PARAMS_pchip_chfe_0(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ierr,PDL_PPSYM_PARAM_ierr,PDL_TYPE_PARAM_skip,PDL_PPSYM_PARAM_skip) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[0]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, f, (__privtrans->pdls[1]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, d, (__privtrans->pdls[2]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, xe, (__privtrans->pdls[3]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, fe, (__privtrans->pdls[4]), 0, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ierr, ierr, (__privtrans->pdls[5]), 0, PDL_PPSYM_PARAM_ierr) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_skip, skip, (__privtrans->pdls[6]), 0, PDL_PPSYM_PARAM_skip)
  #endif
  #define PDL_IF_BAD(t,f) f
  switch (__privtrans->__datatype) { /* Start generic switch */
    case PDL_F: {
      PDL_DECLARE_PARAMS_pchip_chfe_0(PDL_Float,F,PDL_Indx,N,PDL_Long,L)
        {if (__privtrans->ind_sizes[0] < 2)  return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "NUMBER OF DATA POINTS LESS THAN TWO");
        if (__privtrans->ind_sizes[1] < 1) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "NUMBER OF EVALUATION POINTS LESS THAN ONE");
        }
    } break;
    case PDL_D: {
      PDL_DECLARE_PARAMS_pchip_chfe_0(PDL_Double,D,PDL_Indx,N,PDL_Long,L)
        {if (__privtrans->ind_sizes[0] < 2)  return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "NUMBER OF DATA POINTS LESS THAN TWO");
        if (__privtrans->ind_sizes[1] < 1) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "NUMBER OF EVALUATION POINTS LESS THAN ONE");
        }
    } break;
    case PDL_LD: {
      PDL_DECLARE_PARAMS_pchip_chfe_0(PDL_LDouble,E,PDL_Indx,N,PDL_Long,L)
        {if (__privtrans->ind_sizes[0] < 2)  return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "NUMBER OF DATA POINTS LESS THAN TWO");
        if (__privtrans->ind_sizes[1] < 1) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "NUMBER OF EVALUATION POINTS LESS THAN ONE");
        }
    } break;
    default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_chfe: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
  }
  #undef PDL_IF_BAD

  PDL_RETERROR(PDL_err, PDL->redodims_default(__privtrans));
  return PDL_err;
}


#line 1857 "lib/PDL/PP.pm"
pdl_error pdl_pchip_chfe_readdata(pdl_trans *__privtrans) {
  pdl_error PDL_err = {0, NULL, 0};
#line 88 "lib/PDL/Primitive-pp-pchip_chfe.c"
  register PDL_Indx __n_size = __privtrans->ind_sizes[0];
  register PDL_Indx __ne_size = __privtrans->ind_sizes[1];
  if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "broadcast.incs NULL");
  /* broadcastloop declarations */
  int __brcloopval;
  register PDL_Indx __tind0,__tind1; /* counters along dim */
  register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
  /* dims here are how many steps along those dims */
  register PDL_Indx __tinc0_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
  register PDL_Indx __tinc0_f = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
  register PDL_Indx __tinc0_d = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
  register PDL_Indx __tinc0_xe = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
  register PDL_Indx __tinc0_fe = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
  register PDL_Indx __tinc0_ierr = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
  register PDL_Indx __tinc0_skip = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,0);
  register PDL_Indx __tinc1_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
  register PDL_Indx __tinc1_f = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
  register PDL_Indx __tinc1_d = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
  register PDL_Indx __tinc1_xe = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
  register PDL_Indx __tinc1_fe = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
  register PDL_Indx __tinc1_ierr = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
  register PDL_Indx __tinc1_skip = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,1);
  #define PDL_BROADCASTLOOP_START_pchip_chfe_readdata PDL_BROADCASTLOOP_START( \
    readdata, \
    __privtrans->broadcast, \
    __privtrans->vtable, \
    x_datap += __offsp[0]; \
    f_datap += __offsp[1]; \
    d_datap += __offsp[2]; \
    xe_datap += __offsp[3]; \
    fe_datap += __offsp[4]; \
    ierr_datap += __offsp[5]; \
    skip_datap += __offsp[6]; \
    , \
    (  ,x_datap += __tinc1_x - __tinc0_x * __tdims0 \
    ,f_datap += __tinc1_f - __tinc0_f * __tdims0 \
    ,d_datap += __tinc1_d - __tinc0_d * __tdims0 \
    ,xe_datap += __tinc1_xe - __tinc0_xe * __tdims0 \
    ,fe_datap += __tinc1_fe - __tinc0_fe * __tdims0 \
    ,ierr_datap += __tinc1_ierr - __tinc0_ierr * __tdims0 \
    ,skip_datap += __tinc1_skip - __tinc0_skip * __tdims0 \
    ), \
    (  ,x_datap += __tinc0_x \
    ,f_datap += __tinc0_f \
    ,d_datap += __tinc0_d \
    ,xe_datap += __tinc0_xe \
    ,fe_datap += __tinc0_fe \
    ,ierr_datap += __tinc0_ierr \
    ,skip_datap += __tinc0_skip \
    ) \
  )
  #define PDL_BROADCASTLOOP_END_pchip_chfe_readdata PDL_BROADCASTLOOP_END( \
    __privtrans->broadcast, \
    x_datap -= __tinc1_x * __tdims1 + __offsp[0]; \
    f_datap -= __tinc1_f * __tdims1 + __offsp[1]; \
    d_datap -= __tinc1_d * __tdims1 + __offsp[2]; \
    xe_datap -= __tinc1_xe * __tdims1 + __offsp[3]; \
    fe_datap -= __tinc1_fe * __tdims1 + __offsp[4]; \
    ierr_datap -= __tinc1_ierr * __tdims1 + __offsp[5]; \
    skip_datap -= __tinc1_skip * __tdims1 + __offsp[6]; \
  )
  register PDL_Indx __inc_d_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_d_n;
  register PDL_Indx __inc_f_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_f_n;
  register PDL_Indx __inc_fe_ne = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,0)]; (void)__inc_fe_ne;
  register PDL_Indx __inc_x_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_x_n;
  register PDL_Indx __inc_xe_ne = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_xe_ne;
  #ifndef PDL_DECLARE_PARAMS_pchip_chfe_1
  #define PDL_DECLARE_PARAMS_pchip_chfe_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ierr,PDL_PPSYM_PARAM_ierr,PDL_TYPE_PARAM_skip,PDL_PPSYM_PARAM_skip) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, f, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, d, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, xe, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_OP, fe, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ierr, ierr, (__privtrans->pdls[5]), 1, PDL_PPSYM_PARAM_ierr) \
    PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_skip, skip, (__privtrans->pdls[6]), 1, PDL_PPSYM_PARAM_skip)
  #endif
  #define PDL_IF_BAD(t,f) f
  switch (__privtrans->__datatype) { /* Start generic switch */
    case PDL_F: {
      PDL_DECLARE_PARAMS_pchip_chfe_1(PDL_Float,F,PDL_Indx,N,PDL_Long,L)
      PDL_BROADCASTLOOP_START_pchip_chfe_readdata    {
          #line 5591 "lib/PDL/Primitive.pd"

          /*  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. */
          /*  VALIDITY-CHECK ARGUMENTS. */
          if (!(skip_datap)[0]) {
            {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
                if ((x_datap)[0+(__inc_x_n*(n))] > (x_datap)[0+(__inc_x_n*(n-1))]) continue;
                (ierr_datap)[0] = -3;
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "X-ARRAY NOT STRICTLY INCREASING");
              }} /* Close n=1 */
          }
          /*  FUNCTION DEFINITION IS OK, GO ON. */
          (ierr_datap)[0] = 0;
          (skip_datap)[0] = 1;
          /*  LOOP OVER INTERVALS.    (   INTERVAL INDEX IS  IL = IR-1  . ) */
          /*                ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) */
          PDL_Indx n = __privtrans->ind_sizes[0], ne = __privtrans->ind_sizes[1];
          PDL_Indx jfirst = 0, ir;
          for (ir = 1; ir < n && jfirst < ne; ++ir) {
          /*   SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. */
          /*   LOCATE ALL POINTS IN INTERVAL. */
            char located = 0;
            PDL_Indx j = jfirst;
            {/* Open ne=jfirst */ PDL_EXPAND2(register PDL_Indx ne=PDLMAX((jfirst),0), __ne_stop=(__ne_size)); for(; ne<__ne_stop; ne+=1) {
                j = ne;
                if ((xe_datap)[0+(__inc_xe_ne*(ne))] >= (x_datap)[0+(__inc_x_n*(ir))]) {
                  located = 1;
                  break;
                }
              }} /* Close ne=jfirst */
            if (!located || ir == n-1)
              j = ne;
          /*   HAVE LOCATED FIRST POINT BEYOND INTERVAL. */
            PDL_Indx nj = j - jfirst;
          /*   SKIP EVALUATION IF NO POINTS IN INTERVAL. */
            if (nj == 0)
              continue;
          /*   EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 . */
          /*     ---------------------------------------------------------------- */
            PDL_Indx next[] = {0,0};
            do { /* inline dchfdv */
          /* Local variables */
              PDL_Float x1 = (x_datap)[0+(__inc_x_n*(ir-1))], x2 = (x_datap)[0+(__inc_x_n*(ir))];
              PDL_Float f1 = (f_datap)[0+(__inc_f_n*(ir-1))], f2 = (f_datap)[0+(__inc_f_n*(ir))];
              PDL_Float d1 = (d_datap)[0+(__inc_d_n*(ir-1))], d2 = (d_datap)[0+(__inc_d_n*(ir))];
              PDL_Float h = x2 - x1;
              if (h == 0.)
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "INTERVAL ENDPOINTS EQUAL");
          /*  INITIALIZE. */
              PDL_Float xmi = PDLMIN(0.,h);
              PDL_Float xma = PDLMAX(0.,h);
          /*  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). */
              PDL_Float delta = (f2 - f1) / h;
              PDL_Float del1 = (d1 - delta) / h;
              PDL_Float del2 = (d2 - delta) / h;
          /*                       (DELTA IS NO LONGER NEEDED.) */
              PDL_Float c2 = -(del1 + del1 + del2);

              PDL_Float c3 = (del1 + del2) / h;
          /*                 (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) */

          /*  EVALUATION LOOP. */
              {/* Open ne=:nj */ PDL_EXPAND2(register PDL_Indx ne=0, __ne_stop=PDLMIN(nj, (__ne_size))); for(; ne<__ne_stop; ne+=1) {
                  PDL_Float x = (xe_datap)[0+(__inc_xe_ne*(ne+jfirst))] - x1;
                  (fe_datap)[0+(__inc_fe_ne*(ne+jfirst))] = f1 + x * (d1 + x * (c2 + x * c3));

            /*      COUNT EXTRAPOLATION POINTS. */
                  if (x < xmi)
                    ++next[0];
                  if (x > xma)
                    ++next[1];
            /*    (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) */
                }} /* Close ne=:nj */
            } while (0); /* end inline dchfdv */
          /*     ---------------------------------------------------------------- */
            if (next[1] != 0) {
          /*       IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE */
          /*       RIGHT OF X(IR). */
          /*        THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
              (ierr_datap)[0] += 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_datap)[0] += next[0];
                jfirst = j;
                continue;
              }
          /*        XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST */
          /*        EVALUATION INTERVAL. */
          /*        FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). */
              located = 0;
              PDL_Indx i = jfirst;
              {/* Open ne=jfirst:j */ PDL_EXPAND2(register PDL_Indx ne=PDLMAX((jfirst),0), __ne_stop=PDLMIN(j, (__ne_size))); for(; ne<__ne_stop; ne+=1) {
                  i = ne;
                  if ((xe_datap)[0+(__inc_xe_ne*(ne))] < (x_datap)[0+(__inc_x_n*(ir-1))]) {
                    located = 1;
                    break;
                  }
                }} /* Close ne=jfirst:j */
              if (!located) {
          /*        NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR */
          /*           IN DCHFDV. */
                (ierr_datap)[0] = -5;
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "ERROR RETURN FROM DCHFDV -- FATAL");
              }
          /*        RESET J.  (THIS WILL BE THE NEW JFIRST.) */
              j = i;
          /*        NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. */
              {/* Open n=:ir */ PDL_EXPAND2(register PDL_Indx n=0, __n_stop=PDLMIN(ir, (__n_size))); for(; n<__n_stop; n+=1) {
                  i = n;
                  if ((xe_datap)[0+(__inc_xe_ne*(j))] < (x_datap)[0+(__inc_x_n*(n))])
                    break;
                }} /* Close n=:ir */
          /*        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 = PDLMAX(0,i-1);
            }
            jfirst = j;
          /*   END OF IR-LOOP. */
          }
          ;
          #line 310 "lib/PDL/Primitive-pp-pchip_chfe.c"
          }PDL_BROADCASTLOOP_END_pchip_chfe_readdata
    } break;
    case PDL_D: {
      PDL_DECLARE_PARAMS_pchip_chfe_1(PDL_Double,D,PDL_Indx,N,PDL_Long,L)
      PDL_BROADCASTLOOP_START_pchip_chfe_readdata    {
          #line 5591 "lib/PDL/Primitive.pd"

          /*  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. */
          /*  VALIDITY-CHECK ARGUMENTS. */
          if (!(skip_datap)[0]) {
            {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
                if ((x_datap)[0+(__inc_x_n*(n))] > (x_datap)[0+(__inc_x_n*(n-1))]) continue;
                (ierr_datap)[0] = -3;
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "X-ARRAY NOT STRICTLY INCREASING");
              }} /* Close n=1 */
          }
          /*  FUNCTION DEFINITION IS OK, GO ON. */
          (ierr_datap)[0] = 0;
          (skip_datap)[0] = 1;
          /*  LOOP OVER INTERVALS.    (   INTERVAL INDEX IS  IL = IR-1  . ) */
          /*                ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) */
          PDL_Indx n = __privtrans->ind_sizes[0], ne = __privtrans->ind_sizes[1];
          PDL_Indx jfirst = 0, ir;
          for (ir = 1; ir < n && jfirst < ne; ++ir) {
          /*   SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. */
          /*   LOCATE ALL POINTS IN INTERVAL. */
            char located = 0;
            PDL_Indx j = jfirst;
            {/* Open ne=jfirst */ PDL_EXPAND2(register PDL_Indx ne=PDLMAX((jfirst),0), __ne_stop=(__ne_size)); for(; ne<__ne_stop; ne+=1) {
                j = ne;
                if ((xe_datap)[0+(__inc_xe_ne*(ne))] >= (x_datap)[0+(__inc_x_n*(ir))]) {
                  located = 1;
                  break;
                }
              }} /* Close ne=jfirst */
            if (!located || ir == n-1)
              j = ne;
          /*   HAVE LOCATED FIRST POINT BEYOND INTERVAL. */
            PDL_Indx nj = j - jfirst;
          /*   SKIP EVALUATION IF NO POINTS IN INTERVAL. */
            if (nj == 0)
              continue;
          /*   EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 . */
          /*     ---------------------------------------------------------------- */
            PDL_Indx next[] = {0,0};
            do { /* inline dchfdv */
          /* Local variables */
              PDL_Double x1 = (x_datap)[0+(__inc_x_n*(ir-1))], x2 = (x_datap)[0+(__inc_x_n*(ir))];
              PDL_Double f1 = (f_datap)[0+(__inc_f_n*(ir-1))], f2 = (f_datap)[0+(__inc_f_n*(ir))];
              PDL_Double d1 = (d_datap)[0+(__inc_d_n*(ir-1))], d2 = (d_datap)[0+(__inc_d_n*(ir))];
              PDL_Double h = x2 - x1;
              if (h == 0.)
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "INTERVAL ENDPOINTS EQUAL");
          /*  INITIALIZE. */
              PDL_Double xmi = PDLMIN(0.,h);
              PDL_Double xma = PDLMAX(0.,h);
          /*  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). */
              PDL_Double delta = (f2 - f1) / h;
              PDL_Double del1 = (d1 - delta) / h;
              PDL_Double del2 = (d2 - delta) / h;
          /*                       (DELTA IS NO LONGER NEEDED.) */
              PDL_Double c2 = -(del1 + del1 + del2);

              PDL_Double c3 = (del1 + del2) / h;
          /*                 (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) */

          /*  EVALUATION LOOP. */
              {/* Open ne=:nj */ PDL_EXPAND2(register PDL_Indx ne=0, __ne_stop=PDLMIN(nj, (__ne_size))); for(; ne<__ne_stop; ne+=1) {
                  PDL_Double x = (xe_datap)[0+(__inc_xe_ne*(ne+jfirst))] - x1;
                  (fe_datap)[0+(__inc_fe_ne*(ne+jfirst))] = f1 + x * (d1 + x * (c2 + x * c3));

            /*      COUNT EXTRAPOLATION POINTS. */
                  if (x < xmi)
                    ++next[0];
                  if (x > xma)
                    ++next[1];
            /*    (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) */
                }} /* Close ne=:nj */
            } while (0); /* end inline dchfdv */
          /*     ---------------------------------------------------------------- */
            if (next[1] != 0) {
          /*       IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE */
          /*       RIGHT OF X(IR). */
          /*        THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
              (ierr_datap)[0] += 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_datap)[0] += next[0];
                jfirst = j;
                continue;
              }
          /*        XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST */
          /*        EVALUATION INTERVAL. */
          /*        FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). */
              located = 0;
              PDL_Indx i = jfirst;
              {/* Open ne=jfirst:j */ PDL_EXPAND2(register PDL_Indx ne=PDLMAX((jfirst),0), __ne_stop=PDLMIN(j, (__ne_size))); for(; ne<__ne_stop; ne+=1) {
                  i = ne;
                  if ((xe_datap)[0+(__inc_xe_ne*(ne))] < (x_datap)[0+(__inc_x_n*(ir-1))]) {
                    located = 1;
                    break;
                  }
                }} /* Close ne=jfirst:j */
              if (!located) {
          /*        NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR */
          /*           IN DCHFDV. */
                (ierr_datap)[0] = -5;
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "ERROR RETURN FROM DCHFDV -- FATAL");
              }
          /*        RESET J.  (THIS WILL BE THE NEW JFIRST.) */
              j = i;
          /*        NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. */
              {/* Open n=:ir */ PDL_EXPAND2(register PDL_Indx n=0, __n_stop=PDLMIN(ir, (__n_size))); for(; n<__n_stop; n+=1) {
                  i = n;
                  if ((xe_datap)[0+(__inc_xe_ne*(j))] < (x_datap)[0+(__inc_x_n*(n))])
                    break;
                }} /* Close n=:ir */
          /*        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 = PDLMAX(0,i-1);
            }
            jfirst = j;
          /*   END OF IR-LOOP. */
          }
          ;
          #line 456 "lib/PDL/Primitive-pp-pchip_chfe.c"
          }PDL_BROADCASTLOOP_END_pchip_chfe_readdata
    } break;
    case PDL_LD: {
      PDL_DECLARE_PARAMS_pchip_chfe_1(PDL_LDouble,E,PDL_Indx,N,PDL_Long,L)
      PDL_BROADCASTLOOP_START_pchip_chfe_readdata    {
          #line 5591 "lib/PDL/Primitive.pd"

          /*  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. */
          /*  VALIDITY-CHECK ARGUMENTS. */
          if (!(skip_datap)[0]) {
            {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
                if ((x_datap)[0+(__inc_x_n*(n))] > (x_datap)[0+(__inc_x_n*(n-1))]) continue;
                (ierr_datap)[0] = -3;
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "X-ARRAY NOT STRICTLY INCREASING");
              }} /* Close n=1 */
          }
          /*  FUNCTION DEFINITION IS OK, GO ON. */
          (ierr_datap)[0] = 0;
          (skip_datap)[0] = 1;
          /*  LOOP OVER INTERVALS.    (   INTERVAL INDEX IS  IL = IR-1  . ) */
          /*                ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) */
          PDL_Indx n = __privtrans->ind_sizes[0], ne = __privtrans->ind_sizes[1];
          PDL_Indx jfirst = 0, ir;
          for (ir = 1; ir < n && jfirst < ne; ++ir) {
          /*   SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. */
          /*   LOCATE ALL POINTS IN INTERVAL. */
            char located = 0;
            PDL_Indx j = jfirst;
            {/* Open ne=jfirst */ PDL_EXPAND2(register PDL_Indx ne=PDLMAX((jfirst),0), __ne_stop=(__ne_size)); for(; ne<__ne_stop; ne+=1) {
                j = ne;
                if ((xe_datap)[0+(__inc_xe_ne*(ne))] >= (x_datap)[0+(__inc_x_n*(ir))]) {
                  located = 1;
                  break;
                }
              }} /* Close ne=jfirst */
            if (!located || ir == n-1)
              j = ne;
          /*   HAVE LOCATED FIRST POINT BEYOND INTERVAL. */
            PDL_Indx nj = j - jfirst;
          /*   SKIP EVALUATION IF NO POINTS IN INTERVAL. */
            if (nj == 0)
              continue;
          /*   EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 . */
          /*     ---------------------------------------------------------------- */
            PDL_Indx next[] = {0,0};
            do { /* inline dchfdv */
          /* Local variables */
              PDL_LDouble x1 = (x_datap)[0+(__inc_x_n*(ir-1))], x2 = (x_datap)[0+(__inc_x_n*(ir))];
              PDL_LDouble f1 = (f_datap)[0+(__inc_f_n*(ir-1))], f2 = (f_datap)[0+(__inc_f_n*(ir))];
              PDL_LDouble d1 = (d_datap)[0+(__inc_d_n*(ir-1))], d2 = (d_datap)[0+(__inc_d_n*(ir))];
              PDL_LDouble h = x2 - x1;
              if (h == 0.)
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "INTERVAL ENDPOINTS EQUAL");
          /*  INITIALIZE. */
              PDL_LDouble xmi = PDLMIN(0.,h);
              PDL_LDouble xma = PDLMAX(0.,h);
          /*  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). */
              PDL_LDouble delta = (f2 - f1) / h;
              PDL_LDouble del1 = (d1 - delta) / h;
              PDL_LDouble del2 = (d2 - delta) / h;
          /*                       (DELTA IS NO LONGER NEEDED.) */
              PDL_LDouble c2 = -(del1 + del1 + del2);

              PDL_LDouble c3 = (del1 + del2) / h;
          /*                 (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) */

          /*  EVALUATION LOOP. */
              {/* Open ne=:nj */ PDL_EXPAND2(register PDL_Indx ne=0, __ne_stop=PDLMIN(nj, (__ne_size))); for(; ne<__ne_stop; ne+=1) {
                  PDL_LDouble x = (xe_datap)[0+(__inc_xe_ne*(ne+jfirst))] - x1;
                  (fe_datap)[0+(__inc_fe_ne*(ne+jfirst))] = f1 + x * (d1 + x * (c2 + x * c3));

            /*      COUNT EXTRAPOLATION POINTS. */
                  if (x < xmi)
                    ++next[0];
                  if (x > xma)
                    ++next[1];
            /*    (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) */
                }} /* Close ne=:nj */
            } while (0); /* end inline dchfdv */
          /*     ---------------------------------------------------------------- */
            if (next[1] != 0) {
          /*       IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE */
          /*       RIGHT OF X(IR). */
          /*        THESE ARE ACTUALLY EXTRAPOLATION POINTS. */
              (ierr_datap)[0] += 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_datap)[0] += next[0];
                jfirst = j;
                continue;
              }
          /*        XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST */
          /*        EVALUATION INTERVAL. */
          /*        FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). */
              located = 0;
              PDL_Indx i = jfirst;
              {/* Open ne=jfirst:j */ PDL_EXPAND2(register PDL_Indx ne=PDLMAX((jfirst),0), __ne_stop=PDLMIN(j, (__ne_size))); for(; ne<__ne_stop; ne+=1) {
                  i = ne;
                  if ((xe_datap)[0+(__inc_xe_ne*(ne))] < (x_datap)[0+(__inc_x_n*(ir-1))]) {
                    located = 1;
                    break;
                  }
                }} /* Close ne=jfirst:j */
              if (!located) {
          /*        NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR */
          /*           IN DCHFDV. */
                (ierr_datap)[0] = -5;
                return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chfe:" "ERROR RETURN FROM DCHFDV -- FATAL");
              }
          /*        RESET J.  (THIS WILL BE THE NEW JFIRST.) */
              j = i;
          /*        NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. */
              {/* Open n=:ir */ PDL_EXPAND2(register PDL_Indx n=0, __n_stop=PDLMIN(ir, (__n_size))); for(; n<__n_stop; n+=1) {
                  i = n;
                  if ((xe_datap)[0+(__inc_xe_ne*(j))] < (x_datap)[0+(__inc_x_n*(n))])
                    break;
                }} /* Close n=:ir */
          /*        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 = PDLMAX(0,i-1);
            }
            jfirst = j;
          /*   END OF IR-LOOP. */
          }
          ;
          #line 602 "lib/PDL/Primitive-pp-pchip_chfe.c"
          }PDL_BROADCASTLOOP_END_pchip_chfe_readdata
    } break;
    default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_chfe: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
  }
  #undef PDL_IF_BAD
  return PDL_err;
}

static pdl_datatypes pdl_pchip_chfe_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
static PDL_Indx pdl_pchip_chfe_vtable_realdims[] = { 1, 1, 1, 1, 1, 0, 0 };
static char *pdl_pchip_chfe_vtable_parnames[] = { "x","f","d","xe","fe","ierr","skip" };
static short pdl_pchip_chfe_vtable_parflags[] = {
  0,
  0,
  0,
  0,
  PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
  PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE,
  PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE
};
static pdl_datatypes pdl_pchip_chfe_vtable_partypes[] = { -1, -1, -1, -1, -1, PDL_IND, PDL_L };
static PDL_Indx pdl_pchip_chfe_vtable_realdims_starts[] = { 0, 1, 2, 3, 4, 5, 5 };
static PDL_Indx pdl_pchip_chfe_vtable_realdims_ind_ids[] = { 0, 0, 0, 1, 1 };
static char *pdl_pchip_chfe_vtable_indnames[] = { "n","ne" };
pdl_transvtable pdl_pchip_chfe_vtable = {
  PDL_TRANS_DO_BROADCAST, 0, pdl_pchip_chfe_vtable_gentypes, 4, 7, NULL /*CORE21*/,
  pdl_pchip_chfe_vtable_realdims, pdl_pchip_chfe_vtable_parnames,
  pdl_pchip_chfe_vtable_parflags, pdl_pchip_chfe_vtable_partypes,
  pdl_pchip_chfe_vtable_realdims_starts, pdl_pchip_chfe_vtable_realdims_ind_ids, 5,
  2, pdl_pchip_chfe_vtable_indnames,
  pdl_pchip_chfe_redodims, pdl_pchip_chfe_readdata, NULL,
  NULL,
  0,"PDL::Primitive::pchip_chfe"
};


pdl_error pdl_run_pchip_chfe(pdl  *x,pdl  *f,pdl  *d,pdl  *xe,pdl  *fe,pdl  *ierr,pdl  *skip) {
  pdl_error PDL_err = {0, NULL, 0};
  if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
  pdl_trans *__privtrans = PDL->create_trans(&pdl_pchip_chfe_vtable);
  if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
  __privtrans->pdls[0] = x;
  __privtrans->pdls[1] = f;
  __privtrans->pdls[2] = d;
  __privtrans->pdls[3] = xe;
  __privtrans->pdls[4] = fe;
  __privtrans->pdls[5] = ierr;
  __privtrans->pdls[6] = skip;
  PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
  PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
  return PDL_err;
}