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

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 函数的额外参数来选择函数。伪代码如下：

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

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

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

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