Would it be possible to implement a class that receives a C-style pointer as a template argument and somehow resolves into a static Eigen matrix but using the memory provided? Say a declaration would look something like:
EIGEN_ALIGN16 double array[9];
CMatrix<double,3,3,array> :: m;
I do know about maps, but the example code I provide below has proven them to be slower by 20% when compared to static Eigen matrices.
These would be the premises:
- I need to provide my own C pointer. This way I can efficiently reuse C code without incurring copies.
- The resulting matrix should look static to Eigen so that Eigen can optimize as it would with a static array at compile time. Look at my example above where at compile time I would provide both matrix size (static) and the C pointer.
- CMatrix should fall back to Eigen::Matrix. When the additional template parameter for the C array is not provided I would get the normal Eigen matrix.
- I do not intend to make a full Eigen extension. With that I mean is I do not care about all kind of checks to provide a neat extension for other users. I just want the most efficient solution for this problem
Would it be possible to implement a solution by adding a new constructor? Say something like:
EIGEN_ALIGN16 double data[9];
Eigen::Matrix<double,3,3> m(data); //where data is NOT copied but used to replace the static allocation used by default.
Find below my code for benchmarking Map vs. Matrix efficiency. It is self contained and you can compile with:
g++ -Ofast -DNDEBUG -DEIGEN_NO_MALLOC -I/path_to_my_Eigen benchmark1.cpp -o benchmark1 -lrt
Here is the code:
#include <Eigen/Eigen>
#include <bench/BenchTimer.h>
#include <iostream>
using namespace Eigen;
using namespace std;
//#define CLASSIC_METHOD
#define USE_MAPS
EIGEN_DONT_INLINE void classic(double VO[4], double AT[4][4], double VI[4])
{
for (int ii=0; ii<4; ii++)
{
VO[ii] = AT[ii][0] * VI[0] +
AT[ii][1] * VI[1] +
AT[ii][2] * VI[2] +
AT[ii][3] * VI[3];
}
};
template <typename OutputType, typename MatrixType, typename VectorType>
EIGEN_DONT_INLINE void modern(MatrixBase<OutputType>& VOE, const MatrixBase<MatrixType>& A44, const MatrixBase<VectorType>& VIE)
{
VOE.noalias() = A44.transpose()*VIE;
};
int main()
{
EIGEN_ALIGN16 double AT[4][4] = {0.1, 0.2, 0.3, 2.0, 0.2, 0.3, 0.4, 3.0, 0.3, 0.4, 0.5, 4.0, 0.0, 0.0, 0.0, 1.0};
EIGEN_ALIGN16 double VI[4] = {1, 2, 3, 4};
EIGEN_ALIGN16 double VO[4];
//Eigen matrices
#ifndef USE_MAPS
Matrix4d A44 = Matrix4d::MapAligned(AT[0]);
Vector4d VIE = Vector4d::MapAligned(VI);
Vector4d VOE(0,0,0,0);
#else
Map<Matrix4d,Aligned> A44(AT[0]);
Map<Vector4d,Aligned> VIE(VI);
Map<Vector4d,Aligned> VOE(VO);
// Map<Matrix4d> A44(AT[0]);
// Map<Vector4d> VIE(VI);
// Map<Vector4d> VOE(VO);
#endif
#ifdef EIGEN_VECTORIZE
cout << "EIGEN_VECTORIZE defined" << endl;
#else
cout << "EIGEN_VECTORIZE NOT defined" << endl;
#endif
cout << "VIE:" << endl;
cout << VIE << endl;
VI[0] = 3.14;
cout << "VIE:" << endl;
cout << VIE << endl;
BenchTimer timer;
const int num_tries = 5;
const int num_repetitions = 200000000;
#ifdef CLASSIC_METHOD
BENCH(timer, num_tries, num_repetitions, classic(VO, AT, VI));
std::cout << Vector4d::MapAligned(VO) << std::endl;
#else
BENCH(timer, num_tries, num_repetitions, modern(VOE, A44, VIE));
std::cout << VOE << std::endl;
#endif
double elapsed = timer.best();
std::cout << "elapsed time: " << elapsed*1000.0 << " ms" << std::endl;
return 0;
}
