rmeddis@0: /* intpow - Nick Clark - 3/5/2011
rmeddis@0:  *
rmeddis@0:  * Function that raises a double precision vector to an integer power 
rmeddis@0:  * vector. The following call to this mex function... 
rmeddis@0:  *
rmeddis@0:  * >>      z = intpow(x,y)
rmeddis@0:  * 
rmeddis@0:  * ... does the equivalent of the following in Matlab...
rmeddis@0:  *
rmeddis@0:  * >>      z = x .^ ceil(y)
rmeddis@0:  * 
rmeddis@0:  * NOTE: Under most circumstances, this function is slower than MATLAB's 
rmeddis@0:  * built in ^ operator, but for the (generally) small integer powers used 
rmeddis@0:  * in MAP we observe massive performance boosts using this C function.
rmeddis@0:  */
rmeddis@0: 
rmeddis@0: /*************************************************************************/
rmeddis@0: /* Header(s)                                                             */
rmeddis@0: /*************************************************************************/
rmeddis@0: #include "mex.h"
rmeddis@0: 
rmeddis@0: /*************************************************************************/
rmeddis@0: /* Input vars                                                            */
rmeddis@0: /*************************************************************************/
rmeddis@0: #define IN_x 	prhs[0]
rmeddis@0: #define IN_y 	prhs[1]
rmeddis@0: 
rmeddis@0: /*************************************************************************/
rmeddis@0: /* Output vars                                                           */
rmeddis@0: /*************************************************************************/
rmeddis@0: #define OUT_z 	plhs[0]
rmeddis@0: 
rmeddis@0: /*************************************************************************/
rmeddis@0: /* Gateway function and error checking                                   */
rmeddis@0: /*************************************************************************/
rmeddis@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
rmeddis@0: {
rmeddis@0:     /* variable declarations */
rmeddis@0:     int numElements, xM, xN, yM, yN; 
rmeddis@0:     int nn, kk; 
rmeddis@0:     double *x, *y, *z;
rmeddis@0:      
rmeddis@0:     /*  check the number of input and output parameters  */  
rmeddis@0:     if(nrhs!=2)
rmeddis@0:         mexErrMsgTxt("intpow : Two input args expected");
rmeddis@0:     if(nlhs > 1)
rmeddis@0:         mexErrMsgTxt("intpow : Too many outputs");
rmeddis@0:     
rmeddis@0:     /* check that x and y have equal size */ 
rmeddis@0:     x = mxGetPr(IN_x);
rmeddis@0:     xM = mxGetM(IN_x);
rmeddis@0:     xN = mxGetN(IN_x);
rmeddis@0:     
rmeddis@0:     y = mxGetPr(IN_y);
rmeddis@0:     yM = mxGetM(IN_y);
rmeddis@0:     yN = mxGetN(IN_y);
rmeddis@0:     
rmeddis@0:     if  (xM != yM || xN != yN)
rmeddis@0:         mexErrMsgTxt("intpow : x and y must have equal size");
rmeddis@0:     
rmeddis@0:     /* find upper loop boundary */
rmeddis@0:     numElements = xM * xN; 
rmeddis@0:     
rmeddis@0:     /*  allocate memory and pointer for the output array */ 
rmeddis@0:     OUT_z = mxCreateDoubleMatrix(xM,xN,mxREAL);
rmeddis@0:     z = mxGetPr(OUT_z);
rmeddis@0:     
rmeddis@0:     /*  do stuff */ 
rmeddis@0:     for( nn = 0;  nn<numElements; ++nn )
rmeddis@0:     {
rmeddis@0:         z[nn] = 1.0;
rmeddis@0:         for (kk = 0; kk<y[nn]; ++kk)
rmeddis@0:         {
rmeddis@0:             z[nn] = z[nn]*x[nn];
rmeddis@0:         }
rmeddis@0:     }            
rmeddis@0: }