Passing function handle as input for mex for Matlab

Question

I'm recently working on finite element method with MATLAB

I tried to optimize my code in MATLAB.

While searching, I found that matrix assembling could be accelerated by using mex function.

While porting my "PoissonAssembler.m" to mex function, I stuck into some problems.

PoissonAssembler.m has basically this kind of structure.

for every element{

    loc_stiff_assmbl();

    loc_rhs_assmbl(); // this one involves with function evaluation @(x,y)

    for every edges{
        loc_edge_assmbl();

        if boundary{
            loc_rhs_edge_assmbl(); // this one also have function eval
        }
    }
}

In matlab, I have

f = @(x,y) some_math_function;

Since this function will be changed for other numerical simulation,

I wanted to use function handle as input for mex file

I found that there is a way to do this by using mexCallMatlab() and feval.

However, I think it will slow down my code because of overhead caused by calling matlab.

Is there any way to avoid it without compiling mex file every time I change the function handle?

More precise code is like this

for (k=0;k<nE;k++){ // nE : number of elements (about 10^5 to 10^6)
    // other assembling procedure
    blabla

    // function evaluation for rhs assemble
    for (i=0; i<nP; i++){ // nP : number of points in element
        fxy[i] = f(x[i],y[i])};
    }

    // rhs assemble
    for (j=0; j<nP; j++){
        for (i=0; i<nP; i++){ // matrix vector multiplication
            rhs[k*nE + j] += Mass[ i*nP + j]*fxy[i];
        }
    }

    // face loop
    for (f1 = 0; f1 < nF4E; f1++){ // number of face for each elements
        switch(BCType[k1*nF4E + f1]){ // if boundary

            case 1 : // Dirichlet
                // inner product with evaluation of gD = @(x,y) function
                CODE OMITTED
            case 2 : // Neumann
                // inner product with evaluation of gN = @(x,y) function
                CODE OMITTED
        }
    }
}

Answer 1

C (and C++) support the concept of functions as variables and parameters much like MATLAB function handles. If you can limit the functions supported to some reasonable set, you can use extra arguments passed to the MEX function to select the function. Pseudocode follows:

double fn_linear(double x, double y, double *args)
{
    return arg[0] + x*arg[1] + y*arg[2];
}

double fn_quadratic(double x, double y, double *args)
{
    return arg[0] + x*arg[1] + x*x*arg[2] + /// etc
}

typedef double (*EdgeFunction)(double x, double y, double *args);


void computation_function(other parms, EdgeFunction edge_fcn, double *fcn_parms)
{
    for (k=0;k<nE;k++){ // nE : number of elements (about 10^5 to 10^6)
        // other assembling procedure
        blabla

        // function evaluation for rhs assemble
        for (i=0; i<nP; i++){ // nP : number of points in element
            fxy[i] = (*edge_fcn)(x[i], y[i], fcn_parms);
        }
}

void mexFunction()
{
    EdgeFunction edge_fcn;

    if(strcmp(arg3, "linear")==0) {
       edge_fcn = &fn_linear;
       extra_parms = arg4;
    } else {
      ...
    }
}

This way, you call your function with a string that selects your edge function, plus some parameters that the function need.

If you can use C++ in your situation, and especially C++11, this sort of thing has gotten much easier with std::function and bind and lambdas. It's worth figuring that out if possible. You could also in this case define a std::map of string names to functions for easy lookups.

Finally, I'll point out that either of these implementations would let you define one extra function for the MATLAB callback you've already figured out. Fast for the hand-written choices, slow for the really general case. Best of both worlds.