File Coverage

lib/PDL/Primitive-pp-pchip_chic.c

Criterion	Covered	Total	%
statement	60	80	75.0
branch	54	288	18.7
condition			n/a
subroutine			n/a
pod			n/a
total	114	368	30.9

line	stmt	bran	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_chic.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_chic_redodims(pdl_trans *__privtrans) {
43			pdl_error PDL_err = {0, NULL, 0};
44			#line 45 "lib/PDL/Primitive-pp-pchip_chic.c"
45	3		__privtrans->ind_sizes[1] = __privtrans->ind_sizes[0]-1;
46	3		__privtrans->ind_sizes[2] = 2;
47			#ifndef PDL_DECLARE_PARAMS_pchip_chic_0
48			#define PDL_DECLARE_PARAMS_pchip_chic_0(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ic,PDL_PPSYM_PARAM_ic,PDL_TYPE_PARAM_ierr,PDL_PPSYM_PARAM_ierr) \
49			PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ic, ic, (__privtrans->pdls[0]), 0, PDL_PPSYM_PARAM_ic) \
50			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, vc, (__privtrans->pdls[1]), 0, PDL_PPSYM_OP) \
51			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, mflag, (__privtrans->pdls[2]), 0, PDL_PPSYM_OP) \
52			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[3]), 0, PDL_PPSYM_OP) \
53			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, f, (__privtrans->pdls[4]), 0, PDL_PPSYM_OP) \
54			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, d, (__privtrans->pdls[5]), 0, PDL_PPSYM_OP) \
55			PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ierr, ierr, (__privtrans->pdls[6]), 0, PDL_PPSYM_PARAM_ierr) \
56			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, h, (__privtrans->pdls[7]), 0, PDL_PPSYM_OP) \
57			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, slope, (__privtrans->pdls[8]), 0, PDL_PPSYM_OP)
58			#endif
59			#define PDL_IF_BAD(t,f) f
60	3		switch (__privtrans->__datatype) { /* Start generic switch */
61	0		case PDL_F: {
62	0	0	PDL_DECLARE_PARAMS_pchip_chic_0(PDL_Float,F,PDL_SByte,A,PDL_Indx,N)
		0
		0
		0
		0
		0
		0
		0
		0
63	0	0	{if (__privtrans->ind_sizes[0] < 2) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "NUMBER OF DATA POINTS LESS THAN TWO");
64			}
65	0		} break;
66	3		case PDL_D: {
67	3	50	PDL_DECLARE_PARAMS_pchip_chic_0(PDL_Double,D,PDL_SByte,A,PDL_Indx,N)
		50
		50
		50
		50
		50
		50
		50
		50
68	3	50	{if (__privtrans->ind_sizes[0] < 2) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "NUMBER OF DATA POINTS LESS THAN TWO");
69			}
70	3		} break;
71	0		case PDL_LD: {
72	0	0	PDL_DECLARE_PARAMS_pchip_chic_0(PDL_LDouble,E,PDL_SByte,A,PDL_Indx,N)
		0
		0
		0
		0
		0
		0
		0
		0
73	0	0	{if (__privtrans->ind_sizes[0] < 2) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "NUMBER OF DATA POINTS LESS THAN TWO");
74			}
75	0		} break;
76	0		default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_chic: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
77			}
78			#undef PDL_IF_BAD
79
80	3	50	PDL_RETERROR(PDL_err, PDL->redodims_default(__privtrans));
81	3		return PDL_err;
82			}
83
84
85			#line 1857 "lib/PDL/PP.pm"
86			pdl_error pdl_pchip_chic_readdata(pdl_trans *__privtrans) {
87			pdl_error PDL_err = {0, NULL, 0};
88			#line 89 "lib/PDL/Primitive-pp-pchip_chic.c"
89	3		register PDL_Indx __n_size = __privtrans->ind_sizes[0];
90	3		register PDL_Indx __nless1_size = __privtrans->ind_sizes[1];
91	3	50	if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "broadcast.incs NULL");
92			/* broadcastloop declarations */
93			int __brcloopval;
94			register PDL_Indx __tind0,__tind1; /* counters along dim */
95	3		register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
96			/* dims here are how many steps along those dims */
97	3		register PDL_Indx __tinc0_ic = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
98	3		register PDL_Indx __tinc0_vc = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
99	3		register PDL_Indx __tinc0_mflag = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
100	3		register PDL_Indx __tinc0_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,0);
101	3		register PDL_Indx __tinc0_f = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,0);
102	3		register PDL_Indx __tinc0_d = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,0);
103	3		register PDL_Indx __tinc0_ierr = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,0);
104	3		register PDL_Indx __tinc0_h = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,7,0);
105	3		register PDL_Indx __tinc0_slope = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,8,0);
106	3		register PDL_Indx __tinc1_ic = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
107	3		register PDL_Indx __tinc1_vc = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
108	3		register PDL_Indx __tinc1_mflag = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
109	3		register PDL_Indx __tinc1_x = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,3,1);
110	3		register PDL_Indx __tinc1_f = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,4,1);
111	3		register PDL_Indx __tinc1_d = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,5,1);
112	3		register PDL_Indx __tinc1_ierr = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,6,1);
113	3		register PDL_Indx __tinc1_h = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,7,1);
114	3		register PDL_Indx __tinc1_slope = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,8,1);
115			#define PDL_BROADCASTLOOP_START_pchip_chic_readdata PDL_BROADCASTLOOP_START( \
116			readdata, \
117			__privtrans->broadcast, \
118			__privtrans->vtable, \
119			ic_datap += __offsp[0]; \
120			vc_datap += __offsp[1]; \
121			mflag_datap += __offsp[2]; \
122			x_datap += __offsp[3]; \
123			f_datap += __offsp[4]; \
124			d_datap += __offsp[5]; \
125			ierr_datap += __offsp[6]; \
126			h_datap += __offsp[7]; \
127			slope_datap += __offsp[8]; \
128			, \
129			( ,ic_datap += __tinc1_ic - __tinc0_ic * __tdims0 \
130			,vc_datap += __tinc1_vc - __tinc0_vc * __tdims0 \
131			,mflag_datap += __tinc1_mflag - __tinc0_mflag * __tdims0 \
132			,x_datap += __tinc1_x - __tinc0_x * __tdims0 \
133			,f_datap += __tinc1_f - __tinc0_f * __tdims0 \
134			,d_datap += __tinc1_d - __tinc0_d * __tdims0 \
135			,ierr_datap += __tinc1_ierr - __tinc0_ierr * __tdims0 \
136			,h_datap += __tinc1_h - __tinc0_h * __tdims0 \
137			,slope_datap += __tinc1_slope - __tinc0_slope * __tdims0 \
138			), \
139			( ,ic_datap += __tinc0_ic \
140			,vc_datap += __tinc0_vc \
141			,mflag_datap += __tinc0_mflag \
142			,x_datap += __tinc0_x \
143			,f_datap += __tinc0_f \
144			,d_datap += __tinc0_d \
145			,ierr_datap += __tinc0_ierr \
146			,h_datap += __tinc0_h \
147			,slope_datap += __tinc0_slope \
148			) \
149			)
150			#define PDL_BROADCASTLOOP_END_pchip_chic_readdata PDL_BROADCASTLOOP_END( \
151			__privtrans->broadcast, \
152			ic_datap -= __tinc1_ic * __tdims1 + __offsp[0]; \
153			vc_datap -= __tinc1_vc * __tdims1 + __offsp[1]; \
154			mflag_datap -= __tinc1_mflag * __tdims1 + __offsp[2]; \
155			x_datap -= __tinc1_x * __tdims1 + __offsp[3]; \
156			f_datap -= __tinc1_f * __tdims1 + __offsp[4]; \
157			d_datap -= __tinc1_d * __tdims1 + __offsp[5]; \
158			ierr_datap -= __tinc1_ierr * __tdims1 + __offsp[6]; \
159			h_datap -= __tinc1_h * __tdims1 + __offsp[7]; \
160			slope_datap -= __tinc1_slope * __tdims1 + __offsp[8]; \
161			)
162	3		register PDL_Indx __inc_d_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,5,0)]; (void)__inc_d_n;
163	3		register PDL_Indx __inc_f_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,4,0)]; (void)__inc_f_n;
164	3		register PDL_Indx __inc_h_nless1 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,7,0)]; (void)__inc_h_nless1;
165	3		register PDL_Indx __inc_ic_two = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_ic_two;
166	3		register PDL_Indx __inc_slope_nless1 = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,8,0)]; (void)__inc_slope_nless1;
167	3		register PDL_Indx __inc_vc_two = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_vc_two;
168	3		register PDL_Indx __inc_x_n = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,3,0)]; (void)__inc_x_n;
169			#ifndef PDL_DECLARE_PARAMS_pchip_chic_1
170			#define PDL_DECLARE_PARAMS_pchip_chic_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_ic,PDL_PPSYM_PARAM_ic,PDL_TYPE_PARAM_ierr,PDL_PPSYM_PARAM_ierr) \
171			PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ic, ic, (__privtrans->pdls[0]), 1, PDL_PPSYM_PARAM_ic) \
172			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, vc, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
173			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, mflag, (__privtrans->pdls[2]), 1, PDL_PPSYM_OP) \
174			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, x, (__privtrans->pdls[3]), 1, PDL_PPSYM_OP) \
175			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, f, (__privtrans->pdls[4]), 1, PDL_PPSYM_OP) \
176			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, d, (__privtrans->pdls[5]), 1, PDL_PPSYM_OP) \
177			PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_ierr, ierr, (__privtrans->pdls[6]), 1, PDL_PPSYM_PARAM_ierr) \
178			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, h, (__privtrans->pdls[7]), 1, PDL_PPSYM_OP) \
179			PDL_DECLARE_PARAMETER(PDL_TYPE_OP, slope, (__privtrans->pdls[8]), 1, PDL_PPSYM_OP)
180			#endif
181			#define PDL_IF_BAD(t,f) f
182	3		switch (__privtrans->__datatype) { /* Start generic switch */
183	0		case PDL_F: {
184	0	0	PDL_DECLARE_PARAMS_pchip_chic_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
		0
		0
		0
185	0	0	PDL_BROADCASTLOOP_START_pchip_chic_readdata {
		0
		0
		0
		0
		0
		0
186			#line 4574 "lib/PDL/Primitive.pd"
187			const PDL_Float d1mach = FLT_EPSILON;
188			/* VALIDITY-CHECK ARGUMENTS. */
189			{/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
190			#line 4577 "lib/PDL/Primitive.pd"
191			if ((x_datap)[0+(__inc_x_n(n))] > (x_datap)[0+(__inc_x_n(n-1))]) continue;
192			(ierr_datap)[0] = -1;
193			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "X-ARRAY NOT STRICTLY INCREASING");
194			}} /* Close n=1 */
195			#line 4581 "lib/PDL/Primitive.pd"
196			PDL_Indx ibeg = (ic_datap)[0+(__inc_ic_two(0))], iend = (ic_datap)[0+(__inc_ic_two(1))], n = __privtrans->ind_sizes[0];
197			(ierr_datap)[0] = 0;
198			if (PDL_ABS(ibeg) > 5)
199			--((ierr_datap)[0]);
200			if (PDL_ABS(iend) > 5)
201			(ierr_datap)[0] += -2;
202			if ((ierr_datap)[0] < 0) {
203			(ierr_datap)[0] += -3;
204			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "IC OUT OF RANGE");
205			}
206			/* FUNCTION DEFINITION IS OK -- GO ON. */
207			/* SET UP H AND SLOPE ARRAYS. */
208			{/* Open nless1 */ PDL_EXPAND2(register PDL_Indx nless1=0, __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
209			#line 4594 "lib/PDL/Primitive.pd"
210			(h_datap)[0+(__inc_h_nless1(nless1))] = (x_datap)[0+(__inc_x_n(nless1+1))] - (x_datap)[0+(__inc_x_n*(nless1))];
211			(slope_datap)[0+(__inc_slope_nless1(nless1))] = (f_datap)[0+(__inc_f_n(nless1+1))] - (f_datap)[0+(__inc_f_n*(nless1))];
212			}} /* Close nless1 */
213			#line 4597 "lib/PDL/Primitive.pd"
214			/* SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. */
215			if (__privtrans->ind_sizes[1] <= 1) {
216			(d_datap)[0+(__inc_d_n(0))] = (d_datap)[0+(__inc_d_n(1))] = (slope_datap)[0+(__inc_slope_nless1*(0))];
217			} else {
218			/* NORMAL CASE (N .GE. 3) . */
219			/* SET INTERIOR DERIVATIVES AND DEFAULT END CONDITIONS. */
220			do { /* inline dpchci */
221			/* Local variables */
222			PDL_Float del1 = (slope_datap)[0+(__inc_slope_nless1*(0))];
223			/* SPECIAL CASE N=2 is dealt with in separate branch above */
224			/* NORMAL CASE (N .GE. 3). */
225			PDL_Float del2 = (slope_datap)[0+(__inc_slope_nless1*(1))];
226			/* SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */
227			/* SHAPE-PRESERVING. */
228			PDL_Float hsum = (h_datap)[0+(__inc_h_nless1(0))] + (h_datap)[0+(__inc_h_nless1(1))];
229			PDL_Float w1 = ((h_datap)[0+(__inc_h_nless1*(0))] + hsum) / hsum;
230			PDL_Float w2 = -(h_datap)[0+(__inc_h_nless1*(0))] / hsum;
231			(d_datap)[0+(__inc_d_n(0))] = w1 del1 + w2 * del2;
232			if (((((d_datap)[0+(__inc_d_n(0))]) == 0. \|\| (del1) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.) {
233			(d_datap)[0+(__inc_d_n*(0))] = 0.;
234			} else if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
235			/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */
236			PDL_Float dmax = 3. * del1;
237			if (PDL_ABS((d_datap)[0+(__inc_d_n*(0))]) > PDL_ABS(dmax))
238			(d_datap)[0+(__inc_d_n*(0))] = dmax;
239			}
240			/* LOOP THROUGH INTERIOR POINTS. */
241			{/* Open nless1=1 */ PDL_EXPAND2(register PDL_Indx nless1=PDLMAX((1),0), __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
242			#line 4625 "lib/PDL/Primitive.pd"
243			if (nless1 != 1) {
244			hsum = (h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))];
245			del1 = del2;
246			del2 = (slope_datap)[0+(__inc_slope_nless1*(nless1))];
247			}
248			/* SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. */
249			(d_datap)[0+(__inc_d_n*(nless1))] = 0.;
250			if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.)
251			continue;
252			/* USE BRODLIE MODIFICATION OF BUTLAND FORMULA. */
253			PDL_Float hsumt3 = hsum + hsum + hsum;
254			w1 = (hsum + (h_datap)[0+(__inc_h_nless1*(nless1-1))]) / hsumt3;
255			w2 = (hsum + (h_datap)[0+(__inc_h_nless1*(nless1))]) / hsumt3;
256			/* Computing MAX */
257			PDL_Float dmax = PDLMAX(PDL_ABS(del1),PDL_ABS(del2));
258			/* Computing MIN */
259			PDL_Float dmin = PDLMIN(PDL_ABS(del1),PDL_ABS(del2));
260			PDL_Float drat1 = del1 / dmax, drat2 = del2 / dmax;
261			(d_datap)[0+(__inc_d_n(nless1))] = dmin / (w1 drat1 + w2 * drat2);
262			}} /* Close nless1=1 */
263			#line 4645 "lib/PDL/Primitive.pd"
264			/* SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */
265			/* SHAPE-PRESERVING. */
266			w1 = -(h_datap)[0+(__inc_h_nless1*(n-2))] / hsum;
267			w2 = ((h_datap)[0+(__inc_h_nless1*(n-2))] + hsum) / hsum;
268			(d_datap)[0+(__inc_d_n(n-1))] = w1 del1 + w2 * del2;
269			if (((((d_datap)[0+(__inc_d_n(n-1))]) == 0. \|\| (del2) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(n-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.) {
270			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
271			} else if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
272			/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */
273			PDL_Float dmax = 3. * del2;
274			if (PDL_ABS((d_datap)[0+(__inc_d_n*(n-1))]) > PDL_ABS(dmax))
275			(d_datap)[0+(__inc_d_n*(n-1))] = dmax;
276			}
277			} while (0); /* end inline dpchci */
278			/* SET DERIVATIVES AT POINTS WHERE MONOTONICITY SWITCHES DIRECTION. */
279			if ((mflag_datap)[0] != 0.) {
280			do { /* inline dpchcs */
281			/* **PURPOSE Adjusts derivative values for DPCHIC /
282			/* DPCHCS: DPCHIC Monotonicity Switch Derivative Setter. */
283			/* Called by DPCHIC to adjust the values of D in the vicinity of a */
284			/* switch in direction of monotonicity, to produce a more "visually */
285			/* pleasing" curve than that given by DPCHIM . */
286			static const PDL_Float fudge = 4.;
287			/* Local variables */
288			PDL_Indx k;
289			PDL_Float del[3], fact, dfmx;
290			PDL_Float dext, dfloc, slmax, wtave[2];
291			/* INITIALIZE. */
292			/* LOOP OVER SEGMENTS. */
293			{/* Open nless1=1 */ PDL_EXPAND2(register PDL_Indx nless1=PDLMAX((1),0), __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
294			#line 4675 "lib/PDL/Primitive.pd"
295			PDL_Float dtmp = ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)));
296			if (dtmp > 0.) {
297			continue;
298			}
299			if (dtmp != 0.) {
300			/* ....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT ..................... */
301			/* DO NOT CHANGE D IF 'UP-DOWN-UP'. */
302			if (nless1 > 1) {
303			if (((((slope_datap)[0+(__inc_slope_nless1(nless1-2))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-2))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) > 0.)
304			continue;
305			/* -------------------------- */
306			}
307			if (nless1 < __privtrans->ind_sizes[1]-1 && ((((slope_datap)[0+(__inc_slope_nless1(nless1+1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1+1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) > 0.)
308			continue;
309			/* ---------------------------- */
310			/* ....... COMPUTE PROVISIONAL VALUE FOR D(1,I). */
311			dext = (h_datap)[0+(__inc_h_nless1(nless1))] / ((h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))]) (slope_datap)[0+(__inc_slope_nless1*(nless1-1))] +
312			(h_datap)[0+(__inc_h_nless1(nless1-1))] / ((h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))]) (slope_datap)[0+(__inc_slope_nless1*(nless1))];
313			/* ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM. */
314			dtmp = (((dext) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0.) ? 0. : (((dext)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)));
315			if (dtmp == 0) {
316			continue;
317			}
318			if (dtmp < 0.) {
319			/* DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS -- */
320			/* EXTREMUM IS IN (X(I-1),X(I)). */
321			k = nless1;
322			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I). */
323			wtave[1] = dext;
324			if (k > 1) {
325			wtave[0] = (h_datap)[0+(__inc_h_nless1(k-1))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) (slope_datap)[0+(__inc_slope_nless1*(k-2))] +
326			(h_datap)[0+(__inc_h_nless1(k-2))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1*(k-1))];
327			}
328			} else {
329			/* DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS -- */
330			/* EXTREMUM IS IN (X(I),X(I+1)). */
331			k = nless1 + 1;
332			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). */
333			wtave[0] = dext;
334			if (k < nless1) {
335			wtave[1] = (h_datap)[0+(__inc_h_nless1(k))] / ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k-1))]
336			/ ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) * (slope_datap)[0+(__inc_slope_nless1*(k))];
337			}
338			}
339			} else {
340			/* ....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO -- */
341			/* CHECK FOR FLAT-TOPPED PEAK ....................... */
342			if (nless1 == __privtrans->ind_sizes[1]-1 \|\| ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1+1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1+1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) >= 0.)
343			continue;
344			/* ----------------------------- */
345			/* WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)). */
346			k = nless1+1;
347			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). */
348			wtave[0] = (h_datap)[0+(__inc_h_nless1(k-1))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) (slope_datap)[0+(__inc_slope_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-2))]
349			/ ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) * (slope_datap)[0+(__inc_slope_nless1*(k-1))];
350			wtave[1] = (h_datap)[0+(__inc_h_nless1(k))] / ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k-1))] / (
351			(h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) * (slope_datap)[0+(__inc_slope_nless1*(k))];
352			}
353			/* ....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM */
354			/* ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE-- */
355			/* WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K), */
356			/* IF K.GT.1 */
357			/* WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1), */
358			/* IF K.LT.N-1 */
359			slmax = PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-1))]);
360			if (k > 1) {
361			/* Computing MAX */
362			slmax = PDLMAX(slmax,PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-2))]));
363			}
364			if (k < nless1) {
365			/* Computing MAX */
366			slmax = PDLMAX(slmax,PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k))]));
367			}
368			if (k > 1) {
369			del[0] = (slope_datap)[0+(__inc_slope_nless1*(k-2))] / slmax;
370			}
371			del[1] = (slope_datap)[0+(__inc_slope_nless1*(k-1))] / slmax;
372			if (k < nless1) {
373			del[2] = (slope_datap)[0+(__inc_slope_nless1*(k))] / slmax;
374			}
375			if (k > 1 && k < nless1) {
376			/* NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL. */
377			fact = fudge * PDL_ABS(del[2] * (del[0] - del[1]) * (wtave[1] / slmax));
378			(d_datap)[0+(__inc_d_n(k-1))] += PDLMIN(fact,1.) (wtave[0] - (d_datap)[0+(__inc_d_n*(k-1))]);
379			fact = fudge * PDL_ABS(del[0] * (del[2] - del[1]) * (wtave[0] / slmax));
380			(d_datap)[0+(__inc_d_n(k))] += PDLMIN(fact,1.) (wtave[1] - (d_datap)[0+(__inc_d_n*(k))]);
381			} else {
382			/* SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR */
383			/* K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1). */
384			fact = fudge * PDL_ABS(del[1]);
385			(d_datap)[0+(__inc_d_n(nless1))] = PDLMIN(fact,1.) wtave[nless1+1 - k];
386			/* NOTE THAT I-K+1 = 1 IF K=I (=NLESS1), */
387			/* I-K+1 = 2 IF K=I-1(=1). */
388			}
389			/* ....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA. */
390			if ((mflag_datap)[0] <= 0.) {
391			continue;
392			}
393			dfloc = (h_datap)[0+(__inc_h_nless1(k-1))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-1))]);
394			if (k > 1) {
395			/* Computing MAX */
396			dfloc = PDLMAX(dfloc,(h_datap)[0+(__inc_h_nless1(k-2))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-2))]));
397			}
398			if (k < nless1) {
399			/* Computing MAX */
400			dfloc = PDLMAX(dfloc,(h_datap)[0+(__inc_h_nless1(k))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k))]));
401			}
402			dfmx = (mflag_datap)[0] * dfloc;
403			PDL_Indx indx = nless1 - k;
404			/* INDX = 1 IF K=I, 2 IF K=I-1. */
405			/* --------------------------------------------------------------- */
406			do { /* inline dpchsw */
407			/* NOTATION AND GENERAL REMARKS. */
408			/* RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED. */
409			/* LAMBDA IS THE RATIO OF D2 TO D1. */
410			/* THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM. */
411			/* PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO), */
412			/* WHERE THAT = (XHAT - X1)/H . */
413			/* THAT IS, P(XHAT)-F1 = DHPHI, WHERE D=D1 OR D2. */
414			/* SIMILARLY, P(XHAT)-F2 = DH(PHI-RHO) . */
415			/* Local variables */
416			PDL_Float cp, nu, phi, rho, hphi, that, sigma, small;
417			PDL_Float lambda, radcal;
418			PDL_Float d1 = (d_datap)[0+(__inc_d_n(k-1))], d2 = (d_datap)[0+(__inc_d_n(k))], h2 = (h_datap)[0+(__inc_h_nless1(k-1))], slope2 = (slope_datap)[0+(__inc_slope_nless1(k-1))];
419			/* Initialized data */
420			static const PDL_Float fact = 100.;
421			/* THIRD SHOULD BE SLIGHTLY LESS THAN 1/3. */
422			static const PDL_Float third = .33333;
423			/* SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS. */
424			small = fact * d1mach;
425			/* DO MAIN CALCULATION. */
426			if (d1 == 0.) {
427			/* SPECIAL CASE -- D1.EQ.ZERO . */
428			/* IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED. */
429			if (d2 == 0.) {
430			(ierr_datap)[0] = -1;
431			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "D1 AND/OR D2 INVALID");
432			}
433			rho = slope2 / d2;
434			/* EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . */
435			if (rho >= third) {
436			(ierr_datap)[0] = 0; break;
437			}
438			that = 2. * (3. * rho - 1.) / (3. * (2. * rho - 1.));
439			/* Computing 2nd power */
440			phi = that * that * ((3. * rho - 1.) / 3.);
441			/* CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . */
442			if (indx != 3) {
443			phi -= rho;
444			}
445			/* TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. */
446			hphi = h2 * PDL_ABS(phi);
447			if (hphi * PDL_ABS(d2) > dfmx) {
448			/* AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK. */
449			d2 = ((d2) >= 0 ? PDL_ABS(dfmx / hphi) : -PDL_ABS(dfmx / hphi));
450			}
451			} else {
452			rho = slope2 / d1;
453			lambda = -(d2) / d1;
454			if (d2 == 0.) {
455			/* SPECIAL CASE -- D2.EQ.ZERO . */
456			/* EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . */
457			if (rho >= third) {
458			(ierr_datap)[0] = 0; break;
459			}
460			cp = 2. - 3. * rho;
461			nu = 1. - 2. * rho;
462			that = 1. / (3. * nu);
463			} else {
464			if (lambda <= 0.) {
465			(ierr_datap)[0] = -1;
466			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "D1 AND/OR D2 INVALID");
467			}
468			/* NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS. */
469			nu = 1. - lambda - 2. * rho;
470			sigma = 1. - rho;
471			cp = nu + sigma;
472			if (PDL_ABS(nu) > small) {
473			/* Computing 2nd power */
474			radcal = (nu - (2. * rho + 1.)) * nu + sigma * sigma;
475			if (radcal < 0.) {
476			(ierr_datap)[0] = -2;
477			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "NEGATIVE RADICAL");
478			}
479			that = (cp - sqrt(radcal)) / (3. * nu);
480			} else {
481			that = 1. / (2. * sigma);
482			}
483			}
484			phi = that * ((nu * that - cp) * that + 1.);
485			/* CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . */
486			if (indx != 3) {
487			phi -= rho;
488			}
489			/* TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. */
490			hphi = h2 * PDL_ABS(phi);
491			if (hphi * PDL_ABS(d1) > dfmx) {
492			/* AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK. */
493			d1 = ((d1) >= 0 ? PDL_ABS(dfmx / hphi) : -PDL_ABS(dfmx / hphi));
494			d2 = -lambda * d1;
495			}
496			}
497			(ierr_datap)[0] = 0;
498			} while (0); /* end inline dpchsw */
499			/* --------------------------------------------------------------- */
500			if ((ierr_datap)[0] != 0) {
501			break;
502			}
503			}} /* Close nless1=1 / / ....... END OF SEGMENT LOOP. */
504			#line 4884 "lib/PDL/Primitive.pd"
505			} while (0); /* end inline dpchcs */
506			}
507			}
508			/* SET END CONDITIONS. */
509			if (ibeg == 0 && iend == 0)
510			continue;
511			/* ------------------------------------------------------- */
512			do { /* inline dpchce */
513			/* Local variables */
514			PDL_Indx j, k, ibeg = (ic_datap)[0+(__inc_ic_two(0))], iend = (ic_datap)[0+(__inc_ic_two(1))];
515			PDL_Float stemp[3], xtemp[4];
516			/* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */
517			if (PDL_ABS(ibeg) > n)
518			ibeg = 0;
519			if (PDL_ABS(iend) > n)
520			iend = 0;
521			/* TREAT BEGINNING BOUNDARY CONDITION. */
522			if (ibeg != 0) {
523			k = PDL_ABS(ibeg);
524			if (k == 1) {
525			/* BOUNDARY VALUE PROVIDED. */
526			(d_datap)[0+(__inc_d_n(0))] = (vc_datap)[0+(__inc_vc_two(0))];
527			} else if (k == 2) {
528			/* BOUNDARY SECOND DERIVATIVE PROVIDED. */
529			(d_datap)[0+(__inc_d_n(0))] = 0.5 (3. * (slope_datap)[0+(__inc_slope_nless1(0))] - (d_datap)[0+(__inc_d_n(1))] - 0.5 * (vc_datap)[0+(__inc_vc_two(0))] (h_datap)[0+(__inc_h_nless1*(0))]);
530			} else if (k < 5) {
531			/* USE K-POINT DERIVATIVE FORMULA. */
532			/* PICK UP FIRST K POINTS, IN REVERSE ORDER. */
533			for (j = 0; j < k; ++j) {
534			PDL_Indx index = k - j;
535			/* INDEX RUNS FROM K DOWN TO 1. */
536			xtemp[j] = (x_datap)[0+(__inc_x_n*(index+1))];
537			if (j < k-1) {
538			stemp[j] = (slope_datap)[0+(__inc_slope_nless1*(index))];
539			}
540			}
541			/* ----------------------------- */
542
543			/* PDL version: K, X, S are var names, 4th param output */
544			/* **PURPOSE Computes divided differences for DPCHCE and DPCHSP /
545			/* DPCHDF: DPCHIP Finite Difference Formula */
546			/* Uses a divided difference formulation to compute a K-point approx- */
547			/* imation to the derivative at X(K) based on the data in X and S. */
548			/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
549			/* derivative approximations. */
550			/* ---------------------------------------------------------------------- */
551			/* On input: */
552			/* K is the order of the desired derivative approximation. */
553			/* K must be at least 3 (error return if not). */
554			/* X contains the K values of the independent variable. */
555			/* X need not be ordered, but the values MUST be */
556			/* distinct. (Not checked here.) */
557			/* S contains the associated slope values: */
558			/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
559			/* (Note that S need only be of length K-1.) */
560			/* On return: */
561			/* S will be destroyed. */
562			/* IERR will be set to -1 if K.LT.2 . */
563			/* DPCHDF will be set to the desired derivative approximation if */
564			/* IERR=0 or to zero if IERR=-1. */
565			/* ---------------------------------------------------------------------- */
566			/* **SEE ALSO DPCHCE, DPCHSP /
567			/* **REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- /
568			/* Verlag, New York, 1978, pp. 10-16. */
569			/* CHECK FOR LEGAL VALUE OF K. */
570			{
571			/* Local variables */
572			PDL_Indx i, j, k_cached = k;
573			PDL_Float x = xtemp, s = stemp;
574			if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "K LESS THAN THREE");
575			/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
576			for (j = 2; j < k_cached; ++j) {
577			PDL_Indx itmp = k_cached - j;
578			for (i = 0; i < itmp; ++i)
579			s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
580			}
581			/* EVALUATE DERIVATIVE AT X(K). */
582			PDL_Float value = s[0];
583			for (i = 1; i < k_cached-1; ++i)
584			value = s[i] + value * (x[k_cached-1] - x[i]);
585			(d_datap)[0+(__inc_d_n*(0))] = value;
586			}
587			;
588			/* ----------------------------- */
589			} else {
590			/* USE 'NOT A KNOT' CONDITION. */
591			(d_datap)[0+(__inc_d_n(0))] = (3. ((h_datap)[0+(__inc_h_nless1(0))] (slope_datap)[0+(__inc_slope_nless1(1))] + (h_datap)[0+(__inc_h_nless1(1))] * (slope_datap)[0+(__inc_slope_nless1*(0))]) -
592			2. * ((h_datap)[0+(__inc_h_nless1(0))] + (h_datap)[0+(__inc_h_nless1(1))]) * (d_datap)[0+(__inc_d_n(1))] - (h_datap)[0+(__inc_h_nless1(0))] * (d_datap)[0+(__inc_d_n(2))]) / (h_datap)[0+(__inc_h_nless1(1))];
593			}
594			/* CHECK D(1,1) FOR COMPATIBILITY WITH MONOTONICITY. */
595			if (ibeg <= 0) {
596			if ((slope_datap)[0+(__inc_slope_nless1*(0))] == 0.) {
597			if ((d_datap)[0+(__inc_d_n*(0))] != 0.) {
598			(d_datap)[0+(__inc_d_n*(0))] = 0.;
599			++((ierr_datap)[0]);
600			}
601			} else if (((((d_datap)[0+(__inc_d_n(0))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(0))]) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
602			(d_datap)[0+(__inc_d_n*(0))] = 0.;
603			++((ierr_datap)[0]);
604			} else if (PDL_ABS((d_datap)[0+(__inc_d_n(0))]) > 3. PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(0))])) {
605			(d_datap)[0+(__inc_d_n(0))] = 3. (slope_datap)[0+(__inc_slope_nless1*(0))];
606			++((ierr_datap)[0]);
607			}
608			}
609			}
610			/* TREAT END BOUNDARY CONDITION. */
611			if (iend == 0)
612			break;
613			k = PDL_ABS(iend);
614			if (k == 1) {
615			/* BOUNDARY VALUE PROVIDED. */
616			(d_datap)[0+(__inc_d_n(n-1))] = (vc_datap)[0+(__inc_vc_two(1))];
617			} else if (k == 2) {
618			/* BOUNDARY SECOND DERIVATIVE PROVIDED. */
619			(d_datap)[0+(__inc_d_n(n-1))] = 0.5 (3. * (slope_datap)[0+(__inc_slope_nless1(n-2))] - (d_datap)[0+(__inc_d_n(n-2))]
620			+ 0.5 * (vc_datap)[0+(__inc_vc_two(1))] (h_datap)[0+(__inc_h_nless1*(n-2))]);
621			} else if (k < 5) {
622			/* USE K-POINT DERIVATIVE FORMULA. */
623			/* PICK UP LAST K POINTS. */
624			for (j = 0; j < k; ++j) {
625			PDL_Indx index = n - k + j;
626			/* INDEX RUNS FROM N+1-K UP TO N. */
627			xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
628			if (j < k-1) {
629			stemp[j] = (slope_datap)[0+(__inc_slope_nless1*(index))];
630			}
631			}
632			/* ----------------------------- */
633
634			/* PDL version: K, X, S are var names, 4th param output */
635			/* **PURPOSE Computes divided differences for DPCHCE and DPCHSP /
636			/* DPCHDF: DPCHIP Finite Difference Formula */
637			/* Uses a divided difference formulation to compute a K-point approx- */
638			/* imation to the derivative at X(K) based on the data in X and S. */
639			/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
640			/* derivative approximations. */
641			/* ---------------------------------------------------------------------- */
642			/* On input: */
643			/* K is the order of the desired derivative approximation. */
644			/* K must be at least 3 (error return if not). */
645			/* X contains the K values of the independent variable. */
646			/* X need not be ordered, but the values MUST be */
647			/* distinct. (Not checked here.) */
648			/* S contains the associated slope values: */
649			/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
650			/* (Note that S need only be of length K-1.) */
651			/* On return: */
652			/* S will be destroyed. */
653			/* IERR will be set to -1 if K.LT.2 . */
654			/* DPCHDF will be set to the desired derivative approximation if */
655			/* IERR=0 or to zero if IERR=-1. */
656			/* ---------------------------------------------------------------------- */
657			/* **SEE ALSO DPCHCE, DPCHSP /
658			/* **REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- /
659			/* Verlag, New York, 1978, pp. 10-16. */
660			/* CHECK FOR LEGAL VALUE OF K. */
661			{
662			/* Local variables */
663			PDL_Indx i, j, k_cached = k;
664			PDL_Float x = xtemp, s = stemp;
665			if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "K LESS THAN THREE");
666			/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
667			for (j = 2; j < k_cached; ++j) {
668			PDL_Indx itmp = k_cached - j;
669			for (i = 0; i < itmp; ++i)
670			s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
671			}
672			/* EVALUATE DERIVATIVE AT X(K). */
673			PDL_Float value = s[0];
674			for (i = 1; i < k_cached-1; ++i)
675			value = s[i] + value * (x[k_cached-1] - x[i]);
676			(d_datap)[0+(__inc_d_n*(n-1))] = value;
677			}
678			;
679			/* ----------------------------- */
680			} else {
681			/* USE 'NOT A KNOT' CONDITION. */
682			(d_datap)[0+(__inc_d_n(n-1))] = (3. ((h_datap)[0+(__inc_h_nless1(n-2))] (slope_datap)[0+(__inc_slope_nless1*(n-3))] +
683			(h_datap)[0+(__inc_h_nless1(n-3))] (slope_datap)[0+(__inc_slope_nless1(n-2))]) - 2. ((h_datap)[0+(__inc_h_nless1(n-2))] + (h_datap)[0+(__inc_h_nless1(n-3))]) *
684			(d_datap)[0+(__inc_d_n(n-2))] - (h_datap)[0+(__inc_h_nless1(n-2))] * (d_datap)[0+(__inc_d_n(n-3))]) / (h_datap)[0+(__inc_h_nless1(n-3))];
685			}
686			if (iend > 0)
687			break;
688			/* CHECK D(1,N) FOR COMPATIBILITY WITH MONOTONICITY. */
689			if ((slope_datap)[0+(__inc_slope_nless1*(n-2))] == 0.) {
690			if ((d_datap)[0+(__inc_d_n*(n-1))] != 0.) {
691			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
692			(ierr_datap)[0] += 2;
693			}
694			} else if (((((d_datap)[0+(__inc_d_n(n-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(n-2))]) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(n-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(n-2))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
695			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
696			(ierr_datap)[0] += 2;
697			} else if (PDL_ABS((d_datap)[0+(__inc_d_n(n-1))]) > 3. PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(n-2))])) {
698			(d_datap)[0+(__inc_d_n(n-1))] = 3. (slope_datap)[0+(__inc_slope_nless1*(n-2))];
699			(ierr_datap)[0] += 2;
700			}
701			} while (0); /* end inlined dpchce */
702			/* ------------------------------------------------------- */
703			#line 704 "lib/PDL/Primitive-pp-pchip_chic.c"
704	0	0	}PDL_BROADCASTLOOP_END_pchip_chic_readdata
		0
705	0		} break;
706	3		case PDL_D: {
707	3	50	PDL_DECLARE_PARAMS_pchip_chic_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
		50
		50
		50
708	14	50	PDL_BROADCASTLOOP_START_pchip_chic_readdata {
		50
		50
		50
		50
		100
		100
709			#line 4574 "lib/PDL/Primitive.pd"
710			const PDL_Double d1mach = DBL_EPSILON;
711			/* VALIDITY-CHECK ARGUMENTS. */
712			{/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
713			#line 4577 "lib/PDL/Primitive.pd"
714			if ((x_datap)[0+(__inc_x_n(n))] > (x_datap)[0+(__inc_x_n(n-1))]) continue;
715			(ierr_datap)[0] = -1;
716			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "X-ARRAY NOT STRICTLY INCREASING");
717			}} /* Close n=1 */
718			#line 4581 "lib/PDL/Primitive.pd"
719			PDL_Indx ibeg = (ic_datap)[0+(__inc_ic_two(0))], iend = (ic_datap)[0+(__inc_ic_two(1))], n = __privtrans->ind_sizes[0];
720			(ierr_datap)[0] = 0;
721			if (PDL_ABS(ibeg) > 5)
722			--((ierr_datap)[0]);
723			if (PDL_ABS(iend) > 5)
724			(ierr_datap)[0] += -2;
725			if ((ierr_datap)[0] < 0) {
726			(ierr_datap)[0] += -3;
727			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "IC OUT OF RANGE");
728			}
729			/* FUNCTION DEFINITION IS OK -- GO ON. */
730			/* SET UP H AND SLOPE ARRAYS. */
731			{/* Open nless1 */ PDL_EXPAND2(register PDL_Indx nless1=0, __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
732			#line 4594 "lib/PDL/Primitive.pd"
733			(h_datap)[0+(__inc_h_nless1(nless1))] = (x_datap)[0+(__inc_x_n(nless1+1))] - (x_datap)[0+(__inc_x_n*(nless1))];
734			(slope_datap)[0+(__inc_slope_nless1(nless1))] = (f_datap)[0+(__inc_f_n(nless1+1))] - (f_datap)[0+(__inc_f_n*(nless1))];
735			}} /* Close nless1 */
736			#line 4597 "lib/PDL/Primitive.pd"
737			/* SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. */
738			if (__privtrans->ind_sizes[1] <= 1) {
739			(d_datap)[0+(__inc_d_n(0))] = (d_datap)[0+(__inc_d_n(1))] = (slope_datap)[0+(__inc_slope_nless1*(0))];
740			} else {
741			/* NORMAL CASE (N .GE. 3) . */
742			/* SET INTERIOR DERIVATIVES AND DEFAULT END CONDITIONS. */
743			do { /* inline dpchci */
744			/* Local variables */
745			PDL_Double del1 = (slope_datap)[0+(__inc_slope_nless1*(0))];
746			/* SPECIAL CASE N=2 is dealt with in separate branch above */
747			/* NORMAL CASE (N .GE. 3). */
748			PDL_Double del2 = (slope_datap)[0+(__inc_slope_nless1*(1))];
749			/* SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */
750			/* SHAPE-PRESERVING. */
751			PDL_Double hsum = (h_datap)[0+(__inc_h_nless1(0))] + (h_datap)[0+(__inc_h_nless1(1))];
752			PDL_Double w1 = ((h_datap)[0+(__inc_h_nless1*(0))] + hsum) / hsum;
753			PDL_Double w2 = -(h_datap)[0+(__inc_h_nless1*(0))] / hsum;
754			(d_datap)[0+(__inc_d_n(0))] = w1 del1 + w2 * del2;
755			if (((((d_datap)[0+(__inc_d_n(0))]) == 0. \|\| (del1) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.) {
756			(d_datap)[0+(__inc_d_n*(0))] = 0.;
757			} else if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
758			/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */
759			PDL_Double dmax = 3. * del1;
760			if (PDL_ABS((d_datap)[0+(__inc_d_n*(0))]) > PDL_ABS(dmax))
761			(d_datap)[0+(__inc_d_n*(0))] = dmax;
762			}
763			/* LOOP THROUGH INTERIOR POINTS. */
764			{/* Open nless1=1 */ PDL_EXPAND2(register PDL_Indx nless1=PDLMAX((1),0), __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
765			#line 4625 "lib/PDL/Primitive.pd"
766			if (nless1 != 1) {
767			hsum = (h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))];
768			del1 = del2;
769			del2 = (slope_datap)[0+(__inc_slope_nless1*(nless1))];
770			}
771			/* SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. */
772			(d_datap)[0+(__inc_d_n*(nless1))] = 0.;
773			if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.)
774			continue;
775			/* USE BRODLIE MODIFICATION OF BUTLAND FORMULA. */
776			PDL_Double hsumt3 = hsum + hsum + hsum;
777			w1 = (hsum + (h_datap)[0+(__inc_h_nless1*(nless1-1))]) / hsumt3;
778			w2 = (hsum + (h_datap)[0+(__inc_h_nless1*(nless1))]) / hsumt3;
779			/* Computing MAX */
780			PDL_Double dmax = PDLMAX(PDL_ABS(del1),PDL_ABS(del2));
781			/* Computing MIN */
782			PDL_Double dmin = PDLMIN(PDL_ABS(del1),PDL_ABS(del2));
783			PDL_Double drat1 = del1 / dmax, drat2 = del2 / dmax;
784			(d_datap)[0+(__inc_d_n(nless1))] = dmin / (w1 drat1 + w2 * drat2);
785			}} /* Close nless1=1 */
786			#line 4645 "lib/PDL/Primitive.pd"
787			/* SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */
788			/* SHAPE-PRESERVING. */
789			w1 = -(h_datap)[0+(__inc_h_nless1*(n-2))] / hsum;
790			w2 = ((h_datap)[0+(__inc_h_nless1*(n-2))] + hsum) / hsum;
791			(d_datap)[0+(__inc_d_n(n-1))] = w1 del1 + w2 * del2;
792			if (((((d_datap)[0+(__inc_d_n(n-1))]) == 0. \|\| (del2) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(n-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.) {
793			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
794			} else if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
795			/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */
796			PDL_Double dmax = 3. * del2;
797			if (PDL_ABS((d_datap)[0+(__inc_d_n*(n-1))]) > PDL_ABS(dmax))
798			(d_datap)[0+(__inc_d_n*(n-1))] = dmax;
799			}
800			} while (0); /* end inline dpchci */
801			/* SET DERIVATIVES AT POINTS WHERE MONOTONICITY SWITCHES DIRECTION. */
802			if ((mflag_datap)[0] != 0.) {
803			do { /* inline dpchcs */
804			/* **PURPOSE Adjusts derivative values for DPCHIC /
805			/* DPCHCS: DPCHIC Monotonicity Switch Derivative Setter. */
806			/* Called by DPCHIC to adjust the values of D in the vicinity of a */
807			/* switch in direction of monotonicity, to produce a more "visually */
808			/* pleasing" curve than that given by DPCHIM . */
809			static const PDL_Double fudge = 4.;
810			/* Local variables */
811			PDL_Indx k;
812			PDL_Double del[3], fact, dfmx;
813			PDL_Double dext, dfloc, slmax, wtave[2];
814			/* INITIALIZE. */
815			/* LOOP OVER SEGMENTS. */
816			{/* Open nless1=1 */ PDL_EXPAND2(register PDL_Indx nless1=PDLMAX((1),0), __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
817			#line 4675 "lib/PDL/Primitive.pd"
818			PDL_Double dtmp = ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)));
819			if (dtmp > 0.) {
820			continue;
821			}
822			if (dtmp != 0.) {
823			/* ....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT ..................... */
824			/* DO NOT CHANGE D IF 'UP-DOWN-UP'. */
825			if (nless1 > 1) {
826			if (((((slope_datap)[0+(__inc_slope_nless1(nless1-2))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-2))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) > 0.)
827			continue;
828			/* -------------------------- */
829			}
830			if (nless1 < __privtrans->ind_sizes[1]-1 && ((((slope_datap)[0+(__inc_slope_nless1(nless1+1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1+1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) > 0.)
831			continue;
832			/* ---------------------------- */
833			/* ....... COMPUTE PROVISIONAL VALUE FOR D(1,I). */
834			dext = (h_datap)[0+(__inc_h_nless1(nless1))] / ((h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))]) (slope_datap)[0+(__inc_slope_nless1*(nless1-1))] +
835			(h_datap)[0+(__inc_h_nless1(nless1-1))] / ((h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))]) (slope_datap)[0+(__inc_slope_nless1*(nless1))];
836			/* ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM. */
837			dtmp = (((dext) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0.) ? 0. : (((dext)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)));
838			if (dtmp == 0) {
839			continue;
840			}
841			if (dtmp < 0.) {
842			/* DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS -- */
843			/* EXTREMUM IS IN (X(I-1),X(I)). */
844			k = nless1;
845			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I). */
846			wtave[1] = dext;
847			if (k > 1) {
848			wtave[0] = (h_datap)[0+(__inc_h_nless1(k-1))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) (slope_datap)[0+(__inc_slope_nless1*(k-2))] +
849			(h_datap)[0+(__inc_h_nless1(k-2))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1*(k-1))];
850			}
851			} else {
852			/* DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS -- */
853			/* EXTREMUM IS IN (X(I),X(I+1)). */
854			k = nless1 + 1;
855			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). */
856			wtave[0] = dext;
857			if (k < nless1) {
858			wtave[1] = (h_datap)[0+(__inc_h_nless1(k))] / ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k-1))]
859			/ ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) * (slope_datap)[0+(__inc_slope_nless1*(k))];
860			}
861			}
862			} else {
863			/* ....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO -- */
864			/* CHECK FOR FLAT-TOPPED PEAK ....................... */
865			if (nless1 == __privtrans->ind_sizes[1]-1 \|\| ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1+1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1+1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) >= 0.)
866			continue;
867			/* ----------------------------- */
868			/* WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)). */
869			k = nless1+1;
870			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). */
871			wtave[0] = (h_datap)[0+(__inc_h_nless1(k-1))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) (slope_datap)[0+(__inc_slope_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-2))]
872			/ ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) * (slope_datap)[0+(__inc_slope_nless1*(k-1))];
873			wtave[1] = (h_datap)[0+(__inc_h_nless1(k))] / ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k-1))] / (
874			(h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) * (slope_datap)[0+(__inc_slope_nless1*(k))];
875			}
876			/* ....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM */
877			/* ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE-- */
878			/* WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K), */
879			/* IF K.GT.1 */
880			/* WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1), */
881			/* IF K.LT.N-1 */
882			slmax = PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-1))]);
883			if (k > 1) {
884			/* Computing MAX */
885			slmax = PDLMAX(slmax,PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-2))]));
886			}
887			if (k < nless1) {
888			/* Computing MAX */
889			slmax = PDLMAX(slmax,PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k))]));
890			}
891			if (k > 1) {
892			del[0] = (slope_datap)[0+(__inc_slope_nless1*(k-2))] / slmax;
893			}
894			del[1] = (slope_datap)[0+(__inc_slope_nless1*(k-1))] / slmax;
895			if (k < nless1) {
896			del[2] = (slope_datap)[0+(__inc_slope_nless1*(k))] / slmax;
897			}
898			if (k > 1 && k < nless1) {
899			/* NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL. */
900			fact = fudge * PDL_ABS(del[2] * (del[0] - del[1]) * (wtave[1] / slmax));
901			(d_datap)[0+(__inc_d_n(k-1))] += PDLMIN(fact,1.) (wtave[0] - (d_datap)[0+(__inc_d_n*(k-1))]);
902			fact = fudge * PDL_ABS(del[0] * (del[2] - del[1]) * (wtave[0] / slmax));
903			(d_datap)[0+(__inc_d_n(k))] += PDLMIN(fact,1.) (wtave[1] - (d_datap)[0+(__inc_d_n*(k))]);
904			} else {
905			/* SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR */
906			/* K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1). */
907			fact = fudge * PDL_ABS(del[1]);
908			(d_datap)[0+(__inc_d_n(nless1))] = PDLMIN(fact,1.) wtave[nless1+1 - k];
909			/* NOTE THAT I-K+1 = 1 IF K=I (=NLESS1), */
910			/* I-K+1 = 2 IF K=I-1(=1). */
911			}
912			/* ....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA. */
913			if ((mflag_datap)[0] <= 0.) {
914			continue;
915			}
916			dfloc = (h_datap)[0+(__inc_h_nless1(k-1))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-1))]);
917			if (k > 1) {
918			/* Computing MAX */
919			dfloc = PDLMAX(dfloc,(h_datap)[0+(__inc_h_nless1(k-2))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-2))]));
920			}
921			if (k < nless1) {
922			/* Computing MAX */
923			dfloc = PDLMAX(dfloc,(h_datap)[0+(__inc_h_nless1(k))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k))]));
924			}
925			dfmx = (mflag_datap)[0] * dfloc;
926			PDL_Indx indx = nless1 - k;
927			/* INDX = 1 IF K=I, 2 IF K=I-1. */
928			/* --------------------------------------------------------------- */
929			do { /* inline dpchsw */
930			/* NOTATION AND GENERAL REMARKS. */
931			/* RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED. */
932			/* LAMBDA IS THE RATIO OF D2 TO D1. */
933			/* THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM. */
934			/* PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO), */
935			/* WHERE THAT = (XHAT - X1)/H . */
936			/* THAT IS, P(XHAT)-F1 = DHPHI, WHERE D=D1 OR D2. */
937			/* SIMILARLY, P(XHAT)-F2 = DH(PHI-RHO) . */
938			/* Local variables */
939			PDL_Double cp, nu, phi, rho, hphi, that, sigma, small;
940			PDL_Double lambda, radcal;
941			PDL_Double d1 = (d_datap)[0+(__inc_d_n(k-1))], d2 = (d_datap)[0+(__inc_d_n(k))], h2 = (h_datap)[0+(__inc_h_nless1(k-1))], slope2 = (slope_datap)[0+(__inc_slope_nless1(k-1))];
942			/* Initialized data */
943			static const PDL_Double fact = 100.;
944			/* THIRD SHOULD BE SLIGHTLY LESS THAN 1/3. */
945			static const PDL_Double third = .33333;
946			/* SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS. */
947			small = fact * d1mach;
948			/* DO MAIN CALCULATION. */
949			if (d1 == 0.) {
950			/* SPECIAL CASE -- D1.EQ.ZERO . */
951			/* IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED. */
952			if (d2 == 0.) {
953			(ierr_datap)[0] = -1;
954			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "D1 AND/OR D2 INVALID");
955			}
956			rho = slope2 / d2;
957			/* EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . */
958			if (rho >= third) {
959			(ierr_datap)[0] = 0; break;
960			}
961			that = 2. * (3. * rho - 1.) / (3. * (2. * rho - 1.));
962			/* Computing 2nd power */
963			phi = that * that * ((3. * rho - 1.) / 3.);
964			/* CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . */
965			if (indx != 3) {
966			phi -= rho;
967			}
968			/* TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. */
969			hphi = h2 * PDL_ABS(phi);
970			if (hphi * PDL_ABS(d2) > dfmx) {
971			/* AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK. */
972			d2 = ((d2) >= 0 ? PDL_ABS(dfmx / hphi) : -PDL_ABS(dfmx / hphi));
973			}
974			} else {
975			rho = slope2 / d1;
976			lambda = -(d2) / d1;
977			if (d2 == 0.) {
978			/* SPECIAL CASE -- D2.EQ.ZERO . */
979			/* EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . */
980			if (rho >= third) {
981			(ierr_datap)[0] = 0; break;
982			}
983			cp = 2. - 3. * rho;
984			nu = 1. - 2. * rho;
985			that = 1. / (3. * nu);
986			} else {
987			if (lambda <= 0.) {
988			(ierr_datap)[0] = -1;
989			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "D1 AND/OR D2 INVALID");
990			}
991			/* NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS. */
992			nu = 1. - lambda - 2. * rho;
993			sigma = 1. - rho;
994			cp = nu + sigma;
995			if (PDL_ABS(nu) > small) {
996			/* Computing 2nd power */
997			radcal = (nu - (2. * rho + 1.)) * nu + sigma * sigma;
998			if (radcal < 0.) {
999			(ierr_datap)[0] = -2;
1000			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "NEGATIVE RADICAL");
1001			}
1002			that = (cp - sqrt(radcal)) / (3. * nu);
1003			} else {
1004			that = 1. / (2. * sigma);
1005			}
1006			}
1007			phi = that * ((nu * that - cp) * that + 1.);
1008			/* CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . */
1009			if (indx != 3) {
1010			phi -= rho;
1011			}
1012			/* TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. */
1013			hphi = h2 * PDL_ABS(phi);
1014			if (hphi * PDL_ABS(d1) > dfmx) {
1015			/* AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK. */
1016			d1 = ((d1) >= 0 ? PDL_ABS(dfmx / hphi) : -PDL_ABS(dfmx / hphi));
1017			d2 = -lambda * d1;
1018			}
1019			}
1020			(ierr_datap)[0] = 0;
1021			} while (0); /* end inline dpchsw */
1022			/* --------------------------------------------------------------- */
1023			if ((ierr_datap)[0] != 0) {
1024			break;
1025			}
1026			}} /* Close nless1=1 / / ....... END OF SEGMENT LOOP. */
1027			#line 4884 "lib/PDL/Primitive.pd"
1028			} while (0); /* end inline dpchcs */
1029			}
1030			}
1031			/* SET END CONDITIONS. */
1032			if (ibeg == 0 && iend == 0)
1033			continue;
1034			/* ------------------------------------------------------- */
1035			do { /* inline dpchce */
1036			/* Local variables */
1037			PDL_Indx j, k, ibeg = (ic_datap)[0+(__inc_ic_two(0))], iend = (ic_datap)[0+(__inc_ic_two(1))];
1038			PDL_Double stemp[3], xtemp[4];
1039			/* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */
1040			if (PDL_ABS(ibeg) > n)
1041			ibeg = 0;
1042			if (PDL_ABS(iend) > n)
1043			iend = 0;
1044			/* TREAT BEGINNING BOUNDARY CONDITION. */
1045			if (ibeg != 0) {
1046			k = PDL_ABS(ibeg);
1047			if (k == 1) {
1048			/* BOUNDARY VALUE PROVIDED. */
1049			(d_datap)[0+(__inc_d_n(0))] = (vc_datap)[0+(__inc_vc_two(0))];
1050			} else if (k == 2) {
1051			/* BOUNDARY SECOND DERIVATIVE PROVIDED. */
1052			(d_datap)[0+(__inc_d_n(0))] = 0.5 (3. * (slope_datap)[0+(__inc_slope_nless1(0))] - (d_datap)[0+(__inc_d_n(1))] - 0.5 * (vc_datap)[0+(__inc_vc_two(0))] (h_datap)[0+(__inc_h_nless1*(0))]);
1053			} else if (k < 5) {
1054			/* USE K-POINT DERIVATIVE FORMULA. */
1055			/* PICK UP FIRST K POINTS, IN REVERSE ORDER. */
1056			for (j = 0; j < k; ++j) {
1057			PDL_Indx index = k - j;
1058			/* INDEX RUNS FROM K DOWN TO 1. */
1059			xtemp[j] = (x_datap)[0+(__inc_x_n*(index+1))];
1060			if (j < k-1) {
1061			stemp[j] = (slope_datap)[0+(__inc_slope_nless1*(index))];
1062			}
1063			}
1064			/* ----------------------------- */
1065
1066			/* PDL version: K, X, S are var names, 4th param output */
1067			/* **PURPOSE Computes divided differences for DPCHCE and DPCHSP /
1068			/* DPCHDF: DPCHIP Finite Difference Formula */
1069			/* Uses a divided difference formulation to compute a K-point approx- */
1070			/* imation to the derivative at X(K) based on the data in X and S. */
1071			/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
1072			/* derivative approximations. */
1073			/* ---------------------------------------------------------------------- */
1074			/* On input: */
1075			/* K is the order of the desired derivative approximation. */
1076			/* K must be at least 3 (error return if not). */
1077			/* X contains the K values of the independent variable. */
1078			/* X need not be ordered, but the values MUST be */
1079			/* distinct. (Not checked here.) */
1080			/* S contains the associated slope values: */
1081			/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
1082			/* (Note that S need only be of length K-1.) */
1083			/* On return: */
1084			/* S will be destroyed. */
1085			/* IERR will be set to -1 if K.LT.2 . */
1086			/* DPCHDF will be set to the desired derivative approximation if */
1087			/* IERR=0 or to zero if IERR=-1. */
1088			/* ---------------------------------------------------------------------- */
1089			/* **SEE ALSO DPCHCE, DPCHSP /
1090			/* **REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- /
1091			/* Verlag, New York, 1978, pp. 10-16. */
1092			/* CHECK FOR LEGAL VALUE OF K. */
1093			{
1094			/* Local variables */
1095			PDL_Indx i, j, k_cached = k;
1096			PDL_Double x = xtemp, s = stemp;
1097			if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "K LESS THAN THREE");
1098			/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
1099			for (j = 2; j < k_cached; ++j) {
1100			PDL_Indx itmp = k_cached - j;
1101			for (i = 0; i < itmp; ++i)
1102			s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
1103			}
1104			/* EVALUATE DERIVATIVE AT X(K). */
1105			PDL_Double value = s[0];
1106			for (i = 1; i < k_cached-1; ++i)
1107			value = s[i] + value * (x[k_cached-1] - x[i]);
1108			(d_datap)[0+(__inc_d_n*(0))] = value;
1109			}
1110			;
1111			/* ----------------------------- */
1112			} else {
1113			/* USE 'NOT A KNOT' CONDITION. */
1114			(d_datap)[0+(__inc_d_n(0))] = (3. ((h_datap)[0+(__inc_h_nless1(0))] (slope_datap)[0+(__inc_slope_nless1(1))] + (h_datap)[0+(__inc_h_nless1(1))] * (slope_datap)[0+(__inc_slope_nless1*(0))]) -
1115			2. * ((h_datap)[0+(__inc_h_nless1(0))] + (h_datap)[0+(__inc_h_nless1(1))]) * (d_datap)[0+(__inc_d_n(1))] - (h_datap)[0+(__inc_h_nless1(0))] * (d_datap)[0+(__inc_d_n(2))]) / (h_datap)[0+(__inc_h_nless1(1))];
1116			}
1117			/* CHECK D(1,1) FOR COMPATIBILITY WITH MONOTONICITY. */
1118			if (ibeg <= 0) {
1119			if ((slope_datap)[0+(__inc_slope_nless1*(0))] == 0.) {
1120			if ((d_datap)[0+(__inc_d_n*(0))] != 0.) {
1121			(d_datap)[0+(__inc_d_n*(0))] = 0.;
1122			++((ierr_datap)[0]);
1123			}
1124			} else if (((((d_datap)[0+(__inc_d_n(0))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(0))]) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
1125			(d_datap)[0+(__inc_d_n*(0))] = 0.;
1126			++((ierr_datap)[0]);
1127			} else if (PDL_ABS((d_datap)[0+(__inc_d_n(0))]) > 3. PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(0))])) {
1128			(d_datap)[0+(__inc_d_n(0))] = 3. (slope_datap)[0+(__inc_slope_nless1*(0))];
1129			++((ierr_datap)[0]);
1130			}
1131			}
1132			}
1133			/* TREAT END BOUNDARY CONDITION. */
1134			if (iend == 0)
1135			break;
1136			k = PDL_ABS(iend);
1137			if (k == 1) {
1138			/* BOUNDARY VALUE PROVIDED. */
1139			(d_datap)[0+(__inc_d_n(n-1))] = (vc_datap)[0+(__inc_vc_two(1))];
1140			} else if (k == 2) {
1141			/* BOUNDARY SECOND DERIVATIVE PROVIDED. */
1142			(d_datap)[0+(__inc_d_n(n-1))] = 0.5 (3. * (slope_datap)[0+(__inc_slope_nless1(n-2))] - (d_datap)[0+(__inc_d_n(n-2))]
1143			+ 0.5 * (vc_datap)[0+(__inc_vc_two(1))] (h_datap)[0+(__inc_h_nless1*(n-2))]);
1144			} else if (k < 5) {
1145			/* USE K-POINT DERIVATIVE FORMULA. */
1146			/* PICK UP LAST K POINTS. */
1147			for (j = 0; j < k; ++j) {
1148			PDL_Indx index = n - k + j;
1149			/* INDEX RUNS FROM N+1-K UP TO N. */
1150			xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
1151			if (j < k-1) {
1152			stemp[j] = (slope_datap)[0+(__inc_slope_nless1*(index))];
1153			}
1154			}
1155			/* ----------------------------- */
1156
1157			/* PDL version: K, X, S are var names, 4th param output */
1158			/* **PURPOSE Computes divided differences for DPCHCE and DPCHSP /
1159			/* DPCHDF: DPCHIP Finite Difference Formula */
1160			/* Uses a divided difference formulation to compute a K-point approx- */
1161			/* imation to the derivative at X(K) based on the data in X and S. */
1162			/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
1163			/* derivative approximations. */
1164			/* ---------------------------------------------------------------------- */
1165			/* On input: */
1166			/* K is the order of the desired derivative approximation. */
1167			/* K must be at least 3 (error return if not). */
1168			/* X contains the K values of the independent variable. */
1169			/* X need not be ordered, but the values MUST be */
1170			/* distinct. (Not checked here.) */
1171			/* S contains the associated slope values: */
1172			/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
1173			/* (Note that S need only be of length K-1.) */
1174			/* On return: */
1175			/* S will be destroyed. */
1176			/* IERR will be set to -1 if K.LT.2 . */
1177			/* DPCHDF will be set to the desired derivative approximation if */
1178			/* IERR=0 or to zero if IERR=-1. */
1179			/* ---------------------------------------------------------------------- */
1180			/* **SEE ALSO DPCHCE, DPCHSP /
1181			/* **REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- /
1182			/* Verlag, New York, 1978, pp. 10-16. */
1183			/* CHECK FOR LEGAL VALUE OF K. */
1184			{
1185			/* Local variables */
1186			PDL_Indx i, j, k_cached = k;
1187			PDL_Double x = xtemp, s = stemp;
1188			if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "K LESS THAN THREE");
1189			/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
1190			for (j = 2; j < k_cached; ++j) {
1191			PDL_Indx itmp = k_cached - j;
1192			for (i = 0; i < itmp; ++i)
1193			s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
1194			}
1195			/* EVALUATE DERIVATIVE AT X(K). */
1196			PDL_Double value = s[0];
1197			for (i = 1; i < k_cached-1; ++i)
1198			value = s[i] + value * (x[k_cached-1] - x[i]);
1199			(d_datap)[0+(__inc_d_n*(n-1))] = value;
1200			}
1201			;
1202			/* ----------------------------- */
1203			} else {
1204			/* USE 'NOT A KNOT' CONDITION. */
1205			(d_datap)[0+(__inc_d_n(n-1))] = (3. ((h_datap)[0+(__inc_h_nless1(n-2))] (slope_datap)[0+(__inc_slope_nless1*(n-3))] +
1206			(h_datap)[0+(__inc_h_nless1(n-3))] (slope_datap)[0+(__inc_slope_nless1(n-2))]) - 2. ((h_datap)[0+(__inc_h_nless1(n-2))] + (h_datap)[0+(__inc_h_nless1(n-3))]) *
1207			(d_datap)[0+(__inc_d_n(n-2))] - (h_datap)[0+(__inc_h_nless1(n-2))] * (d_datap)[0+(__inc_d_n(n-3))]) / (h_datap)[0+(__inc_h_nless1(n-3))];
1208			}
1209			if (iend > 0)
1210			break;
1211			/* CHECK D(1,N) FOR COMPATIBILITY WITH MONOTONICITY. */
1212			if ((slope_datap)[0+(__inc_slope_nless1*(n-2))] == 0.) {
1213			if ((d_datap)[0+(__inc_d_n*(n-1))] != 0.) {
1214			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
1215			(ierr_datap)[0] += 2;
1216			}
1217			} else if (((((d_datap)[0+(__inc_d_n(n-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(n-2))]) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(n-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(n-2))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
1218			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
1219			(ierr_datap)[0] += 2;
1220			} else if (PDL_ABS((d_datap)[0+(__inc_d_n(n-1))]) > 3. PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(n-2))])) {
1221			(d_datap)[0+(__inc_d_n(n-1))] = 3. (slope_datap)[0+(__inc_slope_nless1*(n-2))];
1222			(ierr_datap)[0] += 2;
1223			}
1224			} while (0); /* end inlined dpchce */
1225			/* ------------------------------------------------------- */
1226			#line 1227 "lib/PDL/Primitive-pp-pchip_chic.c"
1227	3	50	}PDL_BROADCASTLOOP_END_pchip_chic_readdata
		50
1228	3		} break;
1229	0		case PDL_LD: {
1230	0	0	PDL_DECLARE_PARAMS_pchip_chic_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
		0
		0
		0
1231	0	0	PDL_BROADCASTLOOP_START_pchip_chic_readdata {
		0
		0
		0
		0
		0
		0
1232			#line 4574 "lib/PDL/Primitive.pd"
1233			const PDL_LDouble d1mach = LDBL_EPSILON;
1234			/* VALIDITY-CHECK ARGUMENTS. */
1235			{/* Open n=1 */ PDL_EXPAND2(register PDL_Indx n=PDLMAX((1),0), __n_stop=(__n_size)); for(; n<__n_stop; n+=1) {
1236			#line 4577 "lib/PDL/Primitive.pd"
1237			if ((x_datap)[0+(__inc_x_n(n))] > (x_datap)[0+(__inc_x_n(n-1))]) continue;
1238			(ierr_datap)[0] = -1;
1239			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "X-ARRAY NOT STRICTLY INCREASING");
1240			}} /* Close n=1 */
1241			#line 4581 "lib/PDL/Primitive.pd"
1242			PDL_Indx ibeg = (ic_datap)[0+(__inc_ic_two(0))], iend = (ic_datap)[0+(__inc_ic_two(1))], n = __privtrans->ind_sizes[0];
1243			(ierr_datap)[0] = 0;
1244			if (PDL_ABS(ibeg) > 5)
1245			--((ierr_datap)[0]);
1246			if (PDL_ABS(iend) > 5)
1247			(ierr_datap)[0] += -2;
1248			if ((ierr_datap)[0] < 0) {
1249			(ierr_datap)[0] += -3;
1250			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "IC OUT OF RANGE");
1251			}
1252			/* FUNCTION DEFINITION IS OK -- GO ON. */
1253			/* SET UP H AND SLOPE ARRAYS. */
1254			{/* Open nless1 */ PDL_EXPAND2(register PDL_Indx nless1=0, __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
1255			#line 4594 "lib/PDL/Primitive.pd"
1256			(h_datap)[0+(__inc_h_nless1(nless1))] = (x_datap)[0+(__inc_x_n(nless1+1))] - (x_datap)[0+(__inc_x_n*(nless1))];
1257			(slope_datap)[0+(__inc_slope_nless1(nless1))] = (f_datap)[0+(__inc_f_n(nless1+1))] - (f_datap)[0+(__inc_f_n*(nless1))];
1258			}} /* Close nless1 */
1259			#line 4597 "lib/PDL/Primitive.pd"
1260			/* SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION. */
1261			if (__privtrans->ind_sizes[1] <= 1) {
1262			(d_datap)[0+(__inc_d_n(0))] = (d_datap)[0+(__inc_d_n(1))] = (slope_datap)[0+(__inc_slope_nless1*(0))];
1263			} else {
1264			/* NORMAL CASE (N .GE. 3) . */
1265			/* SET INTERIOR DERIVATIVES AND DEFAULT END CONDITIONS. */
1266			do { /* inline dpchci */
1267			/* Local variables */
1268			PDL_LDouble del1 = (slope_datap)[0+(__inc_slope_nless1*(0))];
1269			/* SPECIAL CASE N=2 is dealt with in separate branch above */
1270			/* NORMAL CASE (N .GE. 3). */
1271			PDL_LDouble del2 = (slope_datap)[0+(__inc_slope_nless1*(1))];
1272			/* SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */
1273			/* SHAPE-PRESERVING. */
1274			PDL_LDouble hsum = (h_datap)[0+(__inc_h_nless1(0))] + (h_datap)[0+(__inc_h_nless1(1))];
1275			PDL_LDouble w1 = ((h_datap)[0+(__inc_h_nless1*(0))] + hsum) / hsum;
1276			PDL_LDouble w2 = -(h_datap)[0+(__inc_h_nless1*(0))] / hsum;
1277			(d_datap)[0+(__inc_d_n(0))] = w1 del1 + w2 * del2;
1278			if (((((d_datap)[0+(__inc_d_n(0))]) == 0. \|\| (del1) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.) {
1279			(d_datap)[0+(__inc_d_n*(0))] = 0.;
1280			} else if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
1281			/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */
1282			PDL_LDouble dmax = 3. * del1;
1283			if (PDL_ABS((d_datap)[0+(__inc_d_n*(0))]) > PDL_ABS(dmax))
1284			(d_datap)[0+(__inc_d_n*(0))] = dmax;
1285			}
1286			/* LOOP THROUGH INTERIOR POINTS. */
1287			{/* Open nless1=1 */ PDL_EXPAND2(register PDL_Indx nless1=PDLMAX((1),0), __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
1288			#line 4625 "lib/PDL/Primitive.pd"
1289			if (nless1 != 1) {
1290			hsum = (h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))];
1291			del1 = del2;
1292			del2 = (slope_datap)[0+(__inc_slope_nless1*(nless1))];
1293			}
1294			/* SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC. */
1295			(d_datap)[0+(__inc_d_n*(nless1))] = 0.;
1296			if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.)
1297			continue;
1298			/* USE BRODLIE MODIFICATION OF BUTLAND FORMULA. */
1299			PDL_LDouble hsumt3 = hsum + hsum + hsum;
1300			w1 = (hsum + (h_datap)[0+(__inc_h_nless1*(nless1-1))]) / hsumt3;
1301			w2 = (hsum + (h_datap)[0+(__inc_h_nless1*(nless1))]) / hsumt3;
1302			/* Computing MAX */
1303			PDL_LDouble dmax = PDLMAX(PDL_ABS(del1),PDL_ABS(del2));
1304			/* Computing MIN */
1305			PDL_LDouble dmin = PDLMIN(PDL_ABS(del1),PDL_ABS(del2));
1306			PDL_LDouble drat1 = del1 / dmax, drat2 = del2 / dmax;
1307			(d_datap)[0+(__inc_d_n(nless1))] = dmin / (w1 drat1 + w2 * drat2);
1308			}} /* Close nless1=1 */
1309			#line 4645 "lib/PDL/Primitive.pd"
1310			/* SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE */
1311			/* SHAPE-PRESERVING. */
1312			w1 = -(h_datap)[0+(__inc_h_nless1*(n-2))] / hsum;
1313			w2 = ((h_datap)[0+(__inc_h_nless1*(n-2))] + hsum) / hsum;
1314			(d_datap)[0+(__inc_d_n(n-1))] = w1 del1 + w2 * del2;
1315			if (((((d_datap)[0+(__inc_d_n(n-1))]) == 0. \|\| (del2) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(n-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) <= 0.) {
1316			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
1317			} else if ((((del1) == 0. \|\| (del2) == 0.) ? 0. : (((del1)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) * (((del2)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
1318			/* NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES. */
1319			PDL_LDouble dmax = 3. * del2;
1320			if (PDL_ABS((d_datap)[0+(__inc_d_n*(n-1))]) > PDL_ABS(dmax))
1321			(d_datap)[0+(__inc_d_n*(n-1))] = dmax;
1322			}
1323			} while (0); /* end inline dpchci */
1324			/* SET DERIVATIVES AT POINTS WHERE MONOTONICITY SWITCHES DIRECTION. */
1325			if ((mflag_datap)[0] != 0.) {
1326			do { /* inline dpchcs */
1327			/* **PURPOSE Adjusts derivative values for DPCHIC /
1328			/* DPCHCS: DPCHIC Monotonicity Switch Derivative Setter. */
1329			/* Called by DPCHIC to adjust the values of D in the vicinity of a */
1330			/* switch in direction of monotonicity, to produce a more "visually */
1331			/* pleasing" curve than that given by DPCHIM . */
1332			static const PDL_LDouble fudge = 4.;
1333			/* Local variables */
1334			PDL_Indx k;
1335			PDL_LDouble del[3], fact, dfmx;
1336			PDL_LDouble dext, dfloc, slmax, wtave[2];
1337			/* INITIALIZE. */
1338			/* LOOP OVER SEGMENTS. */
1339			{/* Open nless1=1 */ PDL_EXPAND2(register PDL_Indx nless1=PDLMAX((1),0), __nless1_stop=(__nless1_size)); for(; nless1<__nless1_stop; nless1+=1) {
1340			#line 4675 "lib/PDL/Primitive.pd"
1341			PDL_LDouble dtmp = ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)));
1342			if (dtmp > 0.) {
1343			continue;
1344			}
1345			if (dtmp != 0.) {
1346			/* ....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT ..................... */
1347			/* DO NOT CHANGE D IF 'UP-DOWN-UP'. */
1348			if (nless1 > 1) {
1349			if (((((slope_datap)[0+(__inc_slope_nless1(nless1-2))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-2))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) > 0.)
1350			continue;
1351			/* -------------------------- */
1352			}
1353			if (nless1 < __privtrans->ind_sizes[1]-1 && ((((slope_datap)[0+(__inc_slope_nless1(nless1+1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1+1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) > 0.)
1354			continue;
1355			/* ---------------------------- */
1356			/* ....... COMPUTE PROVISIONAL VALUE FOR D(1,I). */
1357			dext = (h_datap)[0+(__inc_h_nless1(nless1))] / ((h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))]) (slope_datap)[0+(__inc_slope_nless1*(nless1-1))] +
1358			(h_datap)[0+(__inc_h_nless1(nless1-1))] / ((h_datap)[0+(__inc_h_nless1(nless1-1))] + (h_datap)[0+(__inc_h_nless1(nless1))]) (slope_datap)[0+(__inc_slope_nless1*(nless1))];
1359			/* ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM. */
1360			dtmp = (((dext) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0.) ? 0. : (((dext)) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)));
1361			if (dtmp == 0) {
1362			continue;
1363			}
1364			if (dtmp < 0.) {
1365			/* DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS -- */
1366			/* EXTREMUM IS IN (X(I-1),X(I)). */
1367			k = nless1;
1368			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I). */
1369			wtave[1] = dext;
1370			if (k > 1) {
1371			wtave[0] = (h_datap)[0+(__inc_h_nless1(k-1))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) (slope_datap)[0+(__inc_slope_nless1*(k-2))] +
1372			(h_datap)[0+(__inc_h_nless1(k-2))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1*(k-1))];
1373			}
1374			} else {
1375			/* DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS -- */
1376			/* EXTREMUM IS IN (X(I),X(I+1)). */
1377			k = nless1 + 1;
1378			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). */
1379			wtave[0] = dext;
1380			if (k < nless1) {
1381			wtave[1] = (h_datap)[0+(__inc_h_nless1(k))] / ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k-1))]
1382			/ ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) * (slope_datap)[0+(__inc_slope_nless1*(k))];
1383			}
1384			}
1385			} else {
1386			/* ....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO -- */
1387			/* CHECK FOR FLAT-TOPPED PEAK ....................... */
1388			if (nless1 == __privtrans->ind_sizes[1]-1 \|\| ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(nless1+1))]) == 0.) ? 0. : ((((slope_datap)[0+(__inc_slope_nless1(nless1-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(nless1+1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) >= 0.)
1389			continue;
1390			/* ----------------------------- */
1391			/* WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)). */
1392			k = nless1+1;
1393			/* SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1). */
1394			wtave[0] = (h_datap)[0+(__inc_h_nless1(k-1))] / ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) (slope_datap)[0+(__inc_slope_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-2))]
1395			/ ((h_datap)[0+(__inc_h_nless1(k-2))] + (h_datap)[0+(__inc_h_nless1(k-1))]) * (slope_datap)[0+(__inc_slope_nless1*(k-1))];
1396			wtave[1] = (h_datap)[0+(__inc_h_nless1(k))] / ((h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) (slope_datap)[0+(__inc_slope_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k-1))] / (
1397			(h_datap)[0+(__inc_h_nless1(k-1))] + (h_datap)[0+(__inc_h_nless1(k))]) * (slope_datap)[0+(__inc_slope_nless1*(k))];
1398			}
1399			/* ....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM */
1400			/* ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE-- */
1401			/* WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K), */
1402			/* IF K.GT.1 */
1403			/* WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1), */
1404			/* IF K.LT.N-1 */
1405			slmax = PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-1))]);
1406			if (k > 1) {
1407			/* Computing MAX */
1408			slmax = PDLMAX(slmax,PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-2))]));
1409			}
1410			if (k < nless1) {
1411			/* Computing MAX */
1412			slmax = PDLMAX(slmax,PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k))]));
1413			}
1414			if (k > 1) {
1415			del[0] = (slope_datap)[0+(__inc_slope_nless1*(k-2))] / slmax;
1416			}
1417			del[1] = (slope_datap)[0+(__inc_slope_nless1*(k-1))] / slmax;
1418			if (k < nless1) {
1419			del[2] = (slope_datap)[0+(__inc_slope_nless1*(k))] / slmax;
1420			}
1421			if (k > 1 && k < nless1) {
1422			/* NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL. */
1423			fact = fudge * PDL_ABS(del[2] * (del[0] - del[1]) * (wtave[1] / slmax));
1424			(d_datap)[0+(__inc_d_n(k-1))] += PDLMIN(fact,1.) (wtave[0] - (d_datap)[0+(__inc_d_n*(k-1))]);
1425			fact = fudge * PDL_ABS(del[0] * (del[2] - del[1]) * (wtave[0] / slmax));
1426			(d_datap)[0+(__inc_d_n(k))] += PDLMIN(fact,1.) (wtave[1] - (d_datap)[0+(__inc_d_n*(k))]);
1427			} else {
1428			/* SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR */
1429			/* K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1). */
1430			fact = fudge * PDL_ABS(del[1]);
1431			(d_datap)[0+(__inc_d_n(nless1))] = PDLMIN(fact,1.) wtave[nless1+1 - k];
1432			/* NOTE THAT I-K+1 = 1 IF K=I (=NLESS1), */
1433			/* I-K+1 = 2 IF K=I-1(=1). */
1434			}
1435			/* ....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA. */
1436			if ((mflag_datap)[0] <= 0.) {
1437			continue;
1438			}
1439			dfloc = (h_datap)[0+(__inc_h_nless1(k-1))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-1))]);
1440			if (k > 1) {
1441			/* Computing MAX */
1442			dfloc = PDLMAX(dfloc,(h_datap)[0+(__inc_h_nless1(k-2))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k-2))]));
1443			}
1444			if (k < nless1) {
1445			/* Computing MAX */
1446			dfloc = PDLMAX(dfloc,(h_datap)[0+(__inc_h_nless1(k))] PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(k))]));
1447			}
1448			dfmx = (mflag_datap)[0] * dfloc;
1449			PDL_Indx indx = nless1 - k;
1450			/* INDX = 1 IF K=I, 2 IF K=I-1. */
1451			/* --------------------------------------------------------------- */
1452			do { /* inline dpchsw */
1453			/* NOTATION AND GENERAL REMARKS. */
1454			/* RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED. */
1455			/* LAMBDA IS THE RATIO OF D2 TO D1. */
1456			/* THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM. */
1457			/* PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO), */
1458			/* WHERE THAT = (XHAT - X1)/H . */
1459			/* THAT IS, P(XHAT)-F1 = DHPHI, WHERE D=D1 OR D2. */
1460			/* SIMILARLY, P(XHAT)-F2 = DH(PHI-RHO) . */
1461			/* Local variables */
1462			PDL_LDouble cp, nu, phi, rho, hphi, that, sigma, small;
1463			PDL_LDouble lambda, radcal;
1464			PDL_LDouble d1 = (d_datap)[0+(__inc_d_n(k-1))], d2 = (d_datap)[0+(__inc_d_n(k))], h2 = (h_datap)[0+(__inc_h_nless1(k-1))], slope2 = (slope_datap)[0+(__inc_slope_nless1(k-1))];
1465			/* Initialized data */
1466			static const PDL_LDouble fact = 100.;
1467			/* THIRD SHOULD BE SLIGHTLY LESS THAN 1/3. */
1468			static const PDL_LDouble third = .33333;
1469			/* SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS. */
1470			small = fact * d1mach;
1471			/* DO MAIN CALCULATION. */
1472			if (d1 == 0.) {
1473			/* SPECIAL CASE -- D1.EQ.ZERO . */
1474			/* IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED. */
1475			if (d2 == 0.) {
1476			(ierr_datap)[0] = -1;
1477			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "D1 AND/OR D2 INVALID");
1478			}
1479			rho = slope2 / d2;
1480			/* EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . */
1481			if (rho >= third) {
1482			(ierr_datap)[0] = 0; break;
1483			}
1484			that = 2. * (3. * rho - 1.) / (3. * (2. * rho - 1.));
1485			/* Computing 2nd power */
1486			phi = that * that * ((3. * rho - 1.) / 3.);
1487			/* CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . */
1488			if (indx != 3) {
1489			phi -= rho;
1490			}
1491			/* TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. */
1492			hphi = h2 * PDL_ABS(phi);
1493			if (hphi * PDL_ABS(d2) > dfmx) {
1494			/* AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK. */
1495			d2 = ((d2) >= 0 ? PDL_ABS(dfmx / hphi) : -PDL_ABS(dfmx / hphi));
1496			}
1497			} else {
1498			rho = slope2 / d1;
1499			lambda = -(d2) / d1;
1500			if (d2 == 0.) {
1501			/* SPECIAL CASE -- D2.EQ.ZERO . */
1502			/* EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 . */
1503			if (rho >= third) {
1504			(ierr_datap)[0] = 0; break;
1505			}
1506			cp = 2. - 3. * rho;
1507			nu = 1. - 2. * rho;
1508			that = 1. / (3. * nu);
1509			} else {
1510			if (lambda <= 0.) {
1511			(ierr_datap)[0] = -1;
1512			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "D1 AND/OR D2 INVALID");
1513			}
1514			/* NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS. */
1515			nu = 1. - lambda - 2. * rho;
1516			sigma = 1. - rho;
1517			cp = nu + sigma;
1518			if (PDL_ABS(nu) > small) {
1519			/* Computing 2nd power */
1520			radcal = (nu - (2. * rho + 1.)) * nu + sigma * sigma;
1521			if (radcal < 0.) {
1522			(ierr_datap)[0] = -2;
1523			return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "NEGATIVE RADICAL");
1524			}
1525			that = (cp - sqrt(radcal)) / (3. * nu);
1526			} else {
1527			that = 1. / (2. * sigma);
1528			}
1529			}
1530			phi = that * ((nu * that - cp) * that + 1.);
1531			/* CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 . */
1532			if (indx != 3) {
1533			phi -= rho;
1534			}
1535			/* TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY. */
1536			hphi = h2 * PDL_ABS(phi);
1537			if (hphi * PDL_ABS(d1) > dfmx) {
1538			/* AT THIS POINT, HPHI.GT.0, SO DIVIDE IS OK. */
1539			d1 = ((d1) >= 0 ? PDL_ABS(dfmx / hphi) : -PDL_ABS(dfmx / hphi));
1540			d2 = -lambda * d1;
1541			}
1542			}
1543			(ierr_datap)[0] = 0;
1544			} while (0); /* end inline dpchsw */
1545			/* --------------------------------------------------------------- */
1546			if ((ierr_datap)[0] != 0) {
1547			break;
1548			}
1549			}} /* Close nless1=1 / / ....... END OF SEGMENT LOOP. */
1550			#line 4884 "lib/PDL/Primitive.pd"
1551			} while (0); /* end inline dpchcs */
1552			}
1553			}
1554			/* SET END CONDITIONS. */
1555			if (ibeg == 0 && iend == 0)
1556			continue;
1557			/* ------------------------------------------------------- */
1558			do { /* inline dpchce */
1559			/* Local variables */
1560			PDL_Indx j, k, ibeg = (ic_datap)[0+(__inc_ic_two(0))], iend = (ic_datap)[0+(__inc_ic_two(1))];
1561			PDL_LDouble stemp[3], xtemp[4];
1562			/* SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. */
1563			if (PDL_ABS(ibeg) > n)
1564			ibeg = 0;
1565			if (PDL_ABS(iend) > n)
1566			iend = 0;
1567			/* TREAT BEGINNING BOUNDARY CONDITION. */
1568			if (ibeg != 0) {
1569			k = PDL_ABS(ibeg);
1570			if (k == 1) {
1571			/* BOUNDARY VALUE PROVIDED. */
1572			(d_datap)[0+(__inc_d_n(0))] = (vc_datap)[0+(__inc_vc_two(0))];
1573			} else if (k == 2) {
1574			/* BOUNDARY SECOND DERIVATIVE PROVIDED. */
1575			(d_datap)[0+(__inc_d_n(0))] = 0.5 (3. * (slope_datap)[0+(__inc_slope_nless1(0))] - (d_datap)[0+(__inc_d_n(1))] - 0.5 * (vc_datap)[0+(__inc_vc_two(0))] (h_datap)[0+(__inc_h_nless1*(0))]);
1576			} else if (k < 5) {
1577			/* USE K-POINT DERIVATIVE FORMULA. */
1578			/* PICK UP FIRST K POINTS, IN REVERSE ORDER. */
1579			for (j = 0; j < k; ++j) {
1580			PDL_Indx index = k - j;
1581			/* INDEX RUNS FROM K DOWN TO 1. */
1582			xtemp[j] = (x_datap)[0+(__inc_x_n*(index+1))];
1583			if (j < k-1) {
1584			stemp[j] = (slope_datap)[0+(__inc_slope_nless1*(index))];
1585			}
1586			}
1587			/* ----------------------------- */
1588
1589			/* PDL version: K, X, S are var names, 4th param output */
1590			/* **PURPOSE Computes divided differences for DPCHCE and DPCHSP /
1591			/* DPCHDF: DPCHIP Finite Difference Formula */
1592			/* Uses a divided difference formulation to compute a K-point approx- */
1593			/* imation to the derivative at X(K) based on the data in X and S. */
1594			/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
1595			/* derivative approximations. */
1596			/* ---------------------------------------------------------------------- */
1597			/* On input: */
1598			/* K is the order of the desired derivative approximation. */
1599			/* K must be at least 3 (error return if not). */
1600			/* X contains the K values of the independent variable. */
1601			/* X need not be ordered, but the values MUST be */
1602			/* distinct. (Not checked here.) */
1603			/* S contains the associated slope values: */
1604			/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
1605			/* (Note that S need only be of length K-1.) */
1606			/* On return: */
1607			/* S will be destroyed. */
1608			/* IERR will be set to -1 if K.LT.2 . */
1609			/* DPCHDF will be set to the desired derivative approximation if */
1610			/* IERR=0 or to zero if IERR=-1. */
1611			/* ---------------------------------------------------------------------- */
1612			/* **SEE ALSO DPCHCE, DPCHSP /
1613			/* **REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- /
1614			/* Verlag, New York, 1978, pp. 10-16. */
1615			/* CHECK FOR LEGAL VALUE OF K. */
1616			{
1617			/* Local variables */
1618			PDL_Indx i, j, k_cached = k;
1619			PDL_LDouble x = xtemp, s = stemp;
1620			if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "K LESS THAN THREE");
1621			/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
1622			for (j = 2; j < k_cached; ++j) {
1623			PDL_Indx itmp = k_cached - j;
1624			for (i = 0; i < itmp; ++i)
1625			s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
1626			}
1627			/* EVALUATE DERIVATIVE AT X(K). */
1628			PDL_LDouble value = s[0];
1629			for (i = 1; i < k_cached-1; ++i)
1630			value = s[i] + value * (x[k_cached-1] - x[i]);
1631			(d_datap)[0+(__inc_d_n*(0))] = value;
1632			}
1633			;
1634			/* ----------------------------- */
1635			} else {
1636			/* USE 'NOT A KNOT' CONDITION. */
1637			(d_datap)[0+(__inc_d_n(0))] = (3. ((h_datap)[0+(__inc_h_nless1(0))] (slope_datap)[0+(__inc_slope_nless1(1))] + (h_datap)[0+(__inc_h_nless1(1))] * (slope_datap)[0+(__inc_slope_nless1*(0))]) -
1638			2. * ((h_datap)[0+(__inc_h_nless1(0))] + (h_datap)[0+(__inc_h_nless1(1))]) * (d_datap)[0+(__inc_d_n(1))] - (h_datap)[0+(__inc_h_nless1(0))] * (d_datap)[0+(__inc_d_n(2))]) / (h_datap)[0+(__inc_h_nless1(1))];
1639			}
1640			/* CHECK D(1,1) FOR COMPATIBILITY WITH MONOTONICITY. */
1641			if (ibeg <= 0) {
1642			if ((slope_datap)[0+(__inc_slope_nless1*(0))] == 0.) {
1643			if ((d_datap)[0+(__inc_d_n*(0))] != 0.) {
1644			(d_datap)[0+(__inc_d_n*(0))] = 0.;
1645			++((ierr_datap)[0]);
1646			}
1647			} else if (((((d_datap)[0+(__inc_d_n(0))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(0))]) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(0))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
1648			(d_datap)[0+(__inc_d_n*(0))] = 0.;
1649			++((ierr_datap)[0]);
1650			} else if (PDL_ABS((d_datap)[0+(__inc_d_n(0))]) > 3. PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(0))])) {
1651			(d_datap)[0+(__inc_d_n(0))] = 3. (slope_datap)[0+(__inc_slope_nless1*(0))];
1652			++((ierr_datap)[0]);
1653			}
1654			}
1655			}
1656			/* TREAT END BOUNDARY CONDITION. */
1657			if (iend == 0)
1658			break;
1659			k = PDL_ABS(iend);
1660			if (k == 1) {
1661			/* BOUNDARY VALUE PROVIDED. */
1662			(d_datap)[0+(__inc_d_n(n-1))] = (vc_datap)[0+(__inc_vc_two(1))];
1663			} else if (k == 2) {
1664			/* BOUNDARY SECOND DERIVATIVE PROVIDED. */
1665			(d_datap)[0+(__inc_d_n(n-1))] = 0.5 (3. * (slope_datap)[0+(__inc_slope_nless1(n-2))] - (d_datap)[0+(__inc_d_n(n-2))]
1666			+ 0.5 * (vc_datap)[0+(__inc_vc_two(1))] (h_datap)[0+(__inc_h_nless1*(n-2))]);
1667			} else if (k < 5) {
1668			/* USE K-POINT DERIVATIVE FORMULA. */
1669			/* PICK UP LAST K POINTS. */
1670			for (j = 0; j < k; ++j) {
1671			PDL_Indx index = n - k + j;
1672			/* INDEX RUNS FROM N+1-K UP TO N. */
1673			xtemp[j] = (x_datap)[0+(__inc_x_n*(index))];
1674			if (j < k-1) {
1675			stemp[j] = (slope_datap)[0+(__inc_slope_nless1*(index))];
1676			}
1677			}
1678			/* ----------------------------- */
1679
1680			/* PDL version: K, X, S are var names, 4th param output */
1681			/* **PURPOSE Computes divided differences for DPCHCE and DPCHSP /
1682			/* DPCHDF: DPCHIP Finite Difference Formula */
1683			/* Uses a divided difference formulation to compute a K-point approx- */
1684			/* imation to the derivative at X(K) based on the data in X and S. */
1685			/* Called by DPCHCE and DPCHSP to compute 3- and 4-point boundary */
1686			/* derivative approximations. */
1687			/* ---------------------------------------------------------------------- */
1688			/* On input: */
1689			/* K is the order of the desired derivative approximation. */
1690			/* K must be at least 3 (error return if not). */
1691			/* X contains the K values of the independent variable. */
1692			/* X need not be ordered, but the values MUST be */
1693			/* distinct. (Not checked here.) */
1694			/* S contains the associated slope values: */
1695			/* S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. */
1696			/* (Note that S need only be of length K-1.) */
1697			/* On return: */
1698			/* S will be destroyed. */
1699			/* IERR will be set to -1 if K.LT.2 . */
1700			/* DPCHDF will be set to the desired derivative approximation if */
1701			/* IERR=0 or to zero if IERR=-1. */
1702			/* ---------------------------------------------------------------------- */
1703			/* **SEE ALSO DPCHCE, DPCHSP /
1704			/* **REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- /
1705			/* Verlag, New York, 1978, pp. 10-16. */
1706			/* CHECK FOR LEGAL VALUE OF K. */
1707			{
1708			/* Local variables */
1709			PDL_Indx i, j, k_cached = k;
1710			PDL_LDouble x = xtemp, s = stemp;
1711			if (k_cached < 3) return PDL->make_error(PDL_EUSERERROR, "Error in pchip_chic:" "K LESS THAN THREE");
1712			/* COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. */
1713			for (j = 2; j < k_cached; ++j) {
1714			PDL_Indx itmp = k_cached - j;
1715			for (i = 0; i < itmp; ++i)
1716			s[i] = (s[i+1] - s[i]) / (x[i + j] - x[i]);
1717			}
1718			/* EVALUATE DERIVATIVE AT X(K). */
1719			PDL_LDouble value = s[0];
1720			for (i = 1; i < k_cached-1; ++i)
1721			value = s[i] + value * (x[k_cached-1] - x[i]);
1722			(d_datap)[0+(__inc_d_n*(n-1))] = value;
1723			}
1724			;
1725			/* ----------------------------- */
1726			} else {
1727			/* USE 'NOT A KNOT' CONDITION. */
1728			(d_datap)[0+(__inc_d_n(n-1))] = (3. ((h_datap)[0+(__inc_h_nless1(n-2))] (slope_datap)[0+(__inc_slope_nless1*(n-3))] +
1729			(h_datap)[0+(__inc_h_nless1(n-3))] (slope_datap)[0+(__inc_slope_nless1(n-2))]) - 2. ((h_datap)[0+(__inc_h_nless1(n-2))] + (h_datap)[0+(__inc_h_nless1(n-3))]) *
1730			(d_datap)[0+(__inc_d_n(n-2))] - (h_datap)[0+(__inc_h_nless1(n-2))] * (d_datap)[0+(__inc_d_n(n-3))]) / (h_datap)[0+(__inc_h_nless1(n-3))];
1731			}
1732			if (iend > 0)
1733			break;
1734			/* CHECK D(1,N) FOR COMPATIBILITY WITH MONOTONICITY. */
1735			if ((slope_datap)[0+(__inc_slope_nless1*(n-2))] == 0.) {
1736			if ((d_datap)[0+(__inc_d_n*(n-1))] != 0.) {
1737			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
1738			(ierr_datap)[0] += 2;
1739			}
1740			} else if (((((d_datap)[0+(__inc_d_n(n-1))]) == 0. \|\| ((slope_datap)[0+(__inc_slope_nless1(n-2))]) == 0.) ? 0. : ((((d_datap)[0+(__inc_d_n(n-1))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1)) ((((slope_datap)[0+(__inc_slope_nless1*(n-2))])) >= 0 ? PDL_ABS(1) : -PDL_ABS(1))) < 0.) {
1741			(d_datap)[0+(__inc_d_n*(n-1))] = 0.;
1742			(ierr_datap)[0] += 2;
1743			} else if (PDL_ABS((d_datap)[0+(__inc_d_n(n-1))]) > 3. PDL_ABS((slope_datap)[0+(__inc_slope_nless1*(n-2))])) {
1744			(d_datap)[0+(__inc_d_n(n-1))] = 3. (slope_datap)[0+(__inc_slope_nless1*(n-2))];
1745			(ierr_datap)[0] += 2;
1746			}
1747			} while (0); /* end inlined dpchce */
1748			/* ------------------------------------------------------- */
1749			#line 1750 "lib/PDL/Primitive-pp-pchip_chic.c"
1750	0	0	}PDL_BROADCASTLOOP_END_pchip_chic_readdata
		0
1751	0		} break;
1752	0		default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in pchip_chic: unhandled datatype(%d), only handles (FDE)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
1753			}
1754			#undef PDL_IF_BAD
1755	3		return PDL_err;
1756			}
1757
1758			static pdl_datatypes pdl_pchip_chic_vtable_gentypes[] = { PDL_F, PDL_D, PDL_LD, -1 };
1759			static PDL_Indx pdl_pchip_chic_vtable_realdims[] = { 1, 1, 0, 1, 1, 1, 0, 1, 1 };
1760			static char *pdl_pchip_chic_vtable_parnames[] = { "ic","vc","mflag","x","f","d","ierr","h","slope" };
1761			static short pdl_pchip_chic_vtable_parflags[] = {
1762			PDL_PARAM_ISTYPED,
1763			0,
1764			0,
1765			0,
1766			0,
1767			PDL_PARAM_ISCREAT\|PDL_PARAM_ISOUT\|PDL_PARAM_ISWRITE,
1768			PDL_PARAM_ISCREAT\|PDL_PARAM_ISOUT\|PDL_PARAM_ISTYPED\|PDL_PARAM_ISWRITE,
1769			PDL_PARAM_ISCREAT\|PDL_PARAM_ISTEMP\|PDL_PARAM_ISWRITE,
1770			PDL_PARAM_ISCREAT\|PDL_PARAM_ISTEMP\|PDL_PARAM_ISWRITE
1771			};
1772			static pdl_datatypes pdl_pchip_chic_vtable_partypes[] = { PDL_SB, -1, -1, -1, -1, -1, PDL_IND, -1, -1 };
1773			static PDL_Indx pdl_pchip_chic_vtable_realdims_starts[] = { 0, 1, 2, 2, 3, 4, 5, 5, 6 };
1774			static PDL_Indx pdl_pchip_chic_vtable_realdims_ind_ids[] = { 2, 2, 0, 0, 0, 1, 1 };
1775			static char *pdl_pchip_chic_vtable_indnames[] = { "n","nless1","two" };
1776			pdl_transvtable pdl_pchip_chic_vtable = {
1777			PDL_TRANS_DO_BROADCAST, 0, pdl_pchip_chic_vtable_gentypes, 5, 9, NULL /CORE21/,
1778			pdl_pchip_chic_vtable_realdims, pdl_pchip_chic_vtable_parnames,
1779			pdl_pchip_chic_vtable_parflags, pdl_pchip_chic_vtable_partypes,
1780			pdl_pchip_chic_vtable_realdims_starts, pdl_pchip_chic_vtable_realdims_ind_ids, 7,
1781			3, pdl_pchip_chic_vtable_indnames,
1782			pdl_pchip_chic_redodims, pdl_pchip_chic_readdata, NULL,
1783			NULL,
1784			0,"PDL::Primitive::pchip_chic"
1785			};
1786
1787
1788	3		pdl_error pdl_run_pchip_chic(pdl ic,pdl vc,pdl mflag,pdl x,pdl f,pdl d,pdl *ierr) {
1789	3		pdl_error PDL_err = {0, NULL, 0};
1790	3	50	if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
1791	3		pdl_trans *__privtrans = PDL->create_trans(&pdl_pchip_chic_vtable);
1792	3	50	if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
1793	3		__privtrans->pdls[0] = ic;
1794	3		__privtrans->pdls[1] = vc;
1795	3		__privtrans->pdls[2] = mflag;
1796	3		__privtrans->pdls[3] = x;
1797	3		__privtrans->pdls[4] = f;
1798	3		__privtrans->pdls[5] = d;
1799	3		__privtrans->pdls[6] = ierr;
1800	3	50	PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
1801	3	50	PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
1802	3		return PDL_err;
1803			}