File Coverage

lib/PDL/Primitive-pp-pchip_chsp.c
Criterion Covered Total %
statement 55 75 73.3
branch 50 264 18.9
condition n/a
subroutine n/a
pod n/a
total 105 339 30.9


line stmt bran cond sub pod time code
1              
2             #line 453 "lib/PDL/PP.pm"
3             /*
4             * THIS FILE WAS GENERATED BY PDL::PP from lib/PDL/Primitive.pd! Do not modify!
5             */
6              
7             #define PDL_FREE_CODE(trans, destroy, comp_free_code, ntpriv_free_code) \
8             if (destroy) { \
9             comp_free_code \
10             } \
11             if ((trans)->dims_redone) { \
12             ntpriv_free_code \
13             }
14              
15             #include "EXTERN.h"
16             #include "perl.h"
17             #include "XSUB.h"
18             #include "pdl.h"
19             #include "pdlcore.h"
20             #define PDL PDL_Primitive
21             extern Core* PDL; /* Structure hold core C functions */
22             #line 23 "lib/PDL/Primitive-pp-pchip_chsp.c"
23             extern int pdl_srand_threads;
24             extern uint64_t *pdl_rand_state;
25             void pdl_srand(uint64_t **s, uint64_t seed, int n);
26             double pdl_drand(uint64_t *s);
27             #define PDL_MAYBE_SRAND \
28             if (pdl_srand_threads < 0) \
29             pdl_srand(&pdl_rand_state, PDL->pdl_seed(), PDL->online_cpus());
30             #define PDL_RAND_SET_OFFSET(v, thr, pdl) \
31             if (v < 0) { \
32             if (thr.mag_nthr >= 0) { \
33             int thr_no = PDL->magic_get_thread(pdl); \
34             if (thr_no < 0) return PDL->make_error_simple(PDL_EFATAL, "Invalid pdl_magic_get_thread!"); \
35             v = thr_no == 0 ? thr_no : thr_no % PDL->online_cpus(); \
36             } else { \
37             v = 0; \
38             } \
39             }
40              
41             #line 1857 "lib/PDL/PP.pm"
42             pdl_error pdl_pchip_chsp_redodims(pdl_trans *__privtrans) {
43             pdl_error PDL_err = {0, NULL, 0};
44             #line 45 "lib/PDL/Primitive-pp-pchip_chsp.c"
45 3           __privtrans->ind_sizes[1] = 2;
46             #ifndef PDL_DECLARE_PARAMS_pchip_chsp_0
47             #define PDL_DECLARE_PARAMS_pchip_chsp_0(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ic,PDL_PPSYM_PARAM_ic,PDL_TYPE_PARAM_ierr,PDL_PPSYM_PARAM_ierr) \
48             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ic, ic, (__privtrans->pdls[0]), 0, PDL_PPSYM_PARAM_ic) \
49             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, vc, (__privtrans->pdls[1]), 0, PDL_PPSYM_OP) \
50             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[2]), 0, PDL_PPSYM_OP) \
51             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, f, (__privtrans->pdls[3]), 0, PDL_PPSYM_OP) \
52             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, d, (__privtrans->pdls[4]), 0, PDL_PPSYM_OP) \
53             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ierr, ierr, (__privtrans->pdls[5]), 0, PDL_PPSYM_PARAM_ierr) \
54             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, dx, (__privtrans->pdls[6]), 0, PDL_PPSYM_OP) \
55             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, dy_dx, (__privtrans->pdls[7]), 0, PDL_PPSYM_OP)
56             #endif
57             #define PDL_IF_BAD(t,f) f
58 3           switch (__privtrans->__datatype) { /* Start generic switch */
59 0           case PDL_F: {
60 0 0         PDL_DECLARE_PARAMS_pchip_chsp_0(PDL_Float,F,PDL_SByte,A,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
61 0 0         {if (__privtrans->ind_sizes[0] < 2) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "NUMBER OF DATA POINTS LESS THAN TWO");
62             }
63 0           } break;
64 3           case PDL_D: {
65 3 50         PDL_DECLARE_PARAMS_pchip_chsp_0(PDL_Double,D,PDL_SByte,A,PDL_Indx,N)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
66 3 50         {if (__privtrans->ind_sizes[0] < 2) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "NUMBER OF DATA POINTS LESS THAN TWO");
67             }
68 3           } break;
69 0           case PDL_LD: {
70 0 0         PDL_DECLARE_PARAMS_pchip_chsp_0(PDL_LDouble,E,PDL_SByte,A,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
71 0 0         {if (__privtrans->ind_sizes[0] < 2) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "NUMBER OF DATA POINTS LESS THAN TWO");
72             }
73 0           } break;
74 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_chsp: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
75             }
76             #undef PDL_IF_BAD
77              
78 3 50         PDL_RETERROR(PDL_err, PDL->redodims_default(__privtrans));
79 3           return PDL_err;
80             }
81              
82              
83             #line 1857 "lib/PDL/PP.pm"
84             pdl_error pdl_pchip_chsp_readdata(pdl_trans *__privtrans) {
85             pdl_error PDL_err = {0, NULL, 0};
86             #line 87 "lib/PDL/Primitive-pp-pchip_chsp.c"
87 3           register PDL_Indx __n_size = __privtrans->ind_sizes[0];
88 3 50         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "broadcast.incs NULL");
89             /* broadcastloop declarations */
90             int __brcloopval;
91             register PDL_Indx __tind0,__tind1; /* counters along dim */
92 3           register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
93             /* dims here are how many steps along those dims */
94 3           register PDL_Indx __tinc0_ic = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
95 3           register PDL_Indx __tinc0_vc = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
96 3           register PDL_Indx __tinc0_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
97 3           register PDL_Indx __tinc0_f = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
98 3           register PDL_Indx __tinc0_d = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
99 3           register PDL_Indx __tinc0_ierr = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
100 3           register PDL_Indx __tinc0_dx = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,0);
101 3           register PDL_Indx __tinc0_dy_dx = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,7,0);
102 3           register PDL_Indx __tinc1_ic = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
103 3           register PDL_Indx __tinc1_vc = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
104 3           register PDL_Indx __tinc1_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
105 3           register PDL_Indx __tinc1_f = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
106 3           register PDL_Indx __tinc1_d = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
107 3           register PDL_Indx __tinc1_ierr = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
108 3           register PDL_Indx __tinc1_dx = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,1);
109 3           register PDL_Indx __tinc1_dy_dx = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,7,1);
110             #define PDL_BROADCASTLOOP_START_pchip_chsp_readdata PDL_BROADCASTLOOP_START( \
111             readdata, \
112             __privtrans->broadcast, \
113             __privtrans->vtable, \
114             ic_datap += __offsp[0]; \
115             vc_datap += __offsp[1]; \
116             x_datap += __offsp[2]; \
117             f_datap += __offsp[3]; \
118             d_datap += __offsp[4]; \
119             ierr_datap += __offsp[5]; \
120             dx_datap += __offsp[6]; \
121             dy_dx_datap += __offsp[7]; \
122             , \
123             ( ,ic_datap += __tinc1_ic - __tinc0_ic * __tdims0 \
124             ,vc_datap += __tinc1_vc - __tinc0_vc * __tdims0 \
125             ,x_datap += __tinc1_x - __tinc0_x * __tdims0 \
126             ,f_datap += __tinc1_f - __tinc0_f * __tdims0 \
127             ,d_datap += __tinc1_d - __tinc0_d * __tdims0 \
128             ,ierr_datap += __tinc1_ierr - __tinc0_ierr * __tdims0 \
129             ,dx_datap += __tinc1_dx - __tinc0_dx * __tdims0 \
130             ,dy_dx_datap += __tinc1_dy_dx - __tinc0_dy_dx * __tdims0 \
131             ), \
132             ( ,ic_datap += __tinc0_ic \
133             ,vc_datap += __tinc0_vc \
134             ,x_datap += __tinc0_x \
135             ,f_datap += __tinc0_f \
136             ,d_datap += __tinc0_d \
137             ,ierr_datap += __tinc0_ierr \
138             ,dx_datap += __tinc0_dx \
139             ,dy_dx_datap += __tinc0_dy_dx \
140             ) \
141             )
142             #define PDL_BROADCASTLOOP_END_pchip_chsp_readdata PDL_BROADCASTLOOP_END( \
143             __privtrans->broadcast, \
144             ic_datap -= __tinc1_ic * __tdims1 + __offsp[0]; \
145             vc_datap -= __tinc1_vc * __tdims1 + __offsp[1]; \
146             x_datap -= __tinc1_x * __tdims1 + __offsp[2]; \
147             f_datap -= __tinc1_f * __tdims1 + __offsp[3]; \
148             d_datap -= __tinc1_d * __tdims1 + __offsp[4]; \
149             ierr_datap -= __tinc1_ierr * __tdims1 + __offsp[5]; \
150             dx_datap -= __tinc1_dx * __tdims1 + __offsp[6]; \
151             dy_dx_datap -= __tinc1_dy_dx * __tdims1 + __offsp[7]; \
152             )
153 3           register PDL_Indx __inc_d_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,0)]; (void)__inc_d_n;
154 3           register PDL_Indx __inc_dx_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,6,0)]; (void)__inc_dx_n;
155 3           register PDL_Indx __inc_dy_dx_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,7,0)]; (void)__inc_dy_dx_n;
156 3           register PDL_Indx __inc_f_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_f_n;
157 3           register PDL_Indx __inc_ic_two = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_ic_two;
158 3           register PDL_Indx __inc_vc_two = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_vc_two;
159 3           register PDL_Indx __inc_x_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_x_n;
160             #ifndef PDL_DECLARE_PARAMS_pchip_chsp_1
161             #define PDL_DECLARE_PARAMS_pchip_chsp_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ic,PDL_PPSYM_PARAM_ic,PDL_TYPE_PARAM_ierr,PDL_PPSYM_PARAM_ierr) \
162             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ic, ic, (__privtrans->pdls[0]), 1, PDL_PPSYM_PARAM_ic) \
163             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, vc, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
164             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
165             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, f, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
166             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, d, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
167             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ierr, ierr, (__privtrans->pdls[5]), 1, PDL_PPSYM_PARAM_ierr) \
168             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, dx, (__privtrans->pdls[6]), 1, PDL_PPSYM_OP) \
169             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, dy_dx, (__privtrans->pdls[7]), 1, PDL_PPSYM_OP)
170             #endif
171             #define PDL_IF_BAD(t,f) f
172 3           switch (__privtrans->__datatype) { /* Start generic switch */
173 0           case PDL_F: {
174 0 0         PDL_DECLARE_PARAMS_pchip_chsp_1(PDL_Float,F,PDL_SByte,A,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
175 0 0         PDL_BROADCASTLOOP_START_pchip_chsp_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
176             #line 5206 "lib/PDL/Primitive.pd"
177             /* SINGULAR SYSTEM. */
178             /* *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** */
179             /* *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** */
180             #define dpchsp_singular(x, ind) \
181             (ierr_datap)[0] = -8; return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "SINGULAR LINEAR SYSTEM(" #x") at %td",ind);
182             /* Local variables */
183             PDL_Float stemp[3], xtemp[4];
184             PDL_Indx n = __privtrans->ind_sizes[0], nm1 = n - 1;
185             /* VALIDITY-CHECK ARGUMENTS. */
186             {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
187             #line 5216 "lib/PDL/Primitive.pd"
188             if ((x_datap)[0+(__inc_x_n*(n))] > (x_datap)[0+(__inc_x_n*(n-1))]) continue;
189             (ierr_datap)[0] = -1;
190             return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "X-ARRAY NOT STRICTLY INCREASING");
191             }} /* Close n=1 */
192             #line 5220 "lib/PDL/Primitive.pd"
193             PDL_Indx ibeg = (ic_datap)[0+(__inc_ic_two*(0))], iend = (ic_datap)[0+(__inc_ic_two*(1))], j;
194             (ierr_datap)[0] = 0;
195             if (PDL_ABS(ibeg) > 5)
196             --((ierr_datap)[0]);
197             if (PDL_ABS(iend) > 5)
198             (ierr_datap)[0] += -2;
199             if ((ierr_datap)[0] < 0) {
200             (ierr_datap)[0] += -3;
201             return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "IC OUT OF RANGE");
202             }
203             /* FUNCTION DEFINITION IS OK -- GO ON. */
204             /* COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, */
205             /* COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). */
206             {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
207             #line 5234 "lib/PDL/Primitive.pd"
208             (dx_datap)[0+(__inc_dx_n*(n))] = (x_datap)[0+(__inc_x_n*(n))] - (x_datap)[0+(__inc_x_n*(n-1))];
209             (dy_dx_datap)[0+(__inc_dy_dx_n*(n))] = ((f_datap)[0+(__inc_f_n*(n))] - (f_datap)[0+(__inc_f_n*(n-1))]) / (dx_datap)[0+(__inc_dx_n*(n))];
210             }} /* Close n=1 */
211             #line 5237 "lib/PDL/Primitive.pd"
212             /* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */
213             if (ibeg > n) {
214             ibeg = 0;
215             }
216             if (iend > n) {
217             iend = 0;
218             }
219             /* SET UP FOR BOUNDARY CONDITIONS. */
220             if (ibeg == 1 || ibeg == 2) {
221             (d_datap)[0+(__inc_d_n*(0))] = (vc_datap)[0+(__inc_vc_two*(0))];
222             } else if (ibeg > 2) {
223             /* PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. */
224             for (j = 0; j < ibeg; ++j) {
225             PDL_Indx index = ibeg - j + 1;
226             /* INDEX RUNS FROM IBEG DOWN TO 1. */
227             xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
228             if (j < ibeg-1)
229             stemp[j] = (dy_dx_datap)[0+(__inc_dy_dx_n*(index))];
230             }
231             /* -------------------------------- */
232            
233             /* PDL version: K, X, S are var names, 4th param output */
234             /* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */
235             /* DPCHDF: DPCHIP Finite Difference Formula */
236             /* Uses a divided difference formulation to compute a K-point approx- */
237             /* imation to the derivative at X(K) based on the data in X and S. */
238             /* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
239             /* derivative approximations. */
240             /* ---------------------------------------------------------------------- */
241             /* On input: */
242             /* K is the order of the desired derivative approximation. */
243             /* K must be at least 3 (error return if not). */
244             /* X contains the K values of the independent variable. */
245             /* X need not be ordered, but the values **MUST** be */
246             /* distinct. (Not checked here.) */
247             /* S contains the associated slope values: */
248             /* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
249             /* (Note that S need only be of length K-1.) */
250             /* On return: */
251             /* S will be destroyed. */
252             /* IERR will be set to -1 if K.LT.2 . */
253             /* DPCHDF will be set to the desired derivative approximation if */
254             /* IERR=0 or to zero if IERR=-1. */
255             /* ---------------------------------------------------------------------- */
256             /* ***SEE ALSO DPCHCE, DPCHSP */
257             /* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */
258             /* Verlag, New York, 1978, pp. 10-16. */
259             /* CHECK FOR LEGAL VALUE OF K. */
260             {
261             /* Local variables */
262             PDL_Indx i, j, k_cached = ibeg;
263             PDL_Float *x = xtemp, *s = stemp;
264             if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "K LESS THAN THREE");
265             /* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
266             for (j = 2; j < k_cached; ++j) {
267             PDL_Indx itmp = k_cached - j;
268             for (i = 0; i < itmp; ++i)
269             s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
270             }
271             /* EVALUATE DERIVATIVE AT X(K). */
272             PDL_Float value = s[0];
273             for (i = 1; i < k_cached-1; ++i)
274             value = s[i] + value * (x[k_cached-1] - x[i]);
275             (d_datap)[0+(__inc_d_n*(0))] = value;
276             }
277             ;
278             /* -------------------------------- */
279             ibeg = 1;
280             }
281             if (iend == 1 || iend == 2) {
282             (d_datap)[0+(__inc_d_n*(n-1))] = (vc_datap)[0+(__inc_vc_two*(1))];
283             } else if (iend > 2) {
284             /* PICK UP LAST IEND POINTS. */
285             for (j = 0; j < iend; ++j) {
286             PDL_Indx index = n - iend + j;
287             /* INDEX RUNS FROM N+1-IEND UP TO N. */
288             xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
289             if (j < iend-1)
290             stemp[j] = (dy_dx_datap)[0+(__inc_dy_dx_n*(index+1))];
291             }
292             /* -------------------------------- */
293            
294             /* PDL version: K, X, S are var names, 4th param output */
295             /* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */
296             /* DPCHDF: DPCHIP Finite Difference Formula */
297             /* Uses a divided difference formulation to compute a K-point approx- */
298             /* imation to the derivative at X(K) based on the data in X and S. */
299             /* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
300             /* derivative approximations. */
301             /* ---------------------------------------------------------------------- */
302             /* On input: */
303             /* K is the order of the desired derivative approximation. */
304             /* K must be at least 3 (error return if not). */
305             /* X contains the K values of the independent variable. */
306             /* X need not be ordered, but the values **MUST** be */
307             /* distinct. (Not checked here.) */
308             /* S contains the associated slope values: */
309             /* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
310             /* (Note that S need only be of length K-1.) */
311             /* On return: */
312             /* S will be destroyed. */
313             /* IERR will be set to -1 if K.LT.2 . */
314             /* DPCHDF will be set to the desired derivative approximation if */
315             /* IERR=0 or to zero if IERR=-1. */
316             /* ---------------------------------------------------------------------- */
317             /* ***SEE ALSO DPCHCE, DPCHSP */
318             /* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */
319             /* Verlag, New York, 1978, pp. 10-16. */
320             /* CHECK FOR LEGAL VALUE OF K. */
321             {
322             /* Local variables */
323             PDL_Indx i, j, k_cached = iend;
324             PDL_Float *x = xtemp, *s = stemp;
325             if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "K LESS THAN THREE");
326             /* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
327             for (j = 2; j < k_cached; ++j) {
328             PDL_Indx itmp = k_cached - j;
329             for (i = 0; i < itmp; ++i)
330             s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
331             }
332             /* EVALUATE DERIVATIVE AT X(K). */
333             PDL_Float value = s[0];
334             for (i = 1; i < k_cached-1; ++i)
335             value = s[i] + value * (x[k_cached-1] - x[i]);
336             (d_datap)[0+(__inc_d_n*(n-1))] = value;
337             }
338             ;
339             /* -------------------------------- */
340             iend = 1;
341             }
342             /* --------------------( BEGIN CODING FROM CUBSPL )-------------------- */
343             /* **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF */
344             /* F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- */
345             /* INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. */
346             /* WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. */
347             /* CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM */
348             /* WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) */
349             if (ibeg == 0) {
350             if (n == 2) {
351             /* NO CONDITION AT LEFT END AND N = 2. */
352             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 1.;
353             (dx_datap)[0+(__inc_dx_n*(0))] = 1.;
354             (d_datap)[0+(__inc_d_n*(0))] = 2. * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))];
355             } else {
356             /* NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. */
357             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = (dx_datap)[0+(__inc_dx_n*(2))];
358             (dx_datap)[0+(__inc_dx_n*(0))] = (dx_datap)[0+(__inc_dx_n*(1))] + (dx_datap)[0+(__inc_dx_n*(2))];
359             /* Computing 2nd power */
360             (d_datap)[0+(__inc_d_n*(0))] = (((dx_datap)[0+(__inc_dx_n*(1))] + 2. * (dx_datap)[0+(__inc_dx_n*(0))]) * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))] * (dx_datap)[0+(__inc_dx_n*(2))] + (dx_datap)[0+(__inc_dx_n*(1))] *
361             (dx_datap)[0+(__inc_dx_n*(1))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(2))]) / (dx_datap)[0+(__inc_dx_n*(0))];
362             }
363             } else if (ibeg == 1) {
364             /* SLOPE PRESCRIBED AT LEFT END. */
365             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 1.;
366             (dx_datap)[0+(__inc_dx_n*(0))] = 0.;
367             } else {
368             /* SECOND DERIVATIVE PRESCRIBED AT LEFT END. */
369             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 2.;
370             (dx_datap)[0+(__inc_dx_n*(0))] = 1.;
371             (d_datap)[0+(__inc_d_n*(0))] = 3. * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))] - 0.5 * (dx_datap)[0+(__inc_dx_n*(1))] * (d_datap)[0+(__inc_d_n*(0))];
372             }
373             /* IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND */
374             /* CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH */
375             /* EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). */
376             if (n > 2) {
377             {/* Open n=1:-1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size-1)); for(; n<__n_stop; n+=1) {
378             #line 5313 "lib/PDL/Primitive.pd"
379             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] == 0.) {
380             dpchsp_singular(1, n-1);
381             }
382             PDL_Float g = -(dx_datap)[0+(__inc_dx_n*(n+1))] / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
383             (d_datap)[0+(__inc_d_n*(n))] = g * (d_datap)[0+(__inc_d_n*(n-1))] + 3. * ((dx_datap)[0+(__inc_dx_n*(n))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(n+1))] +
384             (dx_datap)[0+(__inc_dx_n*(n+1))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(n))]);
385             (dy_dx_datap)[0+(__inc_dy_dx_n*(n))] = g * (dx_datap)[0+(__inc_dx_n*(n-1))] + 2. * ((dx_datap)[0+(__inc_dx_n*(n))] + (dx_datap)[0+(__inc_dx_n*(n+1))]);
386             }} /* Close n=1:-1 */
387             #line 5321 "lib/PDL/Primitive.pd"
388             }
389             /* CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM */
390             /* (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) */
391             /* IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- */
392             /* SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT */
393             /* AT THIS POINT. */
394             if (iend != 1) {
395             if (iend == 0 && n == 2 && ibeg == 0) {
396             /* NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. */
397             (d_datap)[0+(__inc_d_n*(1))] = (dy_dx_datap)[0+(__inc_dy_dx_n*(1))];
398             } else {
399             PDL_Float g;
400             if (iend == 0) {
401             if (n == 2 || (n == 3 && ibeg == 0)) {
402             /* EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* */
403             /* NOT-A-KNOT AT LEFT END POINT). */
404             (d_datap)[0+(__inc_d_n*(n-1))] = 2. * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
405             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = 1.;
406             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
407             dpchsp_singular(2, n-2);
408             }
409             g = -1. / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
410             } else {
411             /* NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- */
412             /* KNOT AT LEFT END POINT. */
413             g = (dx_datap)[0+(__inc_dx_n*(n-2))] + (dx_datap)[0+(__inc_dx_n*(n-1))];
414             /* DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). */
415             /* Computing 2nd power */
416             PDL_Float dtmp = (dx_datap)[0+(__inc_dx_n*(n-1))];
417             (d_datap)[0+(__inc_d_n*(n-1))] = (((dx_datap)[0+(__inc_dx_n*(n-1))] + 2. * g) * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] * (dx_datap)[0+(__inc_dx_n*(n-2))] + dtmp * dtmp * ((f_datap)[0+(__inc_f_n*(n-2))] - (f_datap)[0+(__inc_f_n*(n-3))]) / (dx_datap)[0+(__inc_dx_n*(n-2))]) / g;
418             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
419             dpchsp_singular(3, n-2);
420             }
421             g /= -(dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
422             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = (dx_datap)[0+(__inc_dx_n*(n-2))];
423             }
424             } else {
425             /* SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. */
426             (d_datap)[0+(__inc_d_n*(n-1))] = 3. * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] + 0.5 * (dx_datap)[0+(__inc_dx_n*(n-1))] * (d_datap)[0+(__inc_d_n*(n-1))];
427             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = 2.;
428             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
429             dpchsp_singular(4, n-2);
430             }
431             g = -1. / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
432             }
433             /* COMPLETE FORWARD PASS OF GAUSS ELIMINATION. */
434             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = g * (dx_datap)[0+(__inc_dx_n*(n-2))] + (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
435             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] == 0.) {
436             dpchsp_singular(5, n-1);
437             }
438             (d_datap)[0+(__inc_d_n*(n-1))] = (g * (d_datap)[0+(__inc_d_n*(n-2))] + (d_datap)[0+(__inc_d_n*(n-1))]) / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
439             }
440             }
441             /* CARRY OUT BACK SUBSTITUTION */
442             {/* Open n=nm1-1::-1 */ PDL_EXPAND2(register PDL_Indx n=PDLMIN(nm1-1, (__n_size-1)), __n_stop=0); for(; n>=__n_stop; n+=-1) {
443             #line 5376 "lib/PDL/Primitive.pd"
444             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n))] == 0.) {
445             dpchsp_singular(6, n);
446             }
447             (d_datap)[0+(__inc_d_n*(n))] = ((d_datap)[0+(__inc_d_n*(n))] - (dx_datap)[0+(__inc_dx_n*(n))] * (d_datap)[0+(__inc_d_n*(n+1))]) / (dy_dx_datap)[0+(__inc_dy_dx_n*(n))];
448             }} /* Close n=nm1-1::-1 */
449             #line 5381 "lib/PDL/Primitive.pd"
450             /* --------------------( END CODING FROM CUBSPL )-------------------- */
451             #undef dpchsp_singular
452             #line 453 "lib/PDL/Primitive-pp-pchip_chsp.c"
453 0 0         }PDL_BROADCASTLOOP_END_pchip_chsp_readdata
    0          
454 0           } break;
455 3           case PDL_D: {
456 3 50         PDL_DECLARE_PARAMS_pchip_chsp_1(PDL_Double,D,PDL_SByte,A,PDL_Indx,N)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
457 14 50         PDL_BROADCASTLOOP_START_pchip_chsp_readdata {
    50          
    50          
    50          
    50          
    100          
    100          
458             #line 5206 "lib/PDL/Primitive.pd"
459             /* SINGULAR SYSTEM. */
460             /* *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** */
461             /* *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** */
462             #define dpchsp_singular(x, ind) \
463             (ierr_datap)[0] = -8; return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "SINGULAR LINEAR SYSTEM(" #x") at %td",ind);
464             /* Local variables */
465             PDL_Double stemp[3], xtemp[4];
466             PDL_Indx n = __privtrans->ind_sizes[0], nm1 = n - 1;
467             /* VALIDITY-CHECK ARGUMENTS. */
468             {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
469             #line 5216 "lib/PDL/Primitive.pd"
470             if ((x_datap)[0+(__inc_x_n*(n))] > (x_datap)[0+(__inc_x_n*(n-1))]) continue;
471             (ierr_datap)[0] = -1;
472             return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "X-ARRAY NOT STRICTLY INCREASING");
473             }} /* Close n=1 */
474             #line 5220 "lib/PDL/Primitive.pd"
475             PDL_Indx ibeg = (ic_datap)[0+(__inc_ic_two*(0))], iend = (ic_datap)[0+(__inc_ic_two*(1))], j;
476             (ierr_datap)[0] = 0;
477             if (PDL_ABS(ibeg) > 5)
478             --((ierr_datap)[0]);
479             if (PDL_ABS(iend) > 5)
480             (ierr_datap)[0] += -2;
481             if ((ierr_datap)[0] < 0) {
482             (ierr_datap)[0] += -3;
483             return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "IC OUT OF RANGE");
484             }
485             /* FUNCTION DEFINITION IS OK -- GO ON. */
486             /* COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, */
487             /* COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). */
488             {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
489             #line 5234 "lib/PDL/Primitive.pd"
490             (dx_datap)[0+(__inc_dx_n*(n))] = (x_datap)[0+(__inc_x_n*(n))] - (x_datap)[0+(__inc_x_n*(n-1))];
491             (dy_dx_datap)[0+(__inc_dy_dx_n*(n))] = ((f_datap)[0+(__inc_f_n*(n))] - (f_datap)[0+(__inc_f_n*(n-1))]) / (dx_datap)[0+(__inc_dx_n*(n))];
492             }} /* Close n=1 */
493             #line 5237 "lib/PDL/Primitive.pd"
494             /* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */
495             if (ibeg > n) {
496             ibeg = 0;
497             }
498             if (iend > n) {
499             iend = 0;
500             }
501             /* SET UP FOR BOUNDARY CONDITIONS. */
502             if (ibeg == 1 || ibeg == 2) {
503             (d_datap)[0+(__inc_d_n*(0))] = (vc_datap)[0+(__inc_vc_two*(0))];
504             } else if (ibeg > 2) {
505             /* PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. */
506             for (j = 0; j < ibeg; ++j) {
507             PDL_Indx index = ibeg - j + 1;
508             /* INDEX RUNS FROM IBEG DOWN TO 1. */
509             xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
510             if (j < ibeg-1)
511             stemp[j] = (dy_dx_datap)[0+(__inc_dy_dx_n*(index))];
512             }
513             /* -------------------------------- */
514            
515             /* PDL version: K, X, S are var names, 4th param output */
516             /* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */
517             /* DPCHDF: DPCHIP Finite Difference Formula */
518             /* Uses a divided difference formulation to compute a K-point approx- */
519             /* imation to the derivative at X(K) based on the data in X and S. */
520             /* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
521             /* derivative approximations. */
522             /* ---------------------------------------------------------------------- */
523             /* On input: */
524             /* K is the order of the desired derivative approximation. */
525             /* K must be at least 3 (error return if not). */
526             /* X contains the K values of the independent variable. */
527             /* X need not be ordered, but the values **MUST** be */
528             /* distinct. (Not checked here.) */
529             /* S contains the associated slope values: */
530             /* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
531             /* (Note that S need only be of length K-1.) */
532             /* On return: */
533             /* S will be destroyed. */
534             /* IERR will be set to -1 if K.LT.2 . */
535             /* DPCHDF will be set to the desired derivative approximation if */
536             /* IERR=0 or to zero if IERR=-1. */
537             /* ---------------------------------------------------------------------- */
538             /* ***SEE ALSO DPCHCE, DPCHSP */
539             /* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */
540             /* Verlag, New York, 1978, pp. 10-16. */
541             /* CHECK FOR LEGAL VALUE OF K. */
542             {
543             /* Local variables */
544             PDL_Indx i, j, k_cached = ibeg;
545             PDL_Double *x = xtemp, *s = stemp;
546             if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "K LESS THAN THREE");
547             /* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
548             for (j = 2; j < k_cached; ++j) {
549             PDL_Indx itmp = k_cached - j;
550             for (i = 0; i < itmp; ++i)
551             s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
552             }
553             /* EVALUATE DERIVATIVE AT X(K). */
554             PDL_Double value = s[0];
555             for (i = 1; i < k_cached-1; ++i)
556             value = s[i] + value * (x[k_cached-1] - x[i]);
557             (d_datap)[0+(__inc_d_n*(0))] = value;
558             }
559             ;
560             /* -------------------------------- */
561             ibeg = 1;
562             }
563             if (iend == 1 || iend == 2) {
564             (d_datap)[0+(__inc_d_n*(n-1))] = (vc_datap)[0+(__inc_vc_two*(1))];
565             } else if (iend > 2) {
566             /* PICK UP LAST IEND POINTS. */
567             for (j = 0; j < iend; ++j) {
568             PDL_Indx index = n - iend + j;
569             /* INDEX RUNS FROM N+1-IEND UP TO N. */
570             xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
571             if (j < iend-1)
572             stemp[j] = (dy_dx_datap)[0+(__inc_dy_dx_n*(index+1))];
573             }
574             /* -------------------------------- */
575            
576             /* PDL version: K, X, S are var names, 4th param output */
577             /* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */
578             /* DPCHDF: DPCHIP Finite Difference Formula */
579             /* Uses a divided difference formulation to compute a K-point approx- */
580             /* imation to the derivative at X(K) based on the data in X and S. */
581             /* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
582             /* derivative approximations. */
583             /* ---------------------------------------------------------------------- */
584             /* On input: */
585             /* K is the order of the desired derivative approximation. */
586             /* K must be at least 3 (error return if not). */
587             /* X contains the K values of the independent variable. */
588             /* X need not be ordered, but the values **MUST** be */
589             /* distinct. (Not checked here.) */
590             /* S contains the associated slope values: */
591             /* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
592             /* (Note that S need only be of length K-1.) */
593             /* On return: */
594             /* S will be destroyed. */
595             /* IERR will be set to -1 if K.LT.2 . */
596             /* DPCHDF will be set to the desired derivative approximation if */
597             /* IERR=0 or to zero if IERR=-1. */
598             /* ---------------------------------------------------------------------- */
599             /* ***SEE ALSO DPCHCE, DPCHSP */
600             /* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */
601             /* Verlag, New York, 1978, pp. 10-16. */
602             /* CHECK FOR LEGAL VALUE OF K. */
603             {
604             /* Local variables */
605             PDL_Indx i, j, k_cached = iend;
606             PDL_Double *x = xtemp, *s = stemp;
607             if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "K LESS THAN THREE");
608             /* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
609             for (j = 2; j < k_cached; ++j) {
610             PDL_Indx itmp = k_cached - j;
611             for (i = 0; i < itmp; ++i)
612             s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
613             }
614             /* EVALUATE DERIVATIVE AT X(K). */
615             PDL_Double value = s[0];
616             for (i = 1; i < k_cached-1; ++i)
617             value = s[i] + value * (x[k_cached-1] - x[i]);
618             (d_datap)[0+(__inc_d_n*(n-1))] = value;
619             }
620             ;
621             /* -------------------------------- */
622             iend = 1;
623             }
624             /* --------------------( BEGIN CODING FROM CUBSPL )-------------------- */
625             /* **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF */
626             /* F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- */
627             /* INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. */
628             /* WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. */
629             /* CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM */
630             /* WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) */
631             if (ibeg == 0) {
632             if (n == 2) {
633             /* NO CONDITION AT LEFT END AND N = 2. */
634             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 1.;
635             (dx_datap)[0+(__inc_dx_n*(0))] = 1.;
636             (d_datap)[0+(__inc_d_n*(0))] = 2. * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))];
637             } else {
638             /* NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. */
639             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = (dx_datap)[0+(__inc_dx_n*(2))];
640             (dx_datap)[0+(__inc_dx_n*(0))] = (dx_datap)[0+(__inc_dx_n*(1))] + (dx_datap)[0+(__inc_dx_n*(2))];
641             /* Computing 2nd power */
642             (d_datap)[0+(__inc_d_n*(0))] = (((dx_datap)[0+(__inc_dx_n*(1))] + 2. * (dx_datap)[0+(__inc_dx_n*(0))]) * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))] * (dx_datap)[0+(__inc_dx_n*(2))] + (dx_datap)[0+(__inc_dx_n*(1))] *
643             (dx_datap)[0+(__inc_dx_n*(1))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(2))]) / (dx_datap)[0+(__inc_dx_n*(0))];
644             }
645             } else if (ibeg == 1) {
646             /* SLOPE PRESCRIBED AT LEFT END. */
647             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 1.;
648             (dx_datap)[0+(__inc_dx_n*(0))] = 0.;
649             } else {
650             /* SECOND DERIVATIVE PRESCRIBED AT LEFT END. */
651             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 2.;
652             (dx_datap)[0+(__inc_dx_n*(0))] = 1.;
653             (d_datap)[0+(__inc_d_n*(0))] = 3. * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))] - 0.5 * (dx_datap)[0+(__inc_dx_n*(1))] * (d_datap)[0+(__inc_d_n*(0))];
654             }
655             /* IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND */
656             /* CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH */
657             /* EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). */
658             if (n > 2) {
659             {/* Open n=1:-1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size-1)); for(; n<__n_stop; n+=1) {
660             #line 5313 "lib/PDL/Primitive.pd"
661             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] == 0.) {
662             dpchsp_singular(1, n-1);
663             }
664             PDL_Double g = -(dx_datap)[0+(__inc_dx_n*(n+1))] / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
665             (d_datap)[0+(__inc_d_n*(n))] = g * (d_datap)[0+(__inc_d_n*(n-1))] + 3. * ((dx_datap)[0+(__inc_dx_n*(n))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(n+1))] +
666             (dx_datap)[0+(__inc_dx_n*(n+1))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(n))]);
667             (dy_dx_datap)[0+(__inc_dy_dx_n*(n))] = g * (dx_datap)[0+(__inc_dx_n*(n-1))] + 2. * ((dx_datap)[0+(__inc_dx_n*(n))] + (dx_datap)[0+(__inc_dx_n*(n+1))]);
668             }} /* Close n=1:-1 */
669             #line 5321 "lib/PDL/Primitive.pd"
670             }
671             /* CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM */
672             /* (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) */
673             /* IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- */
674             /* SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT */
675             /* AT THIS POINT. */
676             if (iend != 1) {
677             if (iend == 0 && n == 2 && ibeg == 0) {
678             /* NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. */
679             (d_datap)[0+(__inc_d_n*(1))] = (dy_dx_datap)[0+(__inc_dy_dx_n*(1))];
680             } else {
681             PDL_Double g;
682             if (iend == 0) {
683             if (n == 2 || (n == 3 && ibeg == 0)) {
684             /* EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* */
685             /* NOT-A-KNOT AT LEFT END POINT). */
686             (d_datap)[0+(__inc_d_n*(n-1))] = 2. * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
687             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = 1.;
688             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
689             dpchsp_singular(2, n-2);
690             }
691             g = -1. / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
692             } else {
693             /* NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- */
694             /* KNOT AT LEFT END POINT. */
695             g = (dx_datap)[0+(__inc_dx_n*(n-2))] + (dx_datap)[0+(__inc_dx_n*(n-1))];
696             /* DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). */
697             /* Computing 2nd power */
698             PDL_Double dtmp = (dx_datap)[0+(__inc_dx_n*(n-1))];
699             (d_datap)[0+(__inc_d_n*(n-1))] = (((dx_datap)[0+(__inc_dx_n*(n-1))] + 2. * g) * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] * (dx_datap)[0+(__inc_dx_n*(n-2))] + dtmp * dtmp * ((f_datap)[0+(__inc_f_n*(n-2))] - (f_datap)[0+(__inc_f_n*(n-3))]) / (dx_datap)[0+(__inc_dx_n*(n-2))]) / g;
700             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
701             dpchsp_singular(3, n-2);
702             }
703             g /= -(dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
704             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = (dx_datap)[0+(__inc_dx_n*(n-2))];
705             }
706             } else {
707             /* SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. */
708             (d_datap)[0+(__inc_d_n*(n-1))] = 3. * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] + 0.5 * (dx_datap)[0+(__inc_dx_n*(n-1))] * (d_datap)[0+(__inc_d_n*(n-1))];
709             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = 2.;
710             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
711             dpchsp_singular(4, n-2);
712             }
713             g = -1. / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
714             }
715             /* COMPLETE FORWARD PASS OF GAUSS ELIMINATION. */
716             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = g * (dx_datap)[0+(__inc_dx_n*(n-2))] + (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
717             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] == 0.) {
718             dpchsp_singular(5, n-1);
719             }
720             (d_datap)[0+(__inc_d_n*(n-1))] = (g * (d_datap)[0+(__inc_d_n*(n-2))] + (d_datap)[0+(__inc_d_n*(n-1))]) / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
721             }
722             }
723             /* CARRY OUT BACK SUBSTITUTION */
724             {/* Open n=nm1-1::-1 */ PDL_EXPAND2(register PDL_Indx n=PDLMIN(nm1-1, (__n_size-1)), __n_stop=0); for(; n>=__n_stop; n+=-1) {
725             #line 5376 "lib/PDL/Primitive.pd"
726             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n))] == 0.) {
727             dpchsp_singular(6, n);
728             }
729             (d_datap)[0+(__inc_d_n*(n))] = ((d_datap)[0+(__inc_d_n*(n))] - (dx_datap)[0+(__inc_dx_n*(n))] * (d_datap)[0+(__inc_d_n*(n+1))]) / (dy_dx_datap)[0+(__inc_dy_dx_n*(n))];
730             }} /* Close n=nm1-1::-1 */
731             #line 5381 "lib/PDL/Primitive.pd"
732             /* --------------------( END CODING FROM CUBSPL )-------------------- */
733             #undef dpchsp_singular
734             #line 735 "lib/PDL/Primitive-pp-pchip_chsp.c"
735 3 50         }PDL_BROADCASTLOOP_END_pchip_chsp_readdata
    50          
736 3           } break;
737 0           case PDL_LD: {
738 0 0         PDL_DECLARE_PARAMS_pchip_chsp_1(PDL_LDouble,E,PDL_SByte,A,PDL_Indx,N)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
739 0 0         PDL_BROADCASTLOOP_START_pchip_chsp_readdata {
    0          
    0          
    0          
    0          
    0          
    0          
740             #line 5206 "lib/PDL/Primitive.pd"
741             /* SINGULAR SYSTEM. */
742             /* *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES *** */
743             /* *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). *** */
744             #define dpchsp_singular(x, ind) \
745             (ierr_datap)[0] = -8; return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "SINGULAR LINEAR SYSTEM(" #x") at %td",ind);
746             /* Local variables */
747             PDL_LDouble stemp[3], xtemp[4];
748             PDL_Indx n = __privtrans->ind_sizes[0], nm1 = n - 1;
749             /* VALIDITY-CHECK ARGUMENTS. */
750             {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
751             #line 5216 "lib/PDL/Primitive.pd"
752             if ((x_datap)[0+(__inc_x_n*(n))] > (x_datap)[0+(__inc_x_n*(n-1))]) continue;
753             (ierr_datap)[0] = -1;
754             return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "X-ARRAY NOT STRICTLY INCREASING");
755             }} /* Close n=1 */
756             #line 5220 "lib/PDL/Primitive.pd"
757             PDL_Indx ibeg = (ic_datap)[0+(__inc_ic_two*(0))], iend = (ic_datap)[0+(__inc_ic_two*(1))], j;
758             (ierr_datap)[0] = 0;
759             if (PDL_ABS(ibeg) > 5)
760             --((ierr_datap)[0]);
761             if (PDL_ABS(iend) > 5)
762             (ierr_datap)[0] += -2;
763             if ((ierr_datap)[0] < 0) {
764             (ierr_datap)[0] += -3;
765             return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "IC OUT OF RANGE");
766             }
767             /* FUNCTION DEFINITION IS OK -- GO ON. */
768             /* COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, */
769             /* COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). */
770             {/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
771             #line 5234 "lib/PDL/Primitive.pd"
772             (dx_datap)[0+(__inc_dx_n*(n))] = (x_datap)[0+(__inc_x_n*(n))] - (x_datap)[0+(__inc_x_n*(n-1))];
773             (dy_dx_datap)[0+(__inc_dy_dx_n*(n))] = ((f_datap)[0+(__inc_f_n*(n))] - (f_datap)[0+(__inc_f_n*(n-1))]) / (dx_datap)[0+(__inc_dx_n*(n))];
774             }} /* Close n=1 */
775             #line 5237 "lib/PDL/Primitive.pd"
776             /* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */
777             if (ibeg > n) {
778             ibeg = 0;
779             }
780             if (iend > n) {
781             iend = 0;
782             }
783             /* SET UP FOR BOUNDARY CONDITIONS. */
784             if (ibeg == 1 || ibeg == 2) {
785             (d_datap)[0+(__inc_d_n*(0))] = (vc_datap)[0+(__inc_vc_two*(0))];
786             } else if (ibeg > 2) {
787             /* PICK UP FIRST IBEG POINTS, IN REVERSE ORDER. */
788             for (j = 0; j < ibeg; ++j) {
789             PDL_Indx index = ibeg - j + 1;
790             /* INDEX RUNS FROM IBEG DOWN TO 1. */
791             xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
792             if (j < ibeg-1)
793             stemp[j] = (dy_dx_datap)[0+(__inc_dy_dx_n*(index))];
794             }
795             /* -------------------------------- */
796            
797             /* PDL version: K, X, S are var names, 4th param output */
798             /* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */
799             /* DPCHDF: DPCHIP Finite Difference Formula */
800             /* Uses a divided difference formulation to compute a K-point approx- */
801             /* imation to the derivative at X(K) based on the data in X and S. */
802             /* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
803             /* derivative approximations. */
804             /* ---------------------------------------------------------------------- */
805             /* On input: */
806             /* K is the order of the desired derivative approximation. */
807             /* K must be at least 3 (error return if not). */
808             /* X contains the K values of the independent variable. */
809             /* X need not be ordered, but the values **MUST** be */
810             /* distinct. (Not checked here.) */
811             /* S contains the associated slope values: */
812             /* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
813             /* (Note that S need only be of length K-1.) */
814             /* On return: */
815             /* S will be destroyed. */
816             /* IERR will be set to -1 if K.LT.2 . */
817             /* DPCHDF will be set to the desired derivative approximation if */
818             /* IERR=0 or to zero if IERR=-1. */
819             /* ---------------------------------------------------------------------- */
820             /* ***SEE ALSO DPCHCE, DPCHSP */
821             /* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */
822             /* Verlag, New York, 1978, pp. 10-16. */
823             /* CHECK FOR LEGAL VALUE OF K. */
824             {
825             /* Local variables */
826             PDL_Indx i, j, k_cached = ibeg;
827             PDL_LDouble *x = xtemp, *s = stemp;
828             if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "K LESS THAN THREE");
829             /* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
830             for (j = 2; j < k_cached; ++j) {
831             PDL_Indx itmp = k_cached - j;
832             for (i = 0; i < itmp; ++i)
833             s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
834             }
835             /* EVALUATE DERIVATIVE AT X(K). */
836             PDL_LDouble value = s[0];
837             for (i = 1; i < k_cached-1; ++i)
838             value = s[i] + value * (x[k_cached-1] - x[i]);
839             (d_datap)[0+(__inc_d_n*(0))] = value;
840             }
841             ;
842             /* -------------------------------- */
843             ibeg = 1;
844             }
845             if (iend == 1 || iend == 2) {
846             (d_datap)[0+(__inc_d_n*(n-1))] = (vc_datap)[0+(__inc_vc_two*(1))];
847             } else if (iend > 2) {
848             /* PICK UP LAST IEND POINTS. */
849             for (j = 0; j < iend; ++j) {
850             PDL_Indx index = n - iend + j;
851             /* INDEX RUNS FROM N+1-IEND UP TO N. */
852             xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
853             if (j < iend-1)
854             stemp[j] = (dy_dx_datap)[0+(__inc_dy_dx_n*(index+1))];
855             }
856             /* -------------------------------- */
857            
858             /* PDL version: K, X, S are var names, 4th param output */
859             /* ***PURPOSE Computes divided differences for DPCHCE and DPCHSP */
860             /* DPCHDF: DPCHIP Finite Difference Formula */
861             /* Uses a divided difference formulation to compute a K-point approx- */
862             /* imation to the derivative at X(K) based on the data in X and S. */
863             /* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
864             /* derivative approximations. */
865             /* ---------------------------------------------------------------------- */
866             /* On input: */
867             /* K is the order of the desired derivative approximation. */
868             /* K must be at least 3 (error return if not). */
869             /* X contains the K values of the independent variable. */
870             /* X need not be ordered, but the values **MUST** be */
871             /* distinct. (Not checked here.) */
872             /* S contains the associated slope values: */
873             /* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
874             /* (Note that S need only be of length K-1.) */
875             /* On return: */
876             /* S will be destroyed. */
877             /* IERR will be set to -1 if K.LT.2 . */
878             /* DPCHDF will be set to the desired derivative approximation if */
879             /* IERR=0 or to zero if IERR=-1. */
880             /* ---------------------------------------------------------------------- */
881             /* ***SEE ALSO DPCHCE, DPCHSP */
882             /* ***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- */
883             /* Verlag, New York, 1978, pp. 10-16. */
884             /* CHECK FOR LEGAL VALUE OF K. */
885             {
886             /* Local variables */
887             PDL_Indx i, j, k_cached = iend;
888             PDL_LDouble *x = xtemp, *s = stemp;
889             if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chsp:" "K LESS THAN THREE");
890             /* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
891             for (j = 2; j < k_cached; ++j) {
892             PDL_Indx itmp = k_cached - j;
893             for (i = 0; i < itmp; ++i)
894             s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
895             }
896             /* EVALUATE DERIVATIVE AT X(K). */
897             PDL_LDouble value = s[0];
898             for (i = 1; i < k_cached-1; ++i)
899             value = s[i] + value * (x[k_cached-1] - x[i]);
900             (d_datap)[0+(__inc_d_n*(n-1))] = value;
901             }
902             ;
903             /* -------------------------------- */
904             iend = 1;
905             }
906             /* --------------------( BEGIN CODING FROM CUBSPL )-------------------- */
907             /* **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF */
908             /* F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- */
909             /* INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. */
910             /* WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. */
911             /* CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM */
912             /* WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) */
913             if (ibeg == 0) {
914             if (n == 2) {
915             /* NO CONDITION AT LEFT END AND N = 2. */
916             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 1.;
917             (dx_datap)[0+(__inc_dx_n*(0))] = 1.;
918             (d_datap)[0+(__inc_d_n*(0))] = 2. * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))];
919             } else {
920             /* NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2. */
921             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = (dx_datap)[0+(__inc_dx_n*(2))];
922             (dx_datap)[0+(__inc_dx_n*(0))] = (dx_datap)[0+(__inc_dx_n*(1))] + (dx_datap)[0+(__inc_dx_n*(2))];
923             /* Computing 2nd power */
924             (d_datap)[0+(__inc_d_n*(0))] = (((dx_datap)[0+(__inc_dx_n*(1))] + 2. * (dx_datap)[0+(__inc_dx_n*(0))]) * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))] * (dx_datap)[0+(__inc_dx_n*(2))] + (dx_datap)[0+(__inc_dx_n*(1))] *
925             (dx_datap)[0+(__inc_dx_n*(1))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(2))]) / (dx_datap)[0+(__inc_dx_n*(0))];
926             }
927             } else if (ibeg == 1) {
928             /* SLOPE PRESCRIBED AT LEFT END. */
929             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 1.;
930             (dx_datap)[0+(__inc_dx_n*(0))] = 0.;
931             } else {
932             /* SECOND DERIVATIVE PRESCRIBED AT LEFT END. */
933             (dy_dx_datap)[0+(__inc_dy_dx_n*(0))] = 2.;
934             (dx_datap)[0+(__inc_dx_n*(0))] = 1.;
935             (d_datap)[0+(__inc_d_n*(0))] = 3. * (dy_dx_datap)[0+(__inc_dy_dx_n*(1))] - 0.5 * (dx_datap)[0+(__inc_dx_n*(1))] * (d_datap)[0+(__inc_d_n*(0))];
936             }
937             /* IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND */
938             /* CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH */
939             /* EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). */
940             if (n > 2) {
941             {/* Open n=1:-1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size-1)); for(; n<__n_stop; n+=1) {
942             #line 5313 "lib/PDL/Primitive.pd"
943             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] == 0.) {
944             dpchsp_singular(1, n-1);
945             }
946             PDL_LDouble g = -(dx_datap)[0+(__inc_dx_n*(n+1))] / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
947             (d_datap)[0+(__inc_d_n*(n))] = g * (d_datap)[0+(__inc_d_n*(n-1))] + 3. * ((dx_datap)[0+(__inc_dx_n*(n))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(n+1))] +
948             (dx_datap)[0+(__inc_dx_n*(n+1))] * (dy_dx_datap)[0+(__inc_dy_dx_n*(n))]);
949             (dy_dx_datap)[0+(__inc_dy_dx_n*(n))] = g * (dx_datap)[0+(__inc_dx_n*(n-1))] + 2. * ((dx_datap)[0+(__inc_dx_n*(n))] + (dx_datap)[0+(__inc_dx_n*(n+1))]);
950             }} /* Close n=1:-1 */
951             #line 5321 "lib/PDL/Primitive.pd"
952             }
953             /* CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM */
954             /* (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) */
955             /* IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- */
956             /* SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT */
957             /* AT THIS POINT. */
958             if (iend != 1) {
959             if (iend == 0 && n == 2 && ibeg == 0) {
960             /* NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2. */
961             (d_datap)[0+(__inc_d_n*(1))] = (dy_dx_datap)[0+(__inc_dy_dx_n*(1))];
962             } else {
963             PDL_LDouble g;
964             if (iend == 0) {
965             if (n == 2 || (n == 3 && ibeg == 0)) {
966             /* EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT* */
967             /* NOT-A-KNOT AT LEFT END POINT). */
968             (d_datap)[0+(__inc_d_n*(n-1))] = 2. * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
969             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = 1.;
970             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
971             dpchsp_singular(2, n-2);
972             }
973             g = -1. / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
974             } else {
975             /* NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR ALSO NOT-A- */
976             /* KNOT AT LEFT END POINT. */
977             g = (dx_datap)[0+(__inc_dx_n*(n-2))] + (dx_datap)[0+(__inc_dx_n*(n-1))];
978             /* DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES). */
979             /* Computing 2nd power */
980             PDL_LDouble dtmp = (dx_datap)[0+(__inc_dx_n*(n-1))];
981             (d_datap)[0+(__inc_d_n*(n-1))] = (((dx_datap)[0+(__inc_dx_n*(n-1))] + 2. * g) * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] * (dx_datap)[0+(__inc_dx_n*(n-2))] + dtmp * dtmp * ((f_datap)[0+(__inc_f_n*(n-2))] - (f_datap)[0+(__inc_f_n*(n-3))]) / (dx_datap)[0+(__inc_dx_n*(n-2))]) / g;
982             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
983             dpchsp_singular(3, n-2);
984             }
985             g /= -(dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
986             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = (dx_datap)[0+(__inc_dx_n*(n-2))];
987             }
988             } else {
989             /* SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT. */
990             (d_datap)[0+(__inc_d_n*(n-1))] = 3. * (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] + 0.5 * (dx_datap)[0+(__inc_dx_n*(n-1))] * (d_datap)[0+(__inc_d_n*(n-1))];
991             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = 2.;
992             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))] == 0.) {
993             dpchsp_singular(4, n-2);
994             }
995             g = -1. / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-2))];
996             }
997             /* COMPLETE FORWARD PASS OF GAUSS ELIMINATION. */
998             (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] = g * (dx_datap)[0+(__inc_dx_n*(n-2))] + (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
999             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))] == 0.) {
1000             dpchsp_singular(5, n-1);
1001             }
1002             (d_datap)[0+(__inc_d_n*(n-1))] = (g * (d_datap)[0+(__inc_d_n*(n-2))] + (d_datap)[0+(__inc_d_n*(n-1))]) / (dy_dx_datap)[0+(__inc_dy_dx_n*(n-1))];
1003             }
1004             }
1005             /* CARRY OUT BACK SUBSTITUTION */
1006             {/* Open n=nm1-1::-1 */ PDL_EXPAND2(register PDL_Indx n=PDLMIN(nm1-1, (__n_size-1)), __n_stop=0); for(; n>=__n_stop; n+=-1) {
1007             #line 5376 "lib/PDL/Primitive.pd"
1008             if ((dy_dx_datap)[0+(__inc_dy_dx_n*(n))] == 0.) {
1009             dpchsp_singular(6, n);
1010             }
1011             (d_datap)[0+(__inc_d_n*(n))] = ((d_datap)[0+(__inc_d_n*(n))] - (dx_datap)[0+(__inc_dx_n*(n))] * (d_datap)[0+(__inc_d_n*(n+1))]) / (dy_dx_datap)[0+(__inc_dy_dx_n*(n))];
1012             }} /* Close n=nm1-1::-1 */
1013             #line 5381 "lib/PDL/Primitive.pd"
1014             /* --------------------( END CODING FROM CUBSPL )-------------------- */
1015             #undef dpchsp_singular
1016             #line 1017 "lib/PDL/Primitive-pp-pchip_chsp.c"
1017 0 0         }PDL_BROADCASTLOOP_END_pchip_chsp_readdata
    0          
1018 0           } break;
1019 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_chsp: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
1020             }
1021             #undef PDL_IF_BAD
1022 3           return PDL_err;
1023             }
1024              
1025             static pdl_datatypes pdl_pchip_chsp_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
1026             static PDL_Indx pdl_pchip_chsp_vtable_realdims[] = { 1, 1, 1, 1, 1, 0, 1, 1 };
1027             static char *pdl_pchip_chsp_vtable_parnames[] = { "ic","vc","x","f","d","ierr","dx","dy_dx" };
1028             static short pdl_pchip_chsp_vtable_parflags[] = {
1029             PDL_PARAM_ISTYPED,
1030             0,
1031             0,
1032             0,
1033             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
1034             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE,
1035             PDL_PARAM_ISCREAT|PDL_PARAM_ISTEMP|PDL_PARAM_ISWRITE,
1036             PDL_PARAM_ISCREAT|PDL_PARAM_ISTEMP|PDL_PARAM_ISWRITE
1037             };
1038             static pdl_datatypes pdl_pchip_chsp_vtable_partypes[] = { PDL_SB, -1, -1, -1, -1, PDL_IND, -1, -1 };
1039             static PDL_Indx pdl_pchip_chsp_vtable_realdims_starts[] = { 0, 1, 2, 3, 4, 5, 5, 6 };
1040             static PDL_Indx pdl_pchip_chsp_vtable_realdims_ind_ids[] = { 1, 1, 0, 0, 0, 0, 0 };
1041             static char *pdl_pchip_chsp_vtable_indnames[] = { "n","two" };
1042             pdl_transvtable pdl_pchip_chsp_vtable = {
1043             PDL_TRANS_DO_BROADCAST, 0, pdl_pchip_chsp_vtable_gentypes, 4, 8, NULL /*CORE21*/,
1044             pdl_pchip_chsp_vtable_realdims, pdl_pchip_chsp_vtable_parnames,
1045             pdl_pchip_chsp_vtable_parflags, pdl_pchip_chsp_vtable_partypes,
1046             pdl_pchip_chsp_vtable_realdims_starts, pdl_pchip_chsp_vtable_realdims_ind_ids, 7,
1047             2, pdl_pchip_chsp_vtable_indnames,
1048             pdl_pchip_chsp_redodims, pdl_pchip_chsp_readdata, NULL,
1049             NULL,
1050             0,"PDL::Primitive::pchip_chsp"
1051             };
1052              
1053              
1054 3           pdl_error pdl_run_pchip_chsp(pdl *ic,pdl *vc,pdl *x,pdl *f,pdl *d,pdl *ierr) {
1055 3           pdl_error PDL_err = {0, NULL, 0};
1056 3 50         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
1057 3           pdl_trans *__privtrans = PDL->create_trans(&pdl_pchip_chsp_vtable);
1058 3 50         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
1059 3           __privtrans->pdls[0] = ic;
1060 3           __privtrans->pdls[1] = vc;
1061 3           __privtrans->pdls[2] = x;
1062 3           __privtrans->pdls[3] = f;
1063 3           __privtrans->pdls[4] = d;
1064 3           __privtrans->pdls[5] = ierr;
1065 3 50         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
1066 3 50         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
1067 3           return PDL_err;
1068             }