wolffd@0: /* mult_by_table.c  ../potential/tables  */
wolffd@0: 
wolffd@0: 
wolffd@0: /******************************************/
wolffd@0: /* 6 input & 1 output                     */
wolffd@0: /* Big table    [0]                       */
wolffd@0: /* Big domain   [1]                       */
wolffd@0: /* big sizes    [2]                       */
wolffd@0: /* Small table  [3]                       */
wolffd@0: /* small domain [4]                       */
wolffd@0: /* small sizes  [5]                       */
wolffd@0: /*                                        */
wolffd@0: /* New big table[0]                       */
wolffd@0: /******************************************/
wolffd@0: 
wolffd@0: #include "mex.h"
wolffd@0: 
wolffd@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0:   int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp;
wolffd@0:   int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2;
wolffd@0:   double  *pbDomain, *psDomain, *sp, *zp, *bs;
wolffd@0:   
wolffd@0:   plhs[0] = mxDuplicateArray(prhs[0]);
wolffd@0:   zp = mxGetPr(plhs[0]);
wolffd@0:   
wolffd@0:   siz_b = mxGetNumberOfElements(prhs[1]);
wolffd@0:   siz_s = mxGetNumberOfElements(prhs[4]);
wolffd@0:   pbDomain = mxGetPr(prhs[1]);
wolffd@0:   psDomain = mxGetPr(prhs[4]);
wolffd@0:   
wolffd@0:   NB = mxGetNumberOfElements(prhs[0]);
wolffd@0:   NS = mxGetNumberOfElements(prhs[3]);
wolffd@0:   sp = mxGetPr(prhs[3]);
wolffd@0:   
wolffd@0:   bs = mxGetPr(prhs[2]);
wolffd@0:   
wolffd@0:   if(NS == 1){
wolffd@0:     for(i=0; i<NB; i++){
wolffd@0:       zp[i] *= *sp;
wolffd@0:     }
wolffd@0:     return;
wolffd@0:   }
wolffd@0:   
wolffd@0:   if(NS == NB){
wolffd@0:     for(i=0; i<NB; i++){
wolffd@0:       zp[i] *= sp[i];
wolffd@0:     }
wolffd@0:     return;
wolffd@0:   }
wolffd@0:   
wolffd@0:   mask = malloc(siz_s * sizeof(int));
wolffd@0:   count = 0;
wolffd@0:   for(i=0; i<siz_s; i++){
wolffd@0:     for(j=0; j<siz_b; j++){
wolffd@0:       if(psDomain[i] == pbDomain[j]){
wolffd@0: 	mask[count] = j;
wolffd@0: 	count++;
wolffd@0: 	break;
wolffd@0:       }
wolffd@0:     }
wolffd@0:   }
wolffd@0:   
wolffd@0:   ndim = siz_b;
wolffd@0:   sx = (int *)malloc(sizeof(int)*ndim);
wolffd@0:   sy = (int *)malloc(sizeof(int)*ndim);
wolffd@0:   for(i=0; i<ndim; i++){
wolffd@0:     sx[i] = (int)bs[i];
wolffd@0:     sy[i] = 1;
wolffd@0:   }
wolffd@0:   for(i=0; i<count; i++){
wolffd@0:     temp = mask[i];
wolffd@0:     sy[temp] = sx[temp];
wolffd@0:   }
wolffd@0:   
wolffd@0:   s = (int *)malloc(sizeof(int)*ndim);
wolffd@0:   *(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
wolffd@0:   subs =   (int *)malloc(sizeof(int)*ndim);
wolffd@0:   cpsy2 =  (int *)malloc(sizeof(int)*ndim);
wolffd@0:   for(i = 0; i < ndim; i++){
wolffd@0:     subs[i] = 0;
wolffd@0:     s[i] = sx[i] - 1;
wolffd@0:   }
wolffd@0:   
wolffd@0:   for(i = 0; i < ndim-1; i++){
wolffd@0:     cpsy[i+1] = cpsy[i]*sy[i]--;
wolffd@0:     cpsy2[i] = cpsy[i]*sy[i];
wolffd@0:   }
wolffd@0:   cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
wolffd@0:   
wolffd@0:   for(j=0; j<NB; j++){
wolffd@0:     *zp++ *= *sp;
wolffd@0:     for(i = 0; i < ndim; i++){
wolffd@0:       if(subs[i] == s[i]){
wolffd@0: 	subs[i] = 0;
wolffd@0: 	if(sy[i])
wolffd@0: 	  sp -= cpsy2[i];
wolffd@0:       }
wolffd@0:       else{
wolffd@0: 	subs[i]++;
wolffd@0: 	if(sy[i])
wolffd@0: 	  sp += cpsy[i];
wolffd@0: 	break;
wolffd@0:       }
wolffd@0:     }
wolffd@0:   }
wolffd@0:   free(sx);
wolffd@0:   free(sy);
wolffd@0:   free(s);
wolffd@0:   free(cpsy);
wolffd@0:   free(subs);
wolffd@0:   free(cpsy2);
wolffd@0:   free(mask);
wolffd@0: }
wolffd@0: