File Coverage

lib/PDL/ImageND-pp-repulse.c
Criterion Covered Total %
statement 84 138 60.8
branch 49 138 35.5
condition n/a
subroutine n/a
pod n/a
total 133 276 48.1


line stmt bran cond sub pod time code
1              
2             #line 453 "lib/PDL/PP.pm"
3             /*
4             * THIS FILE WAS GENERATED BY PDL::PP from lib/PDL/ImageND.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_ImageND
21             extern Core* PDL; /* Structure hold core C functions */
22              
23             #line 1846 "lib/PDL/PP.pm"
24             typedef struct pdl_params_repulse {
25             #line 26 "lib/PDL/ImageND-pp-repulse.c"
26             double boxsize;
27             int dmult;
28             double a;
29             double b;
30             double c;
31             double d;
32             } pdl_params_repulse;
33              
34              
35             #line 1857 "lib/PDL/PP.pm"
36             pdl_error pdl_repulse_readdata(pdl_trans *__privtrans) {
37             pdl_error PDL_err = {0, NULL, 0};
38             #line 39 "lib/PDL/ImageND-pp-repulse.c"
39 1           pdl_params_repulse *__params = __privtrans->params; (void)__params;
40 1           register PDL_Indx __nc_size = __privtrans->ind_sizes[0];
41 1           register PDL_Indx __np_size = __privtrans->ind_sizes[1];
42 1 50         if (!__privtrans->broadcast.incs) return PDL->make_error(PDL_EUSERERROR, "Error in repulse:" "broadcast.incs NULL");
43             /* broadcastloop declarations */
44             int __brcloopval;
45             register PDL_Indx __tind0,__tind1; /* counters along dim */
46 1           register PDL_Indx __tnpdls = __privtrans->broadcast.npdls;
47             /* dims here are how many steps along those dims */
48 1           register PDL_Indx __tinc0_coords = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,0);
49 1           register PDL_Indx __tinc0_vecs = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,0);
50 1           register PDL_Indx __tinc0_links = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,0);
51 1           register PDL_Indx __tinc1_coords = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,0,1);
52 1           register PDL_Indx __tinc1_vecs = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,1,1);
53 1           register PDL_Indx __tinc1_links = PDL_BRC_INC(__privtrans->broadcast.incs,__tnpdls,2,1);
54             #define PDL_BROADCASTLOOP_START_repulse_readdata PDL_BROADCASTLOOP_START( \
55             readdata, \
56             __privtrans->broadcast, \
57             __privtrans->vtable, \
58             coords_datap += __offsp[0]; \
59             vecs_datap += __offsp[1]; \
60             links_datap += __offsp[2]; \
61             , \
62             ( ,coords_datap += __tinc1_coords - __tinc0_coords * __tdims0 \
63             ,vecs_datap += __tinc1_vecs - __tinc0_vecs * __tdims0 \
64             ,links_datap += __tinc1_links - __tinc0_links * __tdims0 \
65             ), \
66             ( ,coords_datap += __tinc0_coords \
67             ,vecs_datap += __tinc0_vecs \
68             ,links_datap += __tinc0_links \
69             ) \
70             )
71             #define PDL_BROADCASTLOOP_END_repulse_readdata PDL_BROADCASTLOOP_END( \
72             __privtrans->broadcast, \
73             coords_datap -= __tinc1_coords * __tdims1 + __offsp[0]; \
74             vecs_datap -= __tinc1_vecs * __tdims1 + __offsp[1]; \
75             links_datap -= __tinc1_links * __tdims1 + __offsp[2]; \
76             )
77 1           register PDL_Indx __inc_coords_nc = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,0)]; (void)__inc_coords_nc;register PDL_Indx __inc_coords_np = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,0,1)]; (void)__inc_coords_np;
78 1           register PDL_Indx __inc_links_np = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,2,0)]; (void)__inc_links_np;
79 1           register PDL_Indx __inc_vecs_nc = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,0)]; (void)__inc_vecs_nc;register PDL_Indx __inc_vecs_np = __privtrans->inc_sizes[PDL_INC_ID(__privtrans->vtable,1,1)]; (void)__inc_vecs_np;
80             #ifndef PDL_DECLARE_PARAMS_repulse_1
81             #define PDL_DECLARE_PARAMS_repulse_1(PDL_TYPE_OP,PDL_PPSYM_OP,PDL_TYPE_PARAM_links,PDL_PPSYM_PARAM_links) \
82             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, coords, (__privtrans->pdls[0]), 1, PDL_PPSYM_OP) \
83             PDL_DECLARE_PARAMETER(PDL_TYPE_OP, vecs, (__privtrans->pdls[1]), 1, PDL_PPSYM_OP) \
84             PDL_DECLARE_PARAMETER(PDL_TYPE_PARAM_links, links, (__privtrans->pdls[2]), 1, PDL_PPSYM_PARAM_links)
85             #endif
86             #define PDL_IF_BAD(t,f) f
87 1           switch (__privtrans->__datatype) { /* Start generic switch */
88 1           case PDL_F: {
89 1 50         PDL_DECLARE_PARAMS_repulse_1(PDL_Float,F,PDL_Long,L)
    50          
    50          
    50          
    50          
    50          
    50          
    50          
    50          
90 4 50         PDL_BROADCASTLOOP_START_repulse_readdata {double a = __params->a;
    50          
    50          
    50          
    50          
    100          
    100          
91 1           double b = __params->b;
92 1           double c = __params->c;
93 1           double d = __params->d;
94             int ind; int x,y,z;
95 1           HV *hv = newHV();
96 1           double boxsize = __params->boxsize;
97 1           int dmult = __params->dmult;
98 7 100         {/* Open np */ PDL_EXPAND2(register PDL_Indx np=0, __np_stop=(__np_size)); for(; np<__np_stop; np+=1) {
99 6           int index = 0;
100 6           (links_datap)[0+(__inc_links_np*(np))] = -1;
101 24 100         {/* Open nc */ PDL_EXPAND2(register PDL_Indx nc=0, __nc_stop=(__nc_size)); for(; nc<__nc_stop; nc+=1) {
102 18           (vecs_datap)[0+(__inc_vecs_nc*(nc))+(__inc_vecs_np*(np))] = 0;
103 18           index *= dmult;
104 18           index += (int)((coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(np))]/boxsize);
105             }} /* Close nc */
106             /* Repulse old (shame to use x,y,z...) */
107 24 100         for (x=-1; x<=1; x++) {
108 72 100         for (y=-1; y<=1; y++) {
109 216 100         for (z=-1; z<=1; z++) {
110 162           int ni = index + x + dmult * y + dmult * dmult * z;
111 162           SV **svp = hv_fetch(hv, (char *)&ni, sizeof(int), 0);
112 162 100         if (svp && *svp) {
    50          
113 11           ind = SvIV(*svp) - 1;
114 26 100         while (ind>=0) {
115 15           double dist = 0;
116             double dist2;
117             double tmp;
118             double func;
119 60 100         {/* Open nc */ PDL_EXPAND2(register PDL_Indx nc=0, __nc_stop=(__nc_size)); for(; nc<__nc_stop; nc+=1) {
120 45           tmp = ((coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(np))] - (coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(ind))]);
121 45           dist += tmp * tmp;
122             }} /* Close nc */
123 15           dist = sqrt(1/(sqrt(dist)+d));
124 15           func = c * dist;
125 15           dist2 = dist * dist;
126 15           func += b * dist2;
127 15           dist2 *= dist2;
128 15           func += a * dist2;
129 60 100         {/* Open nc */ PDL_EXPAND2(register PDL_Indx nc=0, __nc_stop=(__nc_size)); for(; nc<__nc_stop; nc+=1) {
130 45           tmp = ((coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(np))] - (coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(ind))]);
131 45           (vecs_datap)[0+(__inc_vecs_nc*(nc))+(__inc_vecs_np*(np))] -= func * tmp;
132 45           (vecs_datap)[0+(__inc_vecs_nc*(nc))+(__inc_vecs_np*(ind))] += func * tmp;
133             }} /* Close nc */
134 15           ind = (links_datap)[0+(__inc_links_np*(ind))];
135             }
136             }
137             }
138             }
139             }
140             /* Store new */
141 6           SV **svp = hv_fetch(hv, (char *)&index, sizeof(index), 1);
142 6 50         if(!svp || !*svp)
    50          
143 0           return PDL->make_error(PDL_EUSERERROR, "Error in repulse:" "Invalid sv from hvfetch");
144 6           SV *sv = *svp;
145             int npv;
146 6 100         if(SvOK(sv) && (npv = SvIV(sv))) {
    50          
147 2           npv --;
148 2           (links_datap)[0+(__inc_links_np*(np))] = (links_datap)[0+(__inc_links_np*(npv))];
149 2           (links_datap)[0+(__inc_links_np*(npv))] = np;
150             } else {
151 4           sv_setiv(sv,np+1);
152 4           (links_datap)[0+(__inc_links_np*(np))] = -1;
153             }
154             }} /* Close np */
155 1           hv_undef(hv);
156 1 50         }PDL_BROADCASTLOOP_END_repulse_readdata
    50          
157 1           } break;
158 0           case PDL_D: {
159 0 0         PDL_DECLARE_PARAMS_repulse_1(PDL_Double,D,PDL_Long,L)
    0          
    0          
    0          
    0          
    0          
    0          
    0          
    0          
160 0 0         PDL_BROADCASTLOOP_START_repulse_readdata {double a = __params->a;
    0          
    0          
    0          
    0          
    0          
    0          
161 0           double b = __params->b;
162 0           double c = __params->c;
163 0           double d = __params->d;
164             int ind; int x,y,z;
165 0           HV *hv = newHV();
166 0           double boxsize = __params->boxsize;
167 0           int dmult = __params->dmult;
168 0 0         {/* Open np */ PDL_EXPAND2(register PDL_Indx np=0, __np_stop=(__np_size)); for(; np<__np_stop; np+=1) {
169 0           int index = 0;
170 0           (links_datap)[0+(__inc_links_np*(np))] = -1;
171 0 0         {/* Open nc */ PDL_EXPAND2(register PDL_Indx nc=0, __nc_stop=(__nc_size)); for(; nc<__nc_stop; nc+=1) {
172 0           (vecs_datap)[0+(__inc_vecs_nc*(nc))+(__inc_vecs_np*(np))] = 0;
173 0           index *= dmult;
174 0           index += (int)((coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(np))]/boxsize);
175             }} /* Close nc */
176             /* Repulse old (shame to use x,y,z...) */
177 0 0         for (x=-1; x<=1; x++) {
178 0 0         for (y=-1; y<=1; y++) {
179 0 0         for (z=-1; z<=1; z++) {
180 0           int ni = index + x + dmult * y + dmult * dmult * z;
181 0           SV **svp = hv_fetch(hv, (char *)&ni, sizeof(int), 0);
182 0 0         if (svp && *svp) {
    0          
183 0           ind = SvIV(*svp) - 1;
184 0 0         while (ind>=0) {
185 0           double dist = 0;
186             double dist2;
187             double tmp;
188             double func;
189 0 0         {/* Open nc */ PDL_EXPAND2(register PDL_Indx nc=0, __nc_stop=(__nc_size)); for(; nc<__nc_stop; nc+=1) {
190 0           tmp = ((coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(np))] - (coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(ind))]);
191 0           dist += tmp * tmp;
192             }} /* Close nc */
193 0           dist = sqrt(1/(sqrt(dist)+d));
194 0           func = c * dist;
195 0           dist2 = dist * dist;
196 0           func += b * dist2;
197 0           dist2 *= dist2;
198 0           func += a * dist2;
199 0 0         {/* Open nc */ PDL_EXPAND2(register PDL_Indx nc=0, __nc_stop=(__nc_size)); for(; nc<__nc_stop; nc+=1) {
200 0           tmp = ((coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(np))] - (coords_datap)[0+(__inc_coords_nc*(nc))+(__inc_coords_np*(ind))]);
201 0           (vecs_datap)[0+(__inc_vecs_nc*(nc))+(__inc_vecs_np*(np))] -= func * tmp;
202 0           (vecs_datap)[0+(__inc_vecs_nc*(nc))+(__inc_vecs_np*(ind))] += func * tmp;
203             }} /* Close nc */
204 0           ind = (links_datap)[0+(__inc_links_np*(ind))];
205             }
206             }
207             }
208             }
209             }
210             /* Store new */
211 0           SV **svp = hv_fetch(hv, (char *)&index, sizeof(index), 1);
212 0 0         if(!svp || !*svp)
    0          
213 0           return PDL->make_error(PDL_EUSERERROR, "Error in repulse:" "Invalid sv from hvfetch");
214 0           SV *sv = *svp;
215             int npv;
216 0 0         if(SvOK(sv) && (npv = SvIV(sv))) {
    0          
217 0           npv --;
218 0           (links_datap)[0+(__inc_links_np*(np))] = (links_datap)[0+(__inc_links_np*(npv))];
219 0           (links_datap)[0+(__inc_links_np*(npv))] = np;
220             } else {
221 0           sv_setiv(sv,np+1);
222 0           (links_datap)[0+(__inc_links_np*(np))] = -1;
223             }
224             }} /* Close np */
225 0           hv_undef(hv);
226 0 0         }PDL_BROADCASTLOOP_END_repulse_readdata
    0          
227 0           } break;
228 0           default: return PDL->make_error(PDL_EUSERERROR, "PP INTERNAL ERROR in repulse: unhandled datatype(%d), only handles (FD)! PLEASE MAKE A BUG REPORT\n", __privtrans->__datatype);
229             }
230             #undef PDL_IF_BAD
231 1           return PDL_err;
232             }
233              
234             static pdl_datatypes pdl_repulse_vtable_gentypes[] = { PDL_F, PDL_D, -1 };
235             static PDL_Indx pdl_repulse_vtable_realdims[] = { 2, 2, 1 };
236             static char *pdl_repulse_vtable_parnames[] = { "coords","vecs","links" };
237             static short pdl_repulse_vtable_parflags[] = {
238             0,
239             PDL_PARAM_ISCREAT|PDL_PARAM_ISOUT|PDL_PARAM_ISWRITE,
240             PDL_PARAM_ISCREAT|PDL_PARAM_ISTEMP|PDL_PARAM_ISTYPED|PDL_PARAM_ISWRITE
241             };
242             static pdl_datatypes pdl_repulse_vtable_partypes[] = { -1, -1, PDL_L };
243             static PDL_Indx pdl_repulse_vtable_realdims_starts[] = { 0, 2, 4 };
244             static PDL_Indx pdl_repulse_vtable_realdims_ind_ids[] = { 0, 1, 0, 1, 1 };
245             static char *pdl_repulse_vtable_indnames[] = { "nc","np" };
246             pdl_transvtable pdl_repulse_vtable = {
247             PDL_TRANS_DO_BROADCAST, 0, pdl_repulse_vtable_gentypes, 1, 3, NULL /*CORE21*/,
248             pdl_repulse_vtable_realdims, pdl_repulse_vtable_parnames,
249             pdl_repulse_vtable_parflags, pdl_repulse_vtable_partypes,
250             pdl_repulse_vtable_realdims_starts, pdl_repulse_vtable_realdims_ind_ids, 5,
251             2, pdl_repulse_vtable_indnames,
252             NULL, pdl_repulse_readdata, NULL,
253             NULL,
254             sizeof(pdl_params_repulse),"PDL::ImageND::repulse"
255             };
256              
257              
258 1           pdl_error pdl_run_repulse(pdl *coords,pdl *vecs,double boxsize,int dmult,double a,double b,double c,double d) {
259 1           pdl_error PDL_err = {0, NULL, 0};
260 1 50         if (!PDL) return (pdl_error){PDL_EFATAL, "PDL core struct is NULL, can't continue",0};
261 1           pdl_trans *__privtrans = PDL->create_trans(&pdl_repulse_vtable);
262 1 50         if (!__privtrans) return PDL->make_error_simple(PDL_EFATAL, "Couldn't create trans");
263 1           pdl_params_repulse *__params = __privtrans->params;
264 1           __privtrans->pdls[0] = coords;
265 1           __privtrans->pdls[1] = vecs;
266 1 50         PDL_RETERROR(PDL_err, PDL->type_coerce(__privtrans));
267 1           (__params->boxsize) = (boxsize); /* CType.get_copy */
268 1           (__params->dmult) = (dmult); /* CType.get_copy */
269 1           (__params->a) = (a); /* CType.get_copy */
270 1           (__params->b) = (b); /* CType.get_copy */
271 1           (__params->c) = (c); /* CType.get_copy */
272 1           (__params->d) = (d); /* CType.get_copy */
273 1 50         PDL_RETERROR(PDL_err, PDL->make_trans_mutual(__privtrans));
274 1           return PDL_err;
275             }