#include "matrix.h"
void mat_multiply(double *A, double *B, double *C,
   int mA, int nA, int mB, int nB)
{
   int i,j,k,m=0;
   for (i=0; i<mA; i++){
      for (j=0; j<nB; j++){
         C[j*mA+i]=0;
         for (k=0; k<mB; k++){
            C[j*mA+i]+=A[k*mA+i]*B[j*mB+k];
      }}}
}
/* Main interface to MATLAB */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
   double *Ap, *Bp, *Cp;
   int mA,nA,mB,nB,mC,nC;
   Ap=mxGetPr(prhs[0]); Bp=mxGetPr(prhs[1]); 
   mA=mxGetM(prhs[0]); nA=mxGetN(prhs[0]); 
   mB=mxGetM(prhs[1]); nB=mxGetN(prhs[1]);
   plhs[0]=mxCreateDoubleMatrix(mA,nB,mxREAL); 
   Cp=mxGetPr(plhs[0]);
   mat_multiply(Ap, Bp, Cp, mA, nA, mB, nB);
}