使C++特征LU更快(我的测试显示比GSL慢2倍）

我将特征的LU分解/求解与GSL进行比较，发现特征在g++/OSX上与-O3优化相比慢了2倍。我隔离了分解和求解的时间，但发现两者都比它们的GSL同行慢得多。我是在做一些愚蠢的事情，还是Eigen在这个用例中表现不佳（例如，非常小的系统？）使用Eigen 3.2.8和旧的GSL 1.15构建。测试用例非常做作，但反映了我正在编写的一些非线性拟合软件中的结果——对于总的LU/求解操作，Eigen的速度在1.5x-2x+之间。

#define NDEBUG
#include "sys/time.h"
#include "gsl/gsl_linalg.h"
#include <Eigen/LU>
// Ax=b is a 3x3 system for which soln is x=[8,2,3]
//
double avals_col[9] = { 4, 2, 2, 7, 5, 5, 7, 5, 9 };
    // col major
double avals_row[9] = { 4, 7, 7, 2, 5, 5, 2, 5, 9 };
    // row major
double bvals[9] = { 67, 41, 53 };
//----------- helpers
void print_solution( double *x, int dim, char *which ) {
    printf( "%s solve for x:n", which );
    for( int j=0; j<3; j++ ) {
        printf( "%g ", x[j] );
    }
    printf( "n" );
}
struct timeval tv;
struct timezone tz;
double timeNow() {
    gettimeofday( &tv, &tz );
    int _mils = tv.tv_usec/1000;
    int _secs = tv.tv_sec;
    return (double)_secs + ((double)_mils/1000.0);
}
//-----------
void run_gsl( double *A, double *b, double *x, int dim, int reps ) {
    gsl_matrix_view gslA;
    gsl_vector_view gslB;
    gsl_vector_view gslX;
    gsl_permutation *gslP;
    int sign;
    gslA = gsl_matrix_view_array( A, dim, dim );
    gslP = gsl_permutation_alloc( dim );
    gslB = gsl_vector_view_array( b, dim );
    gslX = gsl_vector_view_array( x, dim );
    int err;
    double t, elapsed;
    t = timeNow();
    for( int i=0; i<reps; i++ ) {
        // gsl overwrites A during decompose, so we must copy the initial A each time.
        memcpy( A, avals_row, sizeof(avals_row) );
        err = gsl_linalg_LU_decomp( &gslA.matrix, gslP, &sign );
    }
    elapsed = timeNow() - t;
    printf( "GSL decompose (%dx) time = %gn", reps, elapsed );
    t = timeNow();
    for( int i=0; i<reps; i++ ) {
        err = gsl_linalg_LU_solve( &gslA.matrix, gslP, &gslB.vector, &gslX.vector );
    }
    elapsed = timeNow() - t;
    printf( "GSL solve (%dx) time = %gn", reps, elapsed );
    gsl_permutation_free( gslP );
}
void run_eigen( double *A, double *b, double *x, int dim, int reps ) {
    Eigen::PartialPivLU<Eigen::MatrixXd> eigenA_lu;
    Eigen::Map< Eigen::Matrix < double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor > > ma( A, dim, dim );
    Eigen::Map<Eigen::MatrixXd> mb( b, dim, 1 );
    int err;
    double t, elapsed;
    t = timeNow();
    for( int i=0; i<reps; i++ ) {
        // This memcpy is not necessary for Eigen, which does not overwrite A in the
        // decompose, but do it so that the time is directly comparable to GSL. 
        memcpy( A, avals_col, sizeof(avals_col) );
        eigenA_lu.compute( ma );
    }
    elapsed = timeNow() - t;
    printf( "Eigen decompose (%dx) time = %gn", reps, elapsed );
    t = timeNow();
    Eigen::VectorXd _x;
    for( int i=0; i<reps; i++ ) {
         _x = eigenA_lu.solve( mb );
    }
    elapsed = timeNow() - t;
    printf( "Eigen solve (%dx) time = %gn", reps, elapsed );
    // copy soln to passed x
    for( int i=0; i<dim; i++ ) {
        x[i] = _x(i);
    }
}
int main() {
    // solve a 3x3 system many times
    double A[9], b[3], x[3];
    int dim = 3;
    int reps = 1000000;
    memcpy( b, bvals, sizeof(bvals) );
        // init b vector, A is copied multiple times in run_gsl/run_eigen
    run_eigen( A, b, x, dim, reps );
    print_solution( x, dim, "Eigen" );  
    run_gsl( A, b, x, dim, reps );
    print_solution( x, dim, "GSL" );
    return 0;
}

这会产生，例如：

Eigen decompose (1000000x) time = 0.242
Eigen solve (1000000x) time = 0.108
Eigen solve for x:
8 2 3 
GSL decompose (1000000x) time = 0.049
GSL solve (1000000x) time = 0.075
GSL solve for x:
8 2 3