Matlab and Mex

On this chapter we learn how to use C/C++ (With Cuda/OpenCl) on matlab command line. Matlab has a tool called Mex that creates a module(DLL/so) from your C/C++ code that can be launched from Matlab console.

This is cool because it allows you to gather data from Matlab and exercise our functions with it.

Mex Setup

First we need to configure mex to use the compiler installed on your machine. Just type mex -setup and click on the link mex -setup C++.

Mex Hello World

Create a mex function on C that print some messages and sum the 2 parameters.

#include <math.h>
#include <matrix.h>
#include <mex.h>

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  double *a,*b,*result;

  mexPrintf("Hello World\n");
  mexPrintf("%d Outputs on left side\n",nlhs);
  mexPrintf("%d parameters on right side\n",nrhs);

  /* Delcare pointers to inputs and outputs on matlab */
  mxArray *par_mat_1, *par_mat_2, *out_mat_1;

  /* 2 parameters one output*/
  if (nrhs == 2 && nlhs == 1) {
    /*Make a deep copy of the array elements*/
    par_mat_1 = mxDuplicateArray(prhs[0]);
    par_mat_2 = mxDuplicateArray(prhs[1]);

    /* Real data elements in array of type DOUBLE */
    a = mxGetPr(par_mat_1);
    b = mxGetPr(par_mat_2);

    mexPrintf("A=%f B=%f\n",*a,*b);

    /* Create a 1x1 matrix (scalar)*/
    plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);

    /* Get pointer to output 1x1 matrix (scalar)*/
    result = mxGetPr(plhs[0]);
    *result = *a+*b;
  } else {
    /* Display error message */
    mexErrMsgIdAndTxt( "MATLAB:mexFunction:maxlhs","wrong number of arguments.");
  }
}

Compiling and loading the function on matlab.

>> mex mexHelloSum.c
Building with 'gcc'.
MEX completed successfully.
>> c = mexHelloSum(10,20)
Hello World
1 Outputs on left side
2 parameters on right side
A=10.000000 B=20.000000

c =
    30

If you observe on your file sysem a file with extension mexa64(Depend on your system) will be created. It's just a simple module (DLL on windows or SO on linux)

# Show file properties
$ file mexHelloSum.mexa64 
mexHelloSum.mexa64: ELF 64-bit LSB  shared object, x86-64, version 1 (SYSV), dynamically linked, BuildID[sha1]=889169d9ca5539b1c3bd20d77c5c17109091a54f, not stripped

# Get symbols on shared object (.so)
$ nm -g mexHelloSum.mexa64 
                 w __cxa_finalize@@GLIBC_2.2.5
                 w __gmon_start__
                 w _ITM_deregisterTMCloneTable
                 w _ITM_registerTMCloneTable
                 w _Jv_RegisterClasses
0000000000000000 A MEX
                 U mexErrMsgIdAndTxt
00000000000007e5 T mexFunction
                 U mexPrintf
                 U mxCreateDoubleMatrix_700
                 U mxDuplicateArray
                 U mxGetPr

# Only show the functions exported by this shared object (.so, .dll, etc...)
nm --defined-only -g mexHelloSum.mexa64
0000000000000000 A MEX
00000000000007e5 T mexFunction

Also "mex" is just a wrapper function that will use gcc behind the scenes. You could for instance pass compile flags to help you debug/profile. Also when you are in doubt on what mex is doing just pass the -v option. Also mex can be called from the linux console

mex -v -g CFLAGS="\$CFLAGS -std=c99 -pg -O0" LINKLIBS="\$LINKLIBS -pg" mexHelloSum.c
Building with 'gcc'.
Warning: You are using gcc version '4.8.4'. The version of gcc is not supported. The version currently supported with MEX is '4.7.x'. For a list of currently supported compilers see: http://www.mathworks.com/support/compilers/current_release.
/usr/bin/gcc -c -DMX_COMPAT_32   -D_GNU_SOURCE -DMATLAB_MEX_FILE  -I"/usr/local/MATLAB/R2016a/extern/include" -I"/usr/local/MATLAB/R2016a/simulink/include" -ansi -fexceptions -fPIC -fno-omit-frame-pointer -pthread -std=c99 -pg -O0 -g /home/leo/work/learningOpenCl/chap3/mexHelloSum.c -o /tmp/mex_478236846868418_13989/mexHelloSum.o
/usr/bin/gcc -pthread -Wl,--no-undefined -Wl,-rpath-link,/usr/local/MATLAB/R2016a/bin/glnxa64 -shared  -g -Wl,--version-script,"/usr/local/MATLAB/R2016a/extern/lib/glnxa64/mexFunction.map" /tmp/mex_478236846868418_13989/mexHelloSum.o   -L"/usr/local/MATLAB/R2016a/bin/glnxa64" -lmx -lmex -lmat -lm -lstdc++ -pg -o mexHelloSum.mexa64
rm -f /tmp/mex_478236846868418_13989/mexHelloSum.o

Debugging with Mex

GDB on console

It's possible to use gdb to debug mex functions if you load matlab through gdb. This is not the coolest way (using console) but it show the steps necessary to do the same on eclipse.

1) Run matlab with the -Dgdb option this will open gdb

  matlab -Dgdb

2) From gdb

handle SIGSEGV SIGBUS nostop noprint
run -nojvm

3) From matlab enable mex debugging and call your function

dbmex on
ls
mexHelloSum.c  mexHelloSum.mexa64  texput.log

c = mexHelloSum(10,20)

MEX FILE: /home/leo/work/learningOpenCl/chap3/mexHelloSum.mexa64 entry point located at address 0xbcb18815
Add breakpoints at the debugger prompt and issue a "continue" to resume 
execution of MATLAB.

Program received signal SIGUSR1, User defined signal 1.

4) Now when you execute the function it will switch back to gdb where you can break at your function

(gdb) break mexFunction
Breakpoint 1 at 0x7fffbcb18831: file /home/leo/work/learningOpenCl/chap3/mexHelloSum.c, line 9.
(gdb) cont
Continuing.
[Switching to Thread 0x7fffd6c6c700 (LWP 14216)]

Breakpoint 1, mexFunction (nlhs=1, plhs=0x7fffd6c69588, nrhs=2, prhs=0x7fffd6c694a8) at /home/leo/work/learningOpenCl/chap3/mexHelloSum.c:9
9      mexPrintf("Hello World\n");
(gdb) next
Hello World
10      mexPrintf("%d Outputs on left side\n",nlhs);
(gdb) list
5    void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
6    {
7      double *a,*b,*result;
8    
9      mexPrintf("Hello World\n");
10      mexPrintf("%d Outputs on left side\n",nlhs);
11      mexPrintf("%d parameters on right side\n",nrhs);
12    
13      /* Delcare pointers to inputs and outputs on matlab */
14      mxArray *par_mat_1, *par_mat_2, *out_mat_1;

More realistic use case

Here we're creating a mex function that will serve as a bridge for some matrix 2d implementation. What you can notice here is that we query the input parameters dimensions with the functions mxGetN and mxGetM. Also we check if the inputs has a particular type (float/single) with mxIsSingle and checking if the input is not a complex number with mxIsComplex

/* 
 * mex_matMult2D.c - Matlab bridge for various matrix 2d implementation.
 * This function will allow you to call directly from matlab some matrix
 * 2d implementation
 *
 * Multiplies 2 matrices A[MxN], B[NxO], resulting on C[MxO]
 *
 * The calling syntax is:
 *
 *        C = mex_matMult2D(A, B)
 *
 * This is a MEX file for MATLAB.
*/

/* Mex headers */
#include "mex.h"

/* Define contract (Function prototype for matrix multiplication) */
#include "matrix_2d_mul.h"

/* The gateway function */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
/*verify input and output parameters*/
if(nrhs != 2) {
    mexErrMsgIdAndTxt("matMul2D:nrhs",
                      "Two inputs required.");
}

if(nlhs != 1) {
    mexErrMsgIdAndTxt("matMul2D:nlhs",
                      "One output required.");
}

/* make sure the second input argument is type single */
if( !mxIsSingle(prhs[0]) || !mxIsSingle(prhs[1])
     || mxIsComplex(prhs[0]) || mxIsComplex(prhs[1])) {
    mexErrMsgIdAndTxt("matMul2D:notSingle","Input matrix must be type single.");
}

/* check input matrix is 2D, maximal*/
if( ((mxGetNumberOfDimensions(prhs[0])!=1) && (mxGetNumberOfDimensions(prhs[0])!=2)) || 
    ((mxGetNumberOfDimensions(prhs[1])!=1) && (mxGetNumberOfDimensions(prhs[1])!=2)) ) {
    mexErrMsgIdAndTxt("matMul2D:moreThan2D","Max dimension of input matrix is 2.");
}

/* ncolsA should equal to nrowsB*/
if (mxGetN(prhs[0])!=mxGetM(prhs[1])) {
    mexErrMsgIdAndTxt("matMul2D:mismatchDimension",
                    "For A*B, col num of A should equal to row num of B.");
}

/* variable declarations here */
float * inMatrixA;
float * inMatrixB;
float * outMatrix;   

mwSize ncolsA;
mwSize nrowsA;
mwSize ncolsB;
mwSize nrowsB;

/* Get pointer for data (input parameters) in matlab */
inMatrixA = (float *)mxGetPr(prhs[0]);
inMatrixB = (float *)mxGetPr(prhs[1]);

/* Get matrices dimensions */
ncolsA = mxGetN(prhs[0]);
nrowsA = mxGetM(prhs[0]);
ncolsB = mxGetN(prhs[1]);
nrowsB = mxGetM(prhs[1]);

/* creat ouput matrix */ 
plhs[0] = mxCreateNumericMatrix(nrowsA,ncolsB,mxSINGLE_CLASS,mxREAL);
outMatrix = (float *)mxGetPr(plhs[0]);


/* 
 * Call some implementation of matrix_2d_mul_float. The idea is that you
 * could swap implementations if they follow same contract. Here contract
 * just means same function name and parameters
 */
matrix_2d_mul_float(inMatrixA, inMatrixB, outMatrix,
                    nrowsA, ncolsA,
                    nrowsB, ncolsB,
                    0);
}

This function is wrapping this naive matrix 2d implementation.

void matrix_2d_mul_float(float *A, float *B, float *C, int num_rows_A, int num_cols_A, int num_rows_B, int num_cols_B ,int rowColMajor) {
  float sum = 0;
  int num_rows_C = num_rows_A;
  int num_cols_C = num_cols_B;
  // Iterate on each row of A
  for(int i=0; i<num_rows_A; i++) {
    // Iterate on each collumn of B
    for (int k=0; k<num_cols_B; k++) {
      sum = 0;
      // Do the "multiply add between" row of A and collumn of B
      for (int j=0; j<num_cols_A; j++){
        //sum += A[i][j]*B[j][k];
        if (rowColMajor) {
          sum += A[i*num_cols_A+j]*B[j*num_cols_B+k];
        } else {
          sum += A[j*num_rows_A+i]*B[k*num_rows_B+j];
        }
      }
      // C[i][k] == C[i*num_cols_C+k]
      if (rowColMajor) {
        C[i*num_cols_B+k]=sum;
      } else {
        C[k*num_rows_A+i]=sum;
      }
    }
  }
}