wolffd@0: /****************************************************
wolffd@0: A = mult_by_array(big, small)
wolffd@0: implicitely copies small |big|/|small| times 
wolffd@0: and then does element-wise multiplication.
wolffd@0: 
wolffd@0: i.e.,
wolffd@0: C = repmat(small(:), 1, length(big(:))/length(small(:)))
wolffd@0: A = reshape(big(:) .* C(:), size(big))
wolffd@0: 
wolffd@0: However, this C version avoids the expense of the repmat.
wolffd@0: 
wolffd@0: Written by wei.hu@intel.com, 28 Jan 2002.
wolffd@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: {
wolffd@0:   double		*sp, *zp;
wolffd@0:   int			i, j, NB, NS, xnd, ynd, ndim;
wolffd@0:   const int	*xdim, *ydim;
wolffd@0:   int			*s, *sx, *sy, *cpsy, *subs, *cpsy2;
wolffd@0:   
wolffd@0:   if (nrhs != 2)
wolffd@0:     mexErrMsgTxt("Incorrect number of inputs.");
wolffd@0:   
wolffd@0:   if (nlhs > 1)
wolffd@0:     mexErrMsgTxt("Too many output arguments.");
wolffd@0:   
wolffd@0:   plhs[0] = mxDuplicateArray(prhs[0]);
wolffd@0:   zp = mxGetPr(plhs[0]);
wolffd@0:   sp = mxGetPr(prhs[1]);
wolffd@0:   
wolffd@0:   xnd = mxGetNumberOfDimensions(prhs[0]);
wolffd@0:   ynd = mxGetNumberOfDimensions(prhs[1]);
wolffd@0:   xdim = mxGetDimensions(prhs[0]);
wolffd@0:   ydim = mxGetDimensions(prhs[1]);
wolffd@0:   ndim = xnd;
wolffd@0:   
wolffd@0:   NB = mxGetNumberOfElements(prhs[0]);
wolffd@0:   NS = mxGetNumberOfElements(prhs[1]);
wolffd@0:   
wolffd@0:   if(NS == 1){
wolffd@0:     for(i=0; i<NB; i++){
wolffd@0:       *zp++ *= *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++ *= *sp++;
wolffd@0:     }
wolffd@0:     return;
wolffd@0:   }
wolffd@0:   
wolffd@0:   sx = (int *)malloc(sizeof(int)*ndim);
wolffd@0:   sy = (int *)malloc(sizeof(int)*ndim);
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:     sx[i] = xdim[i];
wolffd@0:     sy[i] = (i < ynd) ? ydim[i] : 1;
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]) sp -= cpsy2[i];
wolffd@0:       }
wolffd@0:       else{
wolffd@0: 	subs[i]++;
wolffd@0: 	if(sy[i]) 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: }