2016-07-08 6 views
1

Я недавно работал на методе конечных элементов с MATLABПередача функции ручки в качестве входных данных для Mex для 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 файла

I обнаружили, что есть способ сделать это, используя 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 
     } 
    } 
} 
+0

Сколько возможных функций вам нужно передать? Один из вариантов - векторизация вашей функции, чтобы вы могли передавать все грани/элементы одновременно, а не вызывать ее в цикле. – Suever

+1

Лучший сценарий здесь, вероятно, создает вашу mex-функцию с корпусом коммутатора для возможных функций .... Преимущество MATLAB в том, что оно очень гибкое, преимущество C в том, что он очень быстрый. Вы не можете иметь оба, Im fear –

ответ

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 этот вид вещей стал намного проще и лямбдами. Стоит подумать об этом, если это возможно. Вы также можете в этом случае определить std::map имен строк для функций для удобного поиска.

Наконец, я укажу, что любая из этих реализаций позволит вам определить одну дополнительную функцию для обратного вызова MATLAB, который вы уже выяснили. Быстро для рукописных вариантов, медленных для действительно общего случая. Лучшее обоих миров.

+0

Спасибо за ваш комментарий! – Dohyun