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