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