File Coverage

lib/PDL/Primitive-pp-interpolate.c
Criterion Covered Total %
statement 35 48 72.9
branch 37 190 19.4
condition n/a
subroutine n/a
pod n/a
total 72 238 30.2


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-interpolate.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_interpolate_readdata(pdl_trans *__privtrans) {
43             pdl_error PDL_err = {0, NULL, 0};
44             #line 45 "lib/PDL/Primitive-pp-interpolate.c"
45 7 50         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "broadcast.incs NULL");
46             /* broadcastloop declarations */
47             int __brcloopval;
48             register PDL_Indx __tind0,__tind1; /* counters along dim */
49 7           register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
50             /* dims here are how many steps along those dims */
51 7           register PDL_Indx __tinc0_xi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
52 7           register PDL_Indx __tinc0_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
53 7           register PDL_Indx __tinc0_y = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
54 7           register PDL_Indx __tinc0_yi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
55 7           register PDL_Indx __tinc0_err = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
56 7           register PDL_Indx __tinc1_xi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
57 7           register PDL_Indx __tinc1_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
58 7           register PDL_Indx __tinc1_y = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
59 7           register PDL_Indx __tinc1_yi = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
60 7           register PDL_Indx __tinc1_err = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
61             #define PDL_BROADCASTLOOP_START_interpolate_readdata PDL_BROADCASTLOOP_START( \
62             readdata, \
63             __privtrans->broadcast, \
64             __privtrans->vtable, \
65             xi_datap += __offsp[0]; \
66             x_datap += __offsp[1]; \
67             y_datap += __offsp[2]; \
68             yi_datap += __offsp[3]; \
69             err_datap += __offsp[4]; \
70             , \
71             ( ,xi_datap += __tinc1_xi - __tinc0_xi * __tdims0 \
72             ,x_datap += __tinc1_x - __tinc0_x * __tdims0 \
73             ,y_datap += __tinc1_y - __tinc0_y * __tdims0 \
74             ,yi_datap += __tinc1_yi - __tinc0_yi * __tdims0 \
75             ,err_datap += __tinc1_err - __tinc0_err * __tdims0 \
76             ), \
77             ( ,xi_datap += __tinc0_xi \
78             ,x_datap += __tinc0_x \
79             ,y_datap += __tinc0_y \
80             ,yi_datap += __tinc0_yi \
81             ,err_datap += __tinc0_err \
82             ) \
83             )
84             #define PDL_BROADCASTLOOP_END_interpolate_readdata PDL_BROADCASTLOOP_END( \
85             __privtrans->broadcast, \
86             xi_datap -= __tinc1_xi * __tdims1 + __offsp[0]; \
87             x_datap -= __tinc1_x * __tdims1 + __offsp[1]; \
88             y_datap -= __tinc1_y * __tdims1 + __offsp[2]; \
89             yi_datap -= __tinc1_yi * __tdims1 + __offsp[3]; \
90             err_datap -= __tinc1_err * __tdims1 + __offsp[4]; \
91             )
92 7           register PDL_Indx __inc_x_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_x_n;
93 7           register PDL_Indx __inc_y_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_y_n;
94             #ifndef PDL_DECLARE_PARAMS_interpolate_1
95             #define PDL_DECLARE_PARAMS_interpolate_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_xi,PDL_PPSYM_PARAM_xi,PDL_TYPE_PARAM_x,PDL_PPSYM_PARAM_x,PDL_TYPE_PARAM_err,PDL_PPSYM_PARAM_err) \
96             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_xi, xi, (__privtrans->pdls[0]), 1, PDL_PPSYM_PARAM_xi) \
97             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_x, x, (__privtrans->pdls[1]), 1, PDL_PPSYM_PARAM_x) \
98             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, y, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
99             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, yi, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
100             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_err, err, (__privtrans->pdls[4]), 1, PDL_PPSYM_PARAM_err)
101             #endif
102             #define PDL_IF_BAD(t,f) f
103 7           switch (__privtrans->__datatype) { /* Start generic switch */
104 0           case PDL_F: {
105 0 0         PDL_DECLARE_PARAMS_interpolate_1(PDL_Float,F,PDL_Float,F,PDL_Float,F,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
106             {
107             #line 2932 "lib/PDL/Primitive.pd"
108             PDL_Indx n = __privtrans->ind_sizes[0], n1 = n-1;
109             PDL_BROADCASTLOOP_START_interpolate_readdata
110             #line 2934 "lib/PDL/Primitive.pd"
111             PDL_Indx jl = -1, jh = n;
112             int carp = 0, up = ((x_datap)[0+(__inc_x_n*(n1))] > (x_datap)[0+(__inc_x_n*(0))]);
113             while (jh-jl > 1) { /* binary search */
114             PDL_Indx m = (jh+jl) >> 1;
115             if (((xi_datap)[0] > (x_datap)[0+(__inc_x_n*(m))]) == up)
116             jl = m;
117             else
118             jh = m;
119             }
120             if (jl == -1) {
121             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(0))]) carp = 1;
122             jl = 0;
123             } else if (jh == n) {
124             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(n1))]) carp = 1;
125             jl = n1-1;
126             }
127             jh = jl+1;
128             PDL_Float d = (x_datap)[0+(__inc_x_n*(jh))] - (x_datap)[0+(__inc_x_n*(jl))];
129             if (d == 0) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "identical abscissas");
130             d = ((x_datap)[0+(__inc_x_n*(jh))] - (xi_datap)[0])/d;
131             (yi_datap)[0] = d*(y_datap)[0+(__inc_y_n*(jl))] + (1-d)*(y_datap)[0+(__inc_y_n*(jh))];
132             (err_datap)[0] = carp;
133             PDL_BROADCASTLOOP_END_interpolate_readdata
134             #line 135 "lib/PDL/Primitive-pp-interpolate.c"
135             }
136 0           } break;
137 6           case PDL_D: {
138 6 50         PDL_DECLARE_PARAMS_interpolate_1(PDL_Double,D,PDL_Double,D,PDL_Double,D,PDL_Long,L)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
139             {
140             #line 2932 "lib/PDL/Primitive.pd"
141             PDL_Indx n = __privtrans->ind_sizes[0], n1 = n-1;
142             PDL_BROADCASTLOOP_START_interpolate_readdata
143             #line 2934 "lib/PDL/Primitive.pd"
144             PDL_Indx jl = -1, jh = n;
145             int carp = 0, up = ((x_datap)[0+(__inc_x_n*(n1))] > (x_datap)[0+(__inc_x_n*(0))]);
146             while (jh-jl > 1) { /* binary search */
147             PDL_Indx m = (jh+jl) >> 1;
148             if (((xi_datap)[0] > (x_datap)[0+(__inc_x_n*(m))]) == up)
149             jl = m;
150             else
151             jh = m;
152             }
153             if (jl == -1) {
154             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(0))]) carp = 1;
155             jl = 0;
156             } else if (jh == n) {
157             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(n1))]) carp = 1;
158             jl = n1-1;
159             }
160             jh = jl+1;
161             PDL_Double d = (x_datap)[0+(__inc_x_n*(jh))] - (x_datap)[0+(__inc_x_n*(jl))];
162             if (d == 0) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "identical abscissas");
163             d = ((x_datap)[0+(__inc_x_n*(jh))] - (xi_datap)[0])/d;
164             (yi_datap)[0] = d*(y_datap)[0+(__inc_y_n*(jl))] + (1-d)*(y_datap)[0+(__inc_y_n*(jh))];
165             (err_datap)[0] = carp;
166             PDL_BROADCASTLOOP_END_interpolate_readdata
167             #line 168 "lib/PDL/Primitive-pp-interpolate.c"
168             }
169 4           } break;
170 0           case PDL_LD: {
171 0 0         PDL_DECLARE_PARAMS_interpolate_1(PDL_LDouble,E,PDL_LDouble,E,PDL_LDouble,E,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
172             {
173             #line 2932 "lib/PDL/Primitive.pd"
174             PDL_Indx n = __privtrans->ind_sizes[0], n1 = n-1;
175             PDL_BROADCASTLOOP_START_interpolate_readdata
176             #line 2934 "lib/PDL/Primitive.pd"
177             PDL_Indx jl = -1, jh = n;
178             int carp = 0, up = ((x_datap)[0+(__inc_x_n*(n1))] > (x_datap)[0+(__inc_x_n*(0))]);
179             while (jh-jl > 1) { /* binary search */
180             PDL_Indx m = (jh+jl) >> 1;
181             if (((xi_datap)[0] > (x_datap)[0+(__inc_x_n*(m))]) == up)
182             jl = m;
183             else
184             jh = m;
185             }
186             if (jl == -1) {
187             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(0))]) carp = 1;
188             jl = 0;
189             } else if (jh == n) {
190             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(n1))]) carp = 1;
191             jl = n1-1;
192             }
193             jh = jl+1;
194             PDL_LDouble d = (x_datap)[0+(__inc_x_n*(jh))] - (x_datap)[0+(__inc_x_n*(jl))];
195             if (d == 0) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "identical abscissas");
196             d = ((x_datap)[0+(__inc_x_n*(jh))] - (xi_datap)[0])/d;
197             (yi_datap)[0] = d*(y_datap)[0+(__inc_y_n*(jl))] + (1-d)*(y_datap)[0+(__inc_y_n*(jh))];
198             (err_datap)[0] = carp;
199             PDL_BROADCASTLOOP_END_interpolate_readdata
200             #line 201 "lib/PDL/Primitive-pp-interpolate.c"
201             }
202 0           } break;
203 0           case PDL_CF: {
204 0 0         PDL_DECLARE_PARAMS_interpolate_1(PDL_CFloat,G,PDL_Float,F,PDL_Float,F,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
205             {
206             #line 2932 "lib/PDL/Primitive.pd"
207             PDL_Indx n = __privtrans->ind_sizes[0], n1 = n-1;
208             PDL_BROADCASTLOOP_START_interpolate_readdata
209             #line 2934 "lib/PDL/Primitive.pd"
210             PDL_Indx jl = -1, jh = n;
211             int carp = 0, up = ((x_datap)[0+(__inc_x_n*(n1))] > (x_datap)[0+(__inc_x_n*(0))]);
212             while (jh-jl > 1) { /* binary search */
213             PDL_Indx m = (jh+jl) >> 1;
214             if (((xi_datap)[0] > (x_datap)[0+(__inc_x_n*(m))]) == up)
215             jl = m;
216             else
217             jh = m;
218             }
219             if (jl == -1) {
220             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(0))]) carp = 1;
221             jl = 0;
222             } else if (jh == n) {
223             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(n1))]) carp = 1;
224             jl = n1-1;
225             }
226             jh = jl+1;
227             PDL_CFloat d = (x_datap)[0+(__inc_x_n*(jh))] - (x_datap)[0+(__inc_x_n*(jl))];
228             if (d == 0) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "identical abscissas");
229             d = ((x_datap)[0+(__inc_x_n*(jh))] - (xi_datap)[0])/d;
230             (yi_datap)[0] = d*(y_datap)[0+(__inc_y_n*(jl))] + (1-d)*(y_datap)[0+(__inc_y_n*(jh))];
231             (err_datap)[0] = carp;
232             PDL_BROADCASTLOOP_END_interpolate_readdata
233             #line 234 "lib/PDL/Primitive-pp-interpolate.c"
234             }
235 0           } break;
236 1           case PDL_CD: {
237 1 50         PDL_DECLARE_PARAMS_interpolate_1(PDL_CDouble,C,PDL_Double,D,PDL_Double,D,PDL_Long,L)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
238             {
239             #line 2932 "lib/PDL/Primitive.pd"
240             PDL_Indx n = __privtrans->ind_sizes[0], n1 = n-1;
241             PDL_BROADCASTLOOP_START_interpolate_readdata
242             #line 2934 "lib/PDL/Primitive.pd"
243             PDL_Indx jl = -1, jh = n;
244             int carp = 0, up = ((x_datap)[0+(__inc_x_n*(n1))] > (x_datap)[0+(__inc_x_n*(0))]);
245             while (jh-jl > 1) { /* binary search */
246             PDL_Indx m = (jh+jl) >> 1;
247             if (((xi_datap)[0] > (x_datap)[0+(__inc_x_n*(m))]) == up)
248             jl = m;
249             else
250             jh = m;
251             }
252             if (jl == -1) {
253             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(0))]) carp = 1;
254             jl = 0;
255             } else if (jh == n) {
256             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(n1))]) carp = 1;
257             jl = n1-1;
258             }
259             jh = jl+1;
260             PDL_CDouble d = (x_datap)[0+(__inc_x_n*(jh))] - (x_datap)[0+(__inc_x_n*(jl))];
261             if (d == 0) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "identical abscissas");
262             d = ((x_datap)[0+(__inc_x_n*(jh))] - (xi_datap)[0])/d;
263             (yi_datap)[0] = d*(y_datap)[0+(__inc_y_n*(jl))] + (1-d)*(y_datap)[0+(__inc_y_n*(jh))];
264             (err_datap)[0] = carp;
265             PDL_BROADCASTLOOP_END_interpolate_readdata
266             #line 267 "lib/PDL/Primitive-pp-interpolate.c"
267             }
268 1           } break;
269 0           case PDL_CLD: {
270 0 0         PDL_DECLARE_PARAMS_interpolate_1(PDL_CLDouble,H,PDL_LDouble,E,PDL_LDouble,E,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
271             {
272             #line 2932 "lib/PDL/Primitive.pd"
273             PDL_Indx n = __privtrans->ind_sizes[0], n1 = n-1;
274             PDL_BROADCASTLOOP_START_interpolate_readdata
275             #line 2934 "lib/PDL/Primitive.pd"
276             PDL_Indx jl = -1, jh = n;
277             int carp = 0, up = ((x_datap)[0+(__inc_x_n*(n1))] > (x_datap)[0+(__inc_x_n*(0))]);
278             while (jh-jl > 1) { /* binary search */
279             PDL_Indx m = (jh+jl) >> 1;
280             if (((xi_datap)[0] > (x_datap)[0+(__inc_x_n*(m))]) == up)
281             jl = m;
282             else
283             jh = m;
284             }
285             if (jl == -1) {
286             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(0))]) carp = 1;
287             jl = 0;
288             } else if (jh == n) {
289             if ((xi_datap)[0] != (x_datap)[0+(__inc_x_n*(n1))]) carp = 1;
290             jl = n1-1;
291             }
292             jh = jl+1;
293             PDL_CLDouble d = (x_datap)[0+(__inc_x_n*(jh))] - (x_datap)[0+(__inc_x_n*(jl))];
294             if (d == 0) return PDL->make_error(PDL_EUSERERROR, "Error in interpolate:" "identical abscissas");
295             d = ((x_datap)[0+(__inc_x_n*(jh))] - (xi_datap)[0])/d;
296             (yi_datap)[0] = d*(y_datap)[0+(__inc_y_n*(jl))] + (1-d)*(y_datap)[0+(__inc_y_n*(jh))];
297             (err_datap)[0] = carp;
298             PDL_BROADCASTLOOP_END_interpolate_readdata
299             #line 300 "lib/PDL/Primitive-pp-interpolate.c"
300             }
301 0           } break;
302 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in interpolate: unhandled datatype(%d), only handles (FDEGCH)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
303             }
304             #undef PDL_IF_BAD
305 5           return PDL_err;
306             }
307              
308             static pdl_datatypes pdl_interpolate_vtable_gentypes[] = { PDL_F, PDL_LD, PDL_CF, PDL_CD, PDL_CLD, PDL_D, -1 };
309             static PDL_Indx pdl_interpolate_vtable_realdims[] = { 0, 1, 1, 0, 0 };
310             static char *pdl_interpolate_vtable_parnames[] = { "xi","x","y","yi","err" };
311             static short pdl_interpolate_vtable_parflags[] = {
312             PDL_PARAM_ISNOTCOMPLEX,
313             PDL_PARAM_ISNOTCOMPLEX,
314             0,
315             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
316             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE
317             };
318             static pdl_datatypes pdl_interpolate_vtable_partypes[] = { -1, -1, -1, -1, PDL_L };
319             static PDL_Indx pdl_interpolate_vtable_realdims_starts[] = { 0, 0, 1, 2, 2 };
320             static PDL_Indx pdl_interpolate_vtable_realdims_ind_ids[] = { 0, 0 };
321             static char *pdl_interpolate_vtable_indnames[] = { "n" };
322             pdl_transvtable pdl_interpolate_vtable = {
323             PDL_TRANS_DO_BROADCAST|PDL_TRANS_BADIGNORE, 0, pdl_interpolate_vtable_gentypes, 3, 5, NULL /*CORE21*/,
324             pdl_interpolate_vtable_realdims, pdl_interpolate_vtable_parnames,
325             pdl_interpolate_vtable_parflags, pdl_interpolate_vtable_partypes,
326             pdl_interpolate_vtable_realdims_starts, pdl_interpolate_vtable_realdims_ind_ids, 2,
327             1, pdl_interpolate_vtable_indnames,
328             NULL, pdl_interpolate_readdata, NULL,
329             NULL,
330             0,"PDL::Primitive::interpolate"
331             };
332              
333              
334 8           pdl_error pdl_run_interpolate(pdl *xi,pdl *x,pdl *y,pdl *yi,pdl *err) {
335 8           pdl_error PDL_err = {0, NULL, 0};
336 8 50         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
337 8           pdl_trans *__privtrans = PDL->create_trans(&pdl_interpolate_vtable);
338 8 50         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
339 8           __privtrans->pdls[0] = xi;
340 8           __privtrans->pdls[1] = x;
341 8           __privtrans->pdls[2] = y;
342 8           __privtrans->pdls[3] = yi;
343 8           __privtrans->pdls[4] = err;
344 8 100         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
345 5 100         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
346 4           return PDL_err;
347             }