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