2016-07-08 2 views
1
MEX

Je travaille récemment sur la méthode des éléments finis avec MatlabPassing fonction poignée comme entrée pour Matlab pour

j'ai essayé d'optimiser mon code en MATLAB.

Lors de la recherche, j'ai trouvé que l'assemblage de la matrice pourrait être accéléré en utilisant la fonction mex. Pendant le portage de mon "PoissonAssembler.m" à la fonction mex, j'ai rencontré quelques problèmes. PoissonAssembler.m a fondamentalement ce genre de 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 
     } 
    } 
} 

Dans Matlab, je

f = @(x,y) some_math_function; 

Étant donné que cette fonction sera modifiée pour d'autres simulation numérique,

Je voulais utiliser la poignée de fonction comme entrée pour le fichier MEX

I trouvé qu'il y a un moyen de le faire en utilisant mexCallMatlab() et feval. Cependant, je pense que cela va ralentir mon code à cause de la surcharge causée par l'appel de matlab.

Y at-il un moyen de l'éviter sans compiler le fichier mex chaque fois que je change la poignée de fonction?

code plus précise est comme ce

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

De combien de fonctions avez-vous vraiment besoin pour le transmettre? Une option consiste à vectoriser votre fonction afin que vous puissiez passer tous les visages/éléments à la fois plutôt que de l'appeler dans une boucle. – Suever

+1

Le meilleur scénario ici est probablement de créer votre fonction mex avec un switch pour les fonctions possibles .... L'avantage de MATLAB est qu'il est très flexible, l'avantage de C est qu'il est très rapide. Vous ne pouvez pas avoir les deux, j'ai peur –

Répondre

1

C (et C++) soutiennent le concept de fonctions en tant que variables et paramètres beaucoup comme les poignées fonction Matlab. Si vous pouvez limiter les fonctions prises en charge à un ensemble raisonnable, vous pouvez utiliser des arguments supplémentaires passés à la fonction MEX pour sélectionner la fonction. Pseudocode suit:

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

De cette façon, vous appelez votre fonction avec une chaîne qui sélectionne votre fonction de bord, ainsi que certains paramètres que le besoin de la fonction.

Si vous pouvez utiliser C++ dans votre situation, et en particulier C++ 11, ce genre de chose est devenu beaucoup plus facile avec std::function et bind et lambdas. Ça vaut le coup de le faire si possible. Vous pouvez également dans ce cas définir un std::map de noms de chaînes à des fonctions pour faciliter les recherches.

Enfin, je ferai remarquer que l'une ou l'autre de ces implémentations vous permettrait de définir une fonction supplémentaire pour le rappel MATLAB que vous avez déjà compris. Rapide pour les choix manuscrits, lent pour le cas très général. Le meilleur des deux mondes.

+0

Merci pour votre commentaire! – Dohyun