传递函数句柄作为 mex for Matlab 的输入

Question

我最近在使用 MATLAB 研究有限元方法

我尝试在 MATLAB 中优化我的代码。

在搜索的过程中，我发现使用 mex 函数可以加速矩阵组装。

在将我的 "PoissonAssembler.m" 移植到 mex 函数时，我遇到了一些问题。

PoissonAssembler.m基本就是这种结构

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
        }
    }
}

在 matlab 中，我有

f = @(x,y) some_math_function;

由于其他数值模拟会更改此功能，

我想使用函数句柄作为 mex 文件的输入

我发现有一种方法可以通过使用 mexCallMatlab() 和 feval 来实现。

但是，我认为它会降低我的代码速度，因为调用 matlab 会产生开销。

有什么方法可以避免每次更改函数句柄时不编译mex文件吗？

更精确的代码是这样的

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（和 C++）支持将函数作为变量和参数的概念，就像 MATLAB 函数句柄一样。如果您可以将支持的函数限制为某个合理的集合，则可以使用传递给 MEX 函数的额外参数来 select 函数。伪代码如下：

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 {
      ...
    }
}

这样，您可以使用 select 边缘函数的字符串以及该函数需要的一些参数来调用您的函数。

如果您可以在您的情况下使用 C++，尤其是 C++11，那么使用 std::function 和 bind 以及 lambda 可以使这类事情变得容易得多。如果可能的话，值得弄清楚这一点。在这种情况下，您还可以为函数定义一个 std::map 字符串名称以便于查找。

最后，我要指出的是，这些实现中的任何一个都可以让您为您已经了解的 MATLAB 回调定义一个额外的函数。手写选择快，真正一般情况下慢。两全其美。

传递函数句柄作为 mex for Matlab 的输入

Passing function handle as input for mex for Matlab

matlab

mex

function-handle

finite-element-analysis