comparison toolboxes/FullBNT-1.0.7/bnt/potentials/Tables/mult_by_tableC.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 /* mult_by_table.c ../potential/tables */
2
3
4 /******************************************/
5 /* 6 input & 1 output */
6 /* Big table [0] */
7 /* Big domain [1] */
8 /* big sizes [2] */
9 /* Small table [3] */
10 /* small domain [4] */
11 /* small sizes [5] */
12 /* */
13 /* New big table[0] */
14 /******************************************/
15
16 #include "mex.h"
17
18 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
19 int i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
20 int *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
21 double *pbDomain, *psDomain, *sp, *zp, *bs;
22
23 plhs[0] = mxDuplicateArray(prhs[0]);
24 zp = mxGetPr(plhs[0]);
25
26 siz_b = mxGetNumberOfElements(prhs[1]);
27 siz_s = mxGetNumberOfElements(prhs[4]);
28 pbDomain = mxGetPr(prhs[1]);
29 psDomain = mxGetPr(prhs[4]);
30
31 NB = mxGetNumberOfElements(prhs[0]);
32 NS = mxGetNumberOfElements(prhs[3]);
33 sp = mxGetPr(prhs[3]);
34
35 bs = mxGetPr(prhs[2]);
36
37 if(NS == 1){
38 for(i=0; i<NB; i++){
39 zp[i] *= *sp;
40 }
41 return;
42 }
43
44 if(NS == NB){
45 for(i=0; i<NB; i++){
46 zp[i] *= sp[i];
47 }
48 return;
49 }
50
51 mask = malloc(siz_s * sizeof(int));
52 count = 0;
53 for(i=0; i<siz_s; i++){
54 for(j=0; j<siz_b; j++){
55 if(psDomain[i] == pbDomain[j]){
56 mask[count] = j;
57 count++;
58 break;
59 }
60 }
61 }
62
63 ndim = siz_b;
64 sx = (int *)malloc(sizeof(int)*ndim);
65 sy = (int *)malloc(sizeof(int)*ndim);
66 for(i=0; i<ndim; i++){
67 sx[i] = (int)bs[i];
68 sy[i] = 1;
69 }
70 for(i=0; i<count; i++){
71 temp = mask[i];
72 sy[temp] = sx[temp];
73 }
74
75 s = (int *)malloc(sizeof(int)*ndim);
76 *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
77 subs = (int *)malloc(sizeof(int)*ndim);
78 cpsy2 = (int *)malloc(sizeof(int)*ndim);
79 for(i = 0; i < ndim; i++){
80 subs[i] = 0;
81 s[i] = sx[i] - 1;
82 }
83
84 for(i = 0; i < ndim-1; i++){
85 cpsy[i+1] = cpsy[i]*sy[i]--;
86 cpsy2[i] = cpsy[i]*sy[i];
87 }
88 cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
89
90 for(j=0; j<NB; j++){
91 *zp++ *= *sp;
92 for(i = 0; i < ndim; i++){
93 if(subs[i] == s[i]){
94 subs[i] = 0;
95 if(sy[i])
96 sp -= cpsy2[i];
97 }
98 else{
99 subs[i]++;
100 if(sy[i])
101 sp += cpsy[i];
102 break;
103 }
104 }
105 }
106 free(sx);
107 free(sy);
108 free(s);
109 free(cpsy);
110 free(subs);
111 free(cpsy2);
112 free(mask);
113 }
114