How Do I use Mexcallmatlab to Call a User-defined function?

I am trying to parallelize a section of my Matlab code using OpenMP in a mex file. The section in theMatlab code that I want to parallelize is:

   for i = 1 : n

        D(:, i) = CALC(A, B(:,i), C(i));

    end

I have written this in order to parallelize it:

 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
    size_t r,n,i,G;
    double *A, *B, *C, *D;
    int nthreads;
    nthreads = 4;
    A = mxGetPr(prhs[0]); /* first input matrix */
    B = mxGetPr(prhs[1]); /* second input matrix */
    C = mxGetPr(prhs[2]);/* third input matrix */
    /* dimensions of input matrices */
    r = mxGetN(prhs[0]);
    n = mxGetN(prhs[1]);
    plhs[0] = mxCreateDoubleMatrix(r,n, mxREAL);
    D = mxGetPr(plhs[0]);
    G=n/nthreads;
    omp_set_num_threads(nthreads);
    #pragma omp parallel for schedule (dynamic, G)
    {
      for i = 1 : n
            D(:, i) = CALC(A, B(:,i), C(i));
    }
 }

CALC is a Matlab function I have written. My challenge is how to use Mexcallmatlab to call in the CALC function to the mex file so that it can execute it in parallel inside my mex file, and return the elements of each column of D (i.e. D(:, i) back to my Matlab code.

Sorry for the lenghty question. Any help I can get on this will be highly appreciated.

Answer 1

[James Tursa]

I doubt you will be able to do this. Inside of mex routines, the API functions that deal with memory management are not (or at least they have not been in the past) thread safe. E.g., all of the mxCreateETC functions dynamically allocate memory for a new mxArray, so those functions are not thread safe. And the memory allocation functions mxMalloc, mxCalloc, and mxRealloc are not thread safe because they allocate memory. mxGetProperty also allocates memory and is not thread safe. Several other functions are in this category as well. So unless things have fundamentally changed in the API function libraries recently, you can't do anything that involves the MATLAB Memory Manager in parallel threads safely.

Functions that are thread safe would be functions like mxGetPr, mxGetPi, mxIsDouble, etc ... i.e., functions that simply return a value or status and do not allocate memory. Those could be safely used inside parallel threads.

Unfortunately, mexCallMATLAB allocates memory for the return value in the plhs array, so mexCallMATLAB is not thread safe in general. Even if mexCallMATLAB did not return anything (e.g., nlhs = 0), it would not surprise me that the memory allocations going on in the MATLAB function you called (i.e., CALC) were not thread safe.

My best advice would be to convert CALC into C code and then put that inside your OpenMP threads. If there is memory allocation involved, do it outside the parallel code.

You could of course ignore all of the above and simply code your CALC call inside the loop and see what happens, but I suspect it will crash MATLAB at some point. Consider this line:

D(:, i) = CALC(A, B(:,i), C(i));

The bulk of the C code for the above line would be getting the mxArray's for B(:,i) and C(i), and then copying the result from mexCallMATLAB into D(:,i). An outline of of how you would do this is as follows:

#include <string.h>
size_t m;
double *d, *b;
mxArray *Ci, *Bi;
mxArray rhs[3], lhs[1];
    :
m = mxGetM(prhs[1]);
    :
for( i=0; i<n; i++ ) {
    :
    Ci = mxCreateDoubleScalar( *C++ );
    Bi = mxCreateDoubleMatrix( m, 1, mxREAL );
    b = mxGetPr(Bi);
    memcpy( b, B, m * sizeof(*b) );
    B += m;
    rhs[0] = prhs[0];
    rhs[1] = Bi;
    rhs[2] = Ci;
    mexCallMATLAB( 1, lhs, 3, rhs, "CALC" );
    mxDestroyArray(Bi);
    mxDestroyArray(Ci);
    d = mxGetPr(lhs[0]);
    memcpy( D, d, r * sizeof(*d) );
    D += r;
    mxDestroyArray(lhs[0]);
}

[Chimalume Atuchukwu]

Hi John,

Thank you very much for your answer and the explanations you gave (I've read a lot of your answers in this forum and I really hoped you would come across my question and assist me with an answer). I was initially hesitant to write the CALC function in C because there's a line in it that calculates pinv(). Please, can you explain to me how I can calculate pseudoinverse in C? The particular line in the CALC function is:

    x = pinv(Z)*y;

where Z is an m by n matrix, y is an n by 1 vector, and x is an m by 1 vector. Many thanks in advance for your help.

[James Tursa]

What are the dimensions involved? Are you up for linking with BLAS and LAPACK libraries inside your C code?

As an alternative, there is an FEX submission called MMX by Yuval that has nD backslash operations multi-threaded (I think) as a mex routine you can call from m-files. So if you can formulate your input arrays to conform to what this expects maybe you can parallelize this computation using MMX called from the m-file level and avoid your own mex routine completely.

http://www.mathworks.com/matlabcentral/fileexchange/37515-mmx-multithreaded-matrix-operations-on-n-d-matrices

Without seeing the CALC code it is not possible for me to advise you on whether MMX will be useful to you or not.

[Chimalume Atuchukwu]

Hi James, the CALC function is outlined below:

function [x] = CALC(Z, y, s)

r = y;

 index_cols = []; 
 atoms      = [];  

   for i = 1 : s

    [max_r, lambda_t] = max(abs(r'*Z));
    index_cols = [index_cols, lambda_t];
    atoms      = [atoms, Z(:,lambda_t)]; 
    x_t = pinv(atoms)*y;
    r = y - atoms*x_t;
  end
x = zeros(n,1);
x(index_cols) = x_t;

end

Answer 2

[Philip Borghesani]

This will not work. MATLAB and mexCallMatlab are not thread safe, there is no way to call a MATLAB function from an OpenMP for loop and expect proper operation. I suggest investigating the parfor MATLAB loop.

Answer 3

[Arwel]

One way I have done this in the past is rather than trying to call back into the origonal Matlab workspace, instead to launch seperate matlab engine sessions for each thread in the loop and use them instead. They seem to stay separate ok (but I havent tested that to destruction!)

So, within the OMP loop, instead call a function such as the one below. This is from an old example I used and the inputs will differ for yours, but it starts a Matlab engine, and then excecutes a function named in the string myFun using 'eval'. In this case the function has four inputs and two outputs - you just need to put them as correct variable names into the engine workspace first.

When you put this in the OMP loop, what happens is that multiple engines open, and it seemed to me to work correctly with no clashes between instances. I had expected to need to do more work to get that, but it just seemed to do it okay. So for four cores, 4 sessions open and so on. They are slow to start on the first pass through, but then seem to stay open until destroyed.

It should be possible to adapt something for this for your purposes, by just making the 'evalstring' what you need. ( Caveat: I haven't tried this in a while, but from memory worked okay!)

Cheers,

Arwel

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "engine.h"


void matlabCallFun(double params[], int arrayLen, char *funName, char *pathCall, double bulkIn, double bulkOut, double contrast, double *sum, double *ans)
{
	static Engine *ep;
	static double engStatus = 0;
	mxArray *result = NULL;
	mxArray *PARAMS = NULL;
	mxArray *BULKIN = NULL;
	mxArray *BULKOUT = NULL;
	mxArray *CONTRAST = NULL;
	mxArray *lays = NULL;

	double* s;
	double* layers;


	/*Contrast = -1 means start engine only, then return*/
	if(contrast == -1) {
		ep = engOpen("");
				if(ep==0) {
					printf("Connecton to Matlab Engine failed\n");
					return;
				}
				else {
					printf("Starting Matlab Engine\n");
					engStatus == 1;
					return;
				}
	}

    /*Contrast = -2 means close engine only (if open), then return*/
	if(contrast == -2) {
				if(ep==0) {
					printf("No engine open\n");
				}
				else {
					printf("Closing Matlab Engine\n");
					engClose(ep);
					engStatus == 0;
					return;
				}
	}


	/*Automatically start engine if closed to catch non-initialised call*/
	if(engStatus == 0) {
		ep = engOpen("");		
		if(ep==0) {
			printf("Connecton to Matlab Engine failed\n");
			return;
		}
		else {
			printf("Starting Matlab Engine\n");
			mxArray* MYPATH = mxCreateString(pathCall);
			engPutVariable(ep, "mypath", MYPATH);
			engEvalString(ep,"eval(mypath)");
			mxDestroyArray(MYPATH);
			engStatus = 1;
		}
	}

	PARAMS = mxCreateDoubleMatrix(1,arrayLen,mxREAL);
	memcpy(mxGetPr(PARAMS), params, arrayLen*sizeof(double));
	engPutVariable(ep,"params",PARAMS);

	BULKIN = mxCreateDoubleMatrix(1,1,mxREAL);
	memcpy((void *)mxGetPr(BULKIN), &bulkIn, 1*sizeof(double));
	engPutVariable(ep,"bulk_in",BULKIN);

	BULKOUT = mxCreateDoubleMatrix(1,1,mxREAL);
	memcpy((void *)mxGetPr(BULKOUT), &bulkOut, 1*sizeof(double));
	engPutVariable(ep,"bulk_out",BULKOUT);

	CONTRAST = mxCreateDoubleMatrix(1,1,mxREAL);
	memcpy((void *)mxGetPr(CONTRAST), &contrast, 1*sizeof(double));
	engPutVariable(ep,"contrast",CONTRAST);

	mxArray* MYFUN = mxCreateString(funName);
	engPutVariable(ep, "myfun", MYFUN);
	engEvalString(ep,"eval(myfun)");
    
    /*Or a full sting instead of eval e.g. */ 
	/*engEvalString(ep,"[total,layers] = debugMfile(params,bulk_in,bulk_out,contrast)");*/
	
    mxDestroyArray(MYFUN);

	result = engGetVariable(ep,"total");
	s = (double *)mxGetData(result);
	memcpy( sum, s, mxGetNumberOfElements(result)*sizeof(double) );

	lays = engGetVariable(ep,"layers");
	layers = (double *)mxGetData(lays);
	memcpy(ans,layers,mxGetNumberOfElements(lays)*sizeof(double) );

	/*engClose(ep)*/;    

	mxDestroyArray(PARAMS);
	mxDestroyArray(BULKIN);
	mxDestroyArray(BULKOUT);
	mxDestroyArray(CONTRAST);

};