matthiasm@8: /* mult_by_sparse_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 <math.h>
matthiasm@8: #include <stdlib.h>
matthiasm@8: #include "mex.h"
matthiasm@8: 
matthiasm@8: int compare(const void* src1, const void* src2){
matthiasm@8: 	int i1 = *(int*)src1 ;
matthiasm@8: 	int i2 = *(int*)src2 ;
matthiasm@8: 	return i1-i2 ;
matthiasm@8: }
matthiasm@8: 
matthiasm@8: void ind_subv(int index, const int *cumprod, int n, int *bsubv){
matthiasm@8: 	int i;
matthiasm@8: 
matthiasm@8: 	for (i = n-1; i >= 0; i--) {
matthiasm@8: 		bsubv[i] = ((int)floor(index / cumprod[i]));
matthiasm@8: 		index = index % cumprod[i];
matthiasm@8: 	}
matthiasm@8: }
matthiasm@8: 
matthiasm@8: int subv_ind(const int n, const int *cumprod, const int *subv){
matthiasm@8: 	int i, index=0;
matthiasm@8: 
matthiasm@8: 	for(i=0; i<n; i++){
matthiasm@8: 		index += subv[i] * cumprod[i];
matthiasm@8: 	}
matthiasm@8: 	return index;
matthiasm@8: }
matthiasm@8: 
matthiasm@8: void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
matthiasm@8: 	double *ptr;
matthiasm@8: 	void   *newptr;
matthiasm@8: 	int    *ir, *jc;
matthiasm@8: 	int    nbytes;
matthiasm@8: 
matthiasm@8: 	if(new_nzmax == old_nzmax) return;
matthiasm@8: 	nbytes = new_nzmax * sizeof(*ptr);
matthiasm@8: 	ptr = mxGetPr(spArray);
matthiasm@8: 	newptr = mxRealloc(ptr, nbytes);
matthiasm@8: 	mxSetPr(spArray, newptr);
matthiasm@8: 	nbytes = new_nzmax * sizeof(*ir);
matthiasm@8: 	ir = mxGetIr(spArray);
matthiasm@8: 	newptr = mxRealloc(ir, nbytes);
matthiasm@8: 	mxSetIr(spArray, newptr);
matthiasm@8: 	jc = mxGetJc(spArray);
matthiasm@8: 	jc[0] = 0;
matthiasm@8: 	jc[1] = new_nzmax;
matthiasm@8: 	mxSetNzmax(spArray, new_nzmax);
matthiasm@8: }
matthiasm@8: 
matthiasm@8: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
matthiasm@8: 	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
matthiasm@8: 	int     *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
matthiasm@8: 	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr;
matthiasm@8: 
matthiasm@8: 	pbDomain = mxGetPr(prhs[1]);
matthiasm@8: 	bdim = mxGetNumberOfElements(prhs[1]);
matthiasm@8: 	psDomain = mxGetPr(prhs[4]);
matthiasm@8: 	sdim = mxGetNumberOfElements(prhs[4]);
matthiasm@8: 
matthiasm@8: 	pbSize = mxGetPr(prhs[2]);
matthiasm@8: 	psSize = mxGetPr(prhs[5]);
matthiasm@8: 
matthiasm@8: 	NB = 1;
matthiasm@8: 	for(i=0; i<bdim; i++){
matthiasm@8: 		NB *= (int)pbSize[i];
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	bpr = mxGetPr(prhs[0]);
matthiasm@8: 	bir = mxGetIr(prhs[0]);
matthiasm@8: 	bjc = mxGetJc(prhs[0]);
matthiasm@8: 	NZB = bjc[1];
matthiasm@8: 
matthiasm@8: 	spr = mxGetPr(prhs[3]);
matthiasm@8: 	sir = mxGetIr(prhs[3]);
matthiasm@8: 	sjc = mxGetJc(prhs[3]);
matthiasm@8: 	NZS = sjc[1];
matthiasm@8: 
matthiasm@8: 	plhs[0] = mxDuplicateArray(prhs[0]);
matthiasm@8: 	rpr = mxGetPr(plhs[0]);
matthiasm@8: 	rir = mxGetIr(plhs[0]);
matthiasm@8: 	rjc = mxGetJc(plhs[0]);
matthiasm@8: 	rjc[0] = 0;
matthiasm@8: 	rjc[1] = NZB;
matthiasm@8: 
matthiasm@8: 	if(sdim == 0){
matthiasm@8: 		for(i=0; i<NZB; i++){
matthiasm@8: 			rpr[i] *= *spr;
matthiasm@8: 		}	
matthiasm@8: 		return;
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	mask = malloc(sdim * sizeof(int));
matthiasm@8: 	bCumprod = malloc(bdim * sizeof(int));
matthiasm@8: 	sCumprod = malloc(sdim * sizeof(int));
matthiasm@8: 	bsubv = malloc(bdim * sizeof(int));
matthiasm@8: 	ssubv = malloc(sdim * sizeof(int));
matthiasm@8: 
matthiasm@8: 	count = 0;
matthiasm@8: 	for(i=0; i<sdim; i++){
matthiasm@8: 		for(j=0; j<bdim; 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: 	bCumprod[0] = 1;
matthiasm@8: 	for(i=0; i<bdim-1; i++){
matthiasm@8: 		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
matthiasm@8: 	}
matthiasm@8: 	sCumprod[0] = 1;
matthiasm@8: 	for(i=0; i<sdim-1; i++){
matthiasm@8: 		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	for(i=0; i<NZB; i++){
matthiasm@8: 		bindex = bir[i];
matthiasm@8: 		ind_subv(bindex, bCumprod, bdim, bsubv);
matthiasm@8: 		for(j=0; j<sdim; j++){
matthiasm@8: 			ssubv[j] = bsubv[mask[j]];
matthiasm@8: 		}
matthiasm@8: 		sindex = subv_ind(sdim, sCumprod, ssubv);
matthiasm@8: 		result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
matthiasm@8: 		if(result){
matthiasm@8: 			position = result - sir;
matthiasm@8: 			rpr[nzCounts] = bpr[i] * spr[position];
matthiasm@8: 			rir[nzCounts] = bindex;
matthiasm@8: 			nzCounts++;
matthiasm@8: 		}
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	reset_nzmax(plhs[0], NZB, nzCounts);
matthiasm@8: 	free(mask);
matthiasm@8: 	free(bCumprod);
matthiasm@8: 	free(sCumprod);
matthiasm@8: 	free(bsubv);
matthiasm@8: 	free(ssubv);
matthiasm@8: }