Following the suggestions from this post Identifying temporary object creation in Eigen, I made some tests to understand when temporary Eigen objects are created. This is the code used to identify temporaries:
static long int nb_temporaries;
inline void on_temporary_creation( long int size ) {
if ( size != 0 ) nb_temporaries++;
}
#define EIGEN_DENSE_STORAGE_CTOR_PLUGIN { on_temporary_creation(size); }
All tests has been perfomed using VisualStudio 2017 (v15.9.20) with Eigen 3.3.7 and the following objects:
const int n = 100;
MatrixXd A( n, n );
VectorXd x( n );
VectorXd y( n );
VectorXd a( n );
VectorXd b( n );
Map<MatrixXd> A_map( A.data(), n, n );
Map<VectorXd> x_map( x.data(), n );
Map<VectorXd> y_map( y.data(), n );
Map<VectorXd> a_map( a.data(), n );
Map<VectorXd> b_map( b.data(), n );
SparseMatrix<double> As( n, n );
As.reserve( nnz );
As.makeCompressed();
Map<SparseMatrix<double>> As_map( n, n, nnz,
As.outerIndexPtr(), As.innerIndexPtr(), As.valuePtr(), nullptr );
Matrix2d A2;
Vector2d x2;
Vector2d b2;
Map<Matrix2d> A2_map( A2.data() );
Map<Vector2d> x2_map( x2.data() );
Map<Vector2d> b2_map( b2.data() );
Test A: coefficient wise operations
None of the following create temporaries:
x = a + b;
x.array() = a.array() / b.array();
x_map.array() = a_map.array() / b_map.array();
Test B: Dense Matrix-Vector product
Both the following create 1 temporary with size n:
y.noalias() = A * x + b;
y_map.noalias() = A_map * x_map + b_map;
To avoid it you should do:
y.noalias() = A.lazyProduct( x ) + b;
y_map.noalias() = A_map.lazyProduct( x_map ) + b_map;
Both the following create 1 temporary with size n:
y.noalias() = A.selfadjointView<Eigen::Lower>() * x + b;
y_map.noalias() = A_map.selfadjointView<Eigen::Lower>() * x_map + b_map;
To avoid it you should do:
y.noalias() = A.selfadjointView<Eigen::Lower>() * x;
y += b;
and
y_map.noalias() = A_map.selfadjointView<Eigen::Lower>() * x_map;
y_map += b_map;
Test C: Sparse Matrix-Vector product
The following creates 1 temporary of size n:
y.noalias() = As * x + b;
To avoid it you should do:
y.noalias() = As * x;
y += b;
The following creates 4 temporary of size (2, 2, n, 2):
y_map.noalias() = As_map * x_map + b_map;
To avoid dynamic allocated temporaries you should do (there are still 2 stack allocated temporaries with size 2):
y_map.noalias() = As_map * x_map;
y_map += b_map;
The following creates 1 temporary of size n:
y.noalias() = As.selfadjointView<Eigen::Lower>() * x + b;
To avoid it you should do:
y.noalias() = As.selfadjointView<Eigen::Lower>() * x;
y += b;
The following creates 6 temporaries with size (2, 2, 2, n, 2, 2):
y_map.noalias() = As_map.selfadjointView<Eigen::Lower>() * x_map + b_map;
To avoid dynamic allocated temporaries you should do (there are still 4 stack allocated temporaries with size 2):
y_map.noalias() = As_map.selfadjointView<Eigen::Lower>() * x_map;
y_map += b_map;
Test D: QR solve with 2x2 fixed matrix
The following creates 5 temporaries with size (4, 2, 2, 2, 1):
x2.noalias() = A2.householderQr().solve( b2 );
The following creates 6 temporaries with size (4, 4, 2, 2, 2, 1):
x2_map.noalias() = A2_map.householderQr().solve( b2_map );
Test E: Sparse Matrix LDLT solve
The following creates 8 temporaries of size (n, n, n+1, n, n, n, n, n) :
SimplicialLDLT<SparseMatrix<double>> solver;
solver.analyzePattern( As );
The following creates 2 temporaries of size (n, n):
solver.factorize( As );
If inplace factorization is used (see Inplace matrix decompositions) then factorize method creates 1 temporary of size n. In this case explicit call of factorize method is useful only if some elements of A have changed but the sparsity pattern is still the same.
The following does not create temporaries:
x = solver.solve( b );
Summing up:
- Coefficients wise operations don't create temporaries
- Dense matrix-vector product doesn't create temporaries if lazyProduct and noalias are used
- Sparse matrix does not provide lazyProduct
- Using Map more temporaries are created
Regarding the above cases are there any tricks to better perform and avoid temporaries?
