Mercurial > hg > camir-aes2014
comparison toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/divide_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 /* divide_by_sparse_table.c ../potential/tables*/ | |
2 | |
3 /******************************************/ | |
4 /* 6 input & 1 output */ | |
5 /* Big table [0] */ | |
6 /* Big domain [1] */ | |
7 /* big sizes [2] */ | |
8 /* Small table [3] */ | |
9 /* small domain [4] */ | |
10 /* small sizes [5] */ | |
11 /* */ | |
12 /* New big table[0] */ | |
13 /******************************************/ | |
14 | |
15 #include <math.h> | |
16 #include <stdlib.h> | |
17 #include "mex.h" | |
18 | |
19 int compare(const void* src1, const void* src2){ | |
20 int i1 = *(int*)src1 ; | |
21 int i2 = *(int*)src2 ; | |
22 return i1-i2 ; | |
23 } | |
24 | |
25 void ind_subv(int index, const int *cumprod, int n, int *bsubv){ | |
26 int i; | |
27 | |
28 for (i = n-1; i >= 0; i--) { | |
29 bsubv[i] = ((int)floor(index / cumprod[i])); | |
30 index = index % cumprod[i]; | |
31 } | |
32 } | |
33 | |
34 int subv_ind(const int n, const int *cumprod, const int *subv){ | |
35 int i, index=0; | |
36 | |
37 for(i=0; i<n; i++){ | |
38 index += subv[i] * cumprod[i]; | |
39 } | |
40 return index; | |
41 } | |
42 | |
43 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ | |
44 int i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex; | |
45 int *mask, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv; | |
46 double *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value; | |
47 | |
48 plhs[0] = mxDuplicateArray(prhs[0]); | |
49 pbDomain = mxGetPr(prhs[1]); | |
50 bdim = mxGetNumberOfElements(prhs[1]); | |
51 psDomain = mxGetPr(prhs[4]); | |
52 sdim = mxGetNumberOfElements(prhs[4]); | |
53 | |
54 pbSize = mxGetPr(prhs[2]); | |
55 psSize = mxGetPr(prhs[5]); | |
56 | |
57 NB = 1; | |
58 for(i=0; i<bdim; i++){ | |
59 NB *= (int)pbSize[i]; | |
60 } | |
61 | |
62 bpr = mxGetPr(plhs[0]); | |
63 bir = mxGetIr(plhs[0]); | |
64 bjc = mxGetJc(plhs[0]); | |
65 NZB = bjc[1]; | |
66 | |
67 spr = mxGetPr(prhs[3]); | |
68 sir = mxGetIr(prhs[3]); | |
69 sjc = mxGetJc(prhs[3]); | |
70 NZS = sjc[1]; | |
71 | |
72 if(sdim == 0){ | |
73 value = *spr; | |
74 if(value == 0)value = 1; | |
75 for(i=0; i<NZB; i++){ | |
76 bpr[i] /= value; | |
77 } | |
78 return; | |
79 } | |
80 | |
81 mask = malloc(sdim * sizeof(int)); | |
82 bCumprod = malloc(bdim * sizeof(int)); | |
83 sCumprod = malloc(sdim * sizeof(int)); | |
84 bsubv = malloc(bdim * sizeof(int)); | |
85 ssubv = malloc(sdim * sizeof(int)); | |
86 | |
87 count = 0; | |
88 for(i=0; i<sdim; i++){ | |
89 for(j=0; j<bdim; j++){ | |
90 if(psDomain[i] == pbDomain[j]){ | |
91 mask[count] = j; | |
92 count++; | |
93 break; | |
94 } | |
95 } | |
96 } | |
97 | |
98 bCumprod[0] = 1; | |
99 for(i=0; i<bdim-1; i++){ | |
100 bCumprod[i+1] = bCumprod[i] * (int)pbSize[i]; | |
101 } | |
102 sCumprod[0] = 1; | |
103 for(i=0; i<sdim-1; i++){ | |
104 sCumprod[i+1] = sCumprod[i] * (int)psSize[i]; | |
105 } | |
106 | |
107 for(i=0; i<NZB; i++){ | |
108 bindex = bir[i]; | |
109 ind_subv(bindex, bCumprod, bdim, bsubv); | |
110 for(j=0; j<sdim; j++){ | |
111 ssubv[j] = bsubv[mask[j]]; | |
112 } | |
113 sindex = subv_ind(sdim, sCumprod, ssubv); | |
114 result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare); | |
115 if(result){ | |
116 position = result - sir; | |
117 bpr[i] /= spr[position]; | |
118 } | |
119 } | |
120 | |
121 free(mask); | |
122 free(bCumprod); | |
123 free(sCumprod); | |
124 free(bsubv); | |
125 free(ssubv); | |
126 } |