Combining GSL and MEX

NB: I assume you already have basic knowledge of both GSL and MEX files. If this is not the case, you can see my previous posts about GSL and about MEX files.

Using GSL in a MEX file is actually pretty straightforward, as MEX files can contain any C/C++ code. It can be useful for the following cases:

• You have some existing C/C++ code using GSL, and you want to use it in Matlab.
• You have to write a Matlab operation as a C/C++ MEX file, and it would be easier with GSL.

Conversion functions

Here are some short functions I use to convert MEX vectors/matrices to their GSL equivalents. You can just copy-paste them where you need, but if you want to use them in a multi-file project, you can get the code and header files in my repository.

Note that they do deep copies, entry after entry, and leave the original matrices untouched. I never had performance problems because I usually convert the input matrices only once, and then do heavy operations in the MEX files. However, if you perform lots of small operations on big matrices, this may not be the best solution. If you don’t need deep copies, you can just copy the data pointers using mexGetPr and gsl_matrix_ptr (see here).

Convert a mxArray to a gsl_matrix

gsl_matrix *
mxmat_to_gslmat (const mxArray * mxmat)
{
  const mwSize *dims = mxGetDimensions (mxmat);
  int nb_row = (int) dims[0];
  int nb_col = (int) dims[1];
  double *rawmat = mxGetPr (mxmat);

  gsl_matrix *m = gsl_matrix_alloc (nb_row, nb_col);
  for (int i = 0; i < nb_row; ++i)
    for (int j = 0; j < nb_col; ++j)
      gsl_matrix_set (m, i, j, rawmat[i + j * nb_row]);

  return m;
}

Convert a mxArray to a gsl_vector

gsl_vector *
mxmat_to_gslvec (const mxArray * mxmat)
{
  const mwSize *dims = mxGetDimensions (mxmat);
  int size = (int) dims[0];
  double *rawvec = mxGetPr (mxmat);

  gsl_vector *m = gsl_vector_alloc (size);
  for (int i = 0; i < size; ++i)
      gsl_vector_set (m, i, rawvec[i]);

  return m;
}

Convert a gsl_matrix to a mxArray

mxArray *
gslmat_to_mxmat (const gsl_matrix * gslmat)
{
  int nb_row = gslmat->size1;
  int nb_col = gslmat->size2;

  mxArray *m = mxCreateDoubleMatrix (nb_row, nb_col, mxREAL);
  double *pm = mxGetPr (m);
  for (int i = 0; i < nb_row; ++i)
    for (int j = 0; j < nb_col; ++j)
      pm[i + j * nb_row] = gsl_matrix_get (gslmat, i, j);

  return m;
}

Convert a gsl_vector to a mxArray

mxArray *
gslvec_to_mxmat (const gsl_vector * gslvec)
{
  int size = gslvec->size;

  mxArray *m = mxCreateDoubleMatrix (size, 1, mxREAL);
  double *pm = mxGetPr (m);
  for (int i = 0; i < size; ++i)
      pm[i] = gsl_vector_get (gslvec, i);

  return m;
}

Examples of use

You can see a one-file example here: acc-hals-mex.

You can also have, for the same code, a MEX interface and a standard C interface. You just have to define one main and one mexFunction that call the same code. See for example my code for sparse-nmf.

Notes

Note that Matlab uses mxArray to store both matrices and vectors. Vectors are considered normal matrices, but with only one column (or one row).

In GSL, vectors are distinct objects. By default, they are column vectors. You can use them for row vectors, but you have to remember they are stored in column before doing operations with them.