|
wolffd@0
|
1 /* mult_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;
|
|
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 for(i=0; i<NB; i++){
|
|
wolffd@0
|
39 zp[i] *= *sp;
|
|
wolffd@0
|
40 }
|
|
wolffd@0
|
41 return;
|
|
wolffd@0
|
42 }
|
|
wolffd@0
|
43
|
|
wolffd@0
|
44 if(NS == NB){
|
|
wolffd@0
|
45 for(i=0; i<NB; i++){
|
|
wolffd@0
|
46 zp[i] *= sp[i];
|
|
wolffd@0
|
47 }
|
|
wolffd@0
|
48 return;
|
|
wolffd@0
|
49 }
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 mask = malloc(siz_s * sizeof(int));
|
|
wolffd@0
|
52 count = 0;
|
|
wolffd@0
|
53 for(i=0; i<siz_s; i++){
|
|
wolffd@0
|
54 for(j=0; j<siz_b; j++){
|
|
wolffd@0
|
55 if(psDomain[i] == pbDomain[j]){
|
|
wolffd@0
|
56 mask[count] = j;
|
|
wolffd@0
|
57 count++;
|
|
wolffd@0
|
58 break;
|
|
wolffd@0
|
59 }
|
|
wolffd@0
|
60 }
|
|
wolffd@0
|
61 }
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 ndim = siz_b;
|
|
wolffd@0
|
64 sx = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
65 sy = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
66 for(i=0; i<ndim; i++){
|
|
wolffd@0
|
67 sx[i] = (int)bs[i];
|
|
wolffd@0
|
68 sy[i] = 1;
|
|
wolffd@0
|
69 }
|
|
wolffd@0
|
70 for(i=0; i<count; i++){
|
|
wolffd@0
|
71 temp = mask[i];
|
|
wolffd@0
|
72 sy[temp] = sx[temp];
|
|
wolffd@0
|
73 }
|
|
wolffd@0
|
74
|
|
wolffd@0
|
75 s = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
76 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
|
|
wolffd@0
|
77 subs = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
78 cpsy2 = (int *)malloc(sizeof(int)*ndim);
|
|
wolffd@0
|
79 for(i = 0; i < ndim; i++){
|
|
wolffd@0
|
80 subs[i] = 0;
|
|
wolffd@0
|
81 s[i] = sx[i] - 1;
|
|
wolffd@0
|
82 }
|
|
wolffd@0
|
83
|
|
wolffd@0
|
84 for(i = 0; i < ndim-1; i++){
|
|
wolffd@0
|
85 cpsy[i+1] = cpsy[i]*sy[i]--;
|
|
wolffd@0
|
86 cpsy2[i] = cpsy[i]*sy[i];
|
|
wolffd@0
|
87 }
|
|
wolffd@0
|
88 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
|
|
wolffd@0
|
89
|
|
wolffd@0
|
90 for(j=0; j<NB; j++){
|
|
wolffd@0
|
91 *zp++ *= *sp;
|
|
wolffd@0
|
92 for(i = 0; i < ndim; i++){
|
|
wolffd@0
|
93 if(subs[i] == s[i]){
|
|
wolffd@0
|
94 subs[i] = 0;
|
|
wolffd@0
|
95 if(sy[i])
|
|
wolffd@0
|
96 sp -= cpsy2[i];
|
|
wolffd@0
|
97 }
|
|
wolffd@0
|
98 else{
|
|
wolffd@0
|
99 subs[i]++;
|
|
wolffd@0
|
100 if(sy[i])
|
|
wolffd@0
|
101 sp += cpsy[i];
|
|
wolffd@0
|
102 break;
|
|
wolffd@0
|
103 }
|
|
wolffd@0
|
104 }
|
|
wolffd@0
|
105 }
|
|
wolffd@0
|
106 free(sx);
|
|
wolffd@0
|
107 free(sy);
|
|
wolffd@0
|
108 free(s);
|
|
wolffd@0
|
109 free(cpsy);
|
|
wolffd@0
|
110 free(subs);
|
|
wolffd@0
|
111 free(cpsy2);
|
|
wolffd@0
|
112 free(mask);
|
|
wolffd@0
|
113 }
|
|
wolffd@0
|
114
|
