|
wolffd@0
|
1 /* divide_by_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 "mex.h"
|
|
wolffd@0
|
17
|
|
wolffd@0
|
18 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
|
|
wolffd@0
|
19 int i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
|
|
wolffd@0
|
20 int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
|
|
wolffd@0
|
21 double *pbDomain, *psDomain, *sp, *zp, *bs, value;
|
|
wolffd@0
|
22
|
|
wolffd@0
|
23 plhs[0] = mxDuplicateArray(prhs[0]);
|
|
wolffd@0
|
24 zp = mxGetPr(plhs[0]);
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 siz_b = mxGetNumberOfElements(prhs[1]);
|
|
wolffd@0
|
27 siz_s = mxGetNumberOfElements(prhs[4]);
|
|
wolffd@0
|
28 pbDomain = mxGetPr(prhs[1]);
|
|
wolffd@0
|
29 psDomain = mxGetPr(prhs[4]);
|
|
wolffd@0
|
30
|
|
wolffd@0
|
31 NB = mxGetNumberOfElements(prhs[0]);
|
|
wolffd@0
|
32 NS = mxGetNumberOfElements(prhs[3]);
|
|
wolffd@0
|
33 sp = mxGetPr(prhs[3]);
|
|
wolffd@0
|
34
|
|
wolffd@0
|
35 bs = mxGetPr(prhs[2]);
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37 if(NS == 1){
|
|
wolffd@0
|
38 value = *sp;
|
|
wolffd@0
|
39 if(value == 0) value = 1;
|
|
wolffd@0
|
40 for(i=0; i<NB; i++){
|
|
wolffd@0
|
41 zp[i] /= value;
|
|
wolffd@0
|
42 }
|
|
wolffd@0
|
43 return;
|
|
wolffd@0
|
44 }
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46 if(NS == NB){
|
|
wolffd@0
|
47 for(i=0; i<NB; i++){
|
|
wolffd@0
|
48 value = sp[i];
|
|
wolffd@0
|
49 if(value == 0) value = 1;
|
|
wolffd@0
|
50 zp[i] /= value;
|
|
wolffd@0
|
51 }
|
|
wolffd@0
|
52 return;
|
|
wolffd@0
|
53 }
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 mask = malloc(siz_s * sizeof(int));
|
|
wolffd@0
|
56 count = 0;
|
|
wolffd@0
|
57 for(i=0; i<siz_s; i++){
|
|
wolffd@0
|
58 for(j=0; j<siz_b; j++){
|
|
wolffd@0
|
59 if(psDomain[i] == pbDomain[j]){
|
|
wolffd@0
|
60 mask[count] = j;
|
|
wolffd@0
|
61 count++;
|
|
wolffd@0
|
62 break;
|
|
wolffd@0
|
63 }
|
|
wolffd@0
|
64 }
|
|
wolffd@0
|
65 }
|
|
wolffd@0
|
66
|
|
wolffd@0
|
67 ndim = siz_b;
|
|
wolffd@0
|
68 sx = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
69 sy = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
70 for(i=0; i<ndim; i++){
|
|
wolffd@0
|
71 sx[i] = (int)bs[i];
|
|
wolffd@0
|
72 sy[i] = 1;
|
|
wolffd@0
|
73 }
|
|
wolffd@0
|
74 for(i=0; i<count; i++){
|
|
wolffd@0
|
75 temp = mask[i];
|
|
wolffd@0
|
76 sy[temp] = sx[temp];
|
|
wolffd@0
|
77 }
|
|
wolffd@0
|
78
|
|
wolffd@0
|
79 s = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
80 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
|
|
wolffd@0
|
81 subs = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
82 cpsy2 = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
83 for(i = 0; i < ndim; i++){
|
|
wolffd@0
|
84 subs[i] = 0;
|
|
wolffd@0
|
85 s[i] = sx[i] - 1;
|
|
wolffd@0
|
86 }
|
|
wolffd@0
|
87
|
|
wolffd@0
|
88 for(i = 0; i < ndim-1; i++){
|
|
wolffd@0
|
89 cpsy[i+1] = cpsy[i]*sy[i]--;
|
|
wolffd@0
|
90 cpsy2[i] = cpsy[i]*sy[i];
|
|
wolffd@0
|
91 }
|
|
wolffd@0
|
92 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
|
|
wolffd@0
|
93
|
|
wolffd@0
|
94 for(j=0; j<NB; j++){
|
|
wolffd@0
|
95 value = *sp;
|
|
wolffd@0
|
96 if(value == 0) value = 1;
|
|
wolffd@0
|
97 *zp++ /= value;
|
|
wolffd@0
|
98 for(i = 0; i < ndim; i++){
|
|
wolffd@0
|
99 if(subs[i] == s[i]){
|
|
wolffd@0
|
100 subs[i] = 0;
|
|
wolffd@0
|
101 if(sy[i])
|
|
wolffd@0
|
102 sp -= cpsy2[i];
|
|
wolffd@0
|
103 }
|
|
wolffd@0
|
104 else{
|
|
wolffd@0
|
105 subs[i]++;
|
|
wolffd@0
|
106 if(sy[i])
|
|
wolffd@0
|
107 sp += cpsy[i];
|
|
wolffd@0
|
108 break;
|
|
wolffd@0
|
109 }
|
|
wolffd@0
|
110 }
|
|
wolffd@0
|
111 }
|
|
wolffd@0
|
112 free(sx);
|
|
wolffd@0
|
113 free(sy);
|
|
wolffd@0
|
114 free(s);
|
|
wolffd@0
|
115 free(cpsy);
|
|
wolffd@0
|
116 free(subs);
|
|
wolffd@0
|
117 free(cpsy2);
|
|
wolffd@0
|
118 free(mask);
|
|
wolffd@0
|
119 }
|
|
wolffd@0
|
120
|
