Daniel@0: /* C mex version of max_mult.m in BPMRF2 directory  */
Daniel@0: /* gcc -Wall -I/mit/matlab_v6.5/distrib/bin/glnx86 -c max_mult.c */
Daniel@0: 
Daniel@0: #include <math.h>
Daniel@0: #include "mex.h"
Daniel@0: 
Daniel@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
Daniel@0: {
Daniel@0:   int rows,cols,common,m,n,p;
Daniel@0:   double y1, y2;
Daniel@0:   double *arr1, *arr2, *arr3;
Daniel@0: 
Daniel@0: 
Daniel@0:   if (nrhs!=2 || nlhs>1)
Daniel@0:     mexErrMsgTxt("max_mult requires two inputs and one output");
Daniel@0:   if (mxIsChar(prhs[0]) || mxIsClass(prhs[0], "sparse") || mxIsComplex(prhs[0])
Daniel@0:       || mxIsChar(prhs[1]) || mxIsClass(prhs[1], "sparse") || mxIsComplex(prhs[1]))
Daniel@0:     mexErrMsgTxt("Inputs must be real, full, and nonstring");
Daniel@0:   if (mxGetN(prhs[0])!=mxGetM(prhs[1]))
Daniel@0:     mexErrMsgTxt("The number of columns of A must be the same as the number of rows of x");
Daniel@0:   
Daniel@0: 
Daniel@0:   arr1=mxGetPr(prhs[0]);
Daniel@0:   arr2=mxGetPr(prhs[1]);
Daniel@0:   p=mxGetN(prhs[0]);
Daniel@0:   m=mxGetM(prhs[0]);
Daniel@0:   n=mxGetN(prhs[1]);
Daniel@0:   plhs[0]=mxCreateDoubleMatrix(m, n, mxREAL);
Daniel@0:   arr3=mxMalloc(m*n*sizeof(double));
Daniel@0: 
Daniel@0:   for (rows=0; rows<m ; rows++)
Daniel@0:     for (cols=0; cols<n ; cols++)
Daniel@0:     {
Daniel@0:       y1=arr1[rows]*arr2[cols*p];  
Daniel@0:       for (common=1; common<p; common++)
Daniel@0:       {
Daniel@0: 	y2=arr1[rows+common*m]*arr2[common+cols*p];
Daniel@0:         if (y2>y1)
Daniel@0:           y1=y2;
Daniel@0:       }
Daniel@0:       arr3[rows+cols*m]=y1;
Daniel@0:     }
Daniel@0: 
Daniel@0:   mxSetPr(plhs[0], arr3);
Daniel@0: 	 
Daniel@0: }