Matrix mal Vektor Betrieb mit rohen Puffer in Eigen

Question

Ich möchte folgende generische Matrix * Vektorfunktion

void MV(double *M, double *x, double *y, int n, int m) 
{ 
    for (int i = 0; i<n; i++) { 
     for (int j = 0; j<m; j++) { 
      y[i] += M[i*m + j] * x[j]; 
     } 
    } 
} 

mit einer ersetzen, die Vektoroperationen Eigens Matrix * verwendet. Das habe ich getan.

void MV2(double *M, double *x, double *y, int n, int m) 
{ 
    typedef Eigen::Matrix<double, -1, -1, Eigen::RowMajor> Rm; 
    Eigen::Map<const Rm> Me(M, n, m); 
    Eigen::Map<const Rm> xe(x, m, 1); 
    Eigen::Map<Rm> ye(y, n, 1); 
    ye = Me*xe; 
} 

Mit einem einfachen "Hallo Welt" Test (mit GCC 5.3, Ubuntu 16.04) das funktioniert gut

int main(void) 
{ 
    double M[12]; 
    double x[4]; 
    double y[3]; 
    for(int i=0; i<12; i++) M[i] = i; 
    for(int i=0; i<4; i++) x[i] = i; 
    int n = 3, m = 4; 

    MV(M,x,y,3,4); 
    printf("%f %f %f\n", y[0], y[1], y[2]); 
    MV2(M,x,y,3,4); 
    printf("%f %f %f\n", y[0], y[1], y[2]); 
} 

Allerdings habe ich jetzt diese Funktionen in meiner Hauptanwendung bin mit Visual Studio 2013 erstellt Ich laufe über viele kleine Matrizen (zB 60x60) verschiedener Größen aus einer spärlichen Matrix aus vielen kleinen dichten Matrizen. Die Funktion MV liefert das korrekte Ergebnis, die Funktion MV2 jedoch nicht. Es ist auch langsamer. Was mache ich falsch?

Answer 1

Sie haben die += in der MV2-Version vergessen. Dann Leistung in Bezug auf, besser lassen Eigen weiß, wann Sie Vektoren bei der Kompilierung haben:

void MV2(const double *M, const double *x, double *y, int n, int m) 
{ 
    typedef Eigen::Matrix<double, -1, -1, Eigen::RowMajor> Rm; 
    Eigen::VectorXd::Map(y,n).noalias() += Rm::Map(M,n,m) 
             * Eigen::VectorXd::Map(x,m); 
} 

Die noalias erlaubt eine vorübergehende zu vermeiden, indem man Eigen weiß, dass das Ergebnis in y geschrieben werden kann, ohne Probleme Aliasing.

EDIT

Wenn M symmetrisch ist, dann können Sie es zu Eigen, wie folgt:

VectorXd::Map(y,n).noalias() += Rm::Map(M,n,m).selfadjointView<Lower>() 
           * VectorXd::Map(x,m);