Mercurial > hg > camir-aes2014
comparison toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/mult_by_sparse_table.c @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:e9a9cd732c1e |
---|---|
1 /* mult_by_sparse_table.c ../potential/tables*/ | |
2 | |
3 | |
4 /******************************************/ | |
5 /* 6 input & 1 output */ | |
6 /* Big table [0] */ | |
7 /* Big domain [1] */ | |
8 /* big sizes [2] */ | |
9 /* Small table [3] */ | |
10 /* small domain [4] */ | |
11 /* small sizes [5] */ | |
12 /* */ | |
13 /* New big table[0] */ | |
14 /******************************************/ | |
15 | |
16 #include <math.h> | |
17 #include <stdlib.h> | |
18 #include "mex.h" | |
19 | |
20 int compare(const void* src1, const void* src2){ | |
21 int i1 = *(int*)src1 ; | |
22 int i2 = *(int*)src2 ; | |
23 return i1-i2 ; | |
24 } | |
25 | |
26 void ind_subv(int index, const int *cumprod, int n, int *bsubv){ | |
27 int i; | |
28 | |
29 for (i = n-1; i >= 0; i--) { | |
30 bsubv[i] = ((int)floor(index / cumprod[i])); | |
31 index = index % cumprod[i]; | |
32 } | |
33 } | |
34 | |
35 int subv_ind(const int n, const int *cumprod, const int *subv){ | |
36 int i, index=0; | |
37 | |
38 for(i=0; i<n; i++){ | |
39 index += subv[i] * cumprod[i]; | |
40 } | |
41 return index; | |
42 } | |
43 | |
44 void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){ | |
45 double *ptr; | |
46 void *newptr; | |
47 int *ir, *jc; | |
48 int nbytes; | |
49 | |
50 if(new_nzmax == old_nzmax) return; | |
51 nbytes = new_nzmax * sizeof(*ptr); | |
52 ptr = mxGetPr(spArray); | |
53 newptr = mxRealloc(ptr, nbytes); | |
54 mxSetPr(spArray, newptr); | |
55 nbytes = new_nzmax * sizeof(*ir); | |
56 ir = mxGetIr(spArray); | |
57 newptr = mxRealloc(ir, nbytes); | |
58 mxSetIr(spArray, newptr); | |
59 jc = mxGetJc(spArray); | |
60 jc[0] = 0; | |
61 jc[1] = new_nzmax; | |
62 mxSetNzmax(spArray, new_nzmax); | |
63 } | |
64 | |
65 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ | |
66 int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0; | |
67 int *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv; | |
68 double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr; | |
69 | |
70 pbDomain = mxGetPr(prhs[1]); | |
71 bdim = mxGetNumberOfElements(prhs[1]); | |
72 psDomain = mxGetPr(prhs[4]); | |
73 sdim = mxGetNumberOfElements(prhs[4]); | |
74 | |
75 pbSize = mxGetPr(prhs[2]); | |
76 psSize = mxGetPr(prhs[5]); | |
77 | |
78 NB = 1; | |
79 for(i=0; i<bdim; i++){ | |
80 NB *= (int)pbSize[i]; | |
81 } | |
82 | |
83 bpr = mxGetPr(prhs[0]); | |
84 bir = mxGetIr(prhs[0]); | |
85 bjc = mxGetJc(prhs[0]); | |
86 NZB = bjc[1]; | |
87 | |
88 spr = mxGetPr(prhs[3]); | |
89 sir = mxGetIr(prhs[3]); | |
90 sjc = mxGetJc(prhs[3]); | |
91 NZS = sjc[1]; | |
92 | |
93 plhs[0] = mxDuplicateArray(prhs[0]); | |
94 rpr = mxGetPr(plhs[0]); | |
95 rir = mxGetIr(plhs[0]); | |
96 rjc = mxGetJc(plhs[0]); | |
97 rjc[0] = 0; | |
98 rjc[1] = NZB; | |
99 | |
100 if(sdim == 0){ | |
101 for(i=0; i<NZB; i++){ | |
102 rpr[i] *= *spr; | |
103 } | |
104 return; | |
105 } | |
106 | |
107 mask = malloc(sdim * sizeof(int)); | |
108 bCumprod = malloc(bdim * sizeof(int)); | |
109 sCumprod = malloc(sdim * sizeof(int)); | |
110 bsubv = malloc(bdim * sizeof(int)); | |
111 ssubv = malloc(sdim * sizeof(int)); | |
112 | |
113 count = 0; | |
114 for(i=0; i<sdim; i++){ | |
115 for(j=0; j<bdim; j++){ | |
116 if(psDomain[i] == pbDomain[j]){ | |
117 mask[count] = j; | |
118 count++; | |
119 break; | |
120 } | |
121 } | |
122 } | |
123 | |
124 bCumprod[0] = 1; | |
125 for(i=0; i<bdim-1; i++){ | |
126 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; | |
127 } | |
128 sCumprod[0] = 1; | |
129 for(i=0; i<sdim-1; i++){ | |
130 sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; | |
131 } | |
132 | |
133 for(i=0; i<NZB; i++){ | |
134 bindex = bir[i]; | |
135 ind_subv(bindex, bCumprod, bdim, bsubv); | |
136 for(j=0; j<sdim; j++){ | |
137 ssubv[j] = bsubv[mask[j]]; | |
138 } | |
139 sindex = subv_ind(sdim, sCumprod, ssubv); | |
140 result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); | |
141 if(result){ | |
142 position = result - sir; | |
143 rpr[nzCounts] = bpr[i] * spr[position]; | |
144 rir[nzCounts] = bindex; | |
145 nzCounts++; | |
146 } | |
147 } | |
148 | |
149 reset_nzmax(plhs[0], NZB, nzCounts); | |
150 free(mask); | |
151 free(bCumprod); | |
152 free(sCumprod); | |
153 free(bsubv); | |
154 free(ssubv); | |
155 } |