开发笔记：Eigen开源库的应用

发表于 更新于 分类于 阅读数： 评论数： 本文字数： 2k 阅读时长 ≈ 3 分钟
Eigen开源库的入门

Eigen与C之间数据转换

to Eigen

1
2
3
4
float* raw_data = malloc(...);
Map<MatrixXd> M(raw_data, rows, cols);
// use M as a MatrixXd
M = M.inverse();

from Eigen

1
2
3
4
MatrixXd M;
float* raw_data = M.data();
int stride = M.outerStride();
raw_data[i+j*stride]

一些预备知识

template programming

4 levels 并行

  • cluster of PCs --MPI
  • multi/many-cores -- OpenMP
  • SIMD -- intrinsics for vector instructions
  • pipelining -- needs non dependent instructions

Peak Performance

Example： Intel Core2 Quad CPU Q9400 @ 2.66GHz (x86_64)

* pipelining → 1 mul + 1 add / cycle (ideal case)
* SSE → x 4 single precision ops at once
* frequency → x 2.66G
* peak performance: 21,790 Mflops (for 1 core)

这就是我们优化的目标

Fused operation: Expression Template

Expression Template是个好东西。通过编译融合嵌入的方式，减少了大量的读写和计算。

Curiously Recurring Template Pattern

Eigen常见的坑

编译的时候最好-DEIGEN_MPL2_ONLY(详见: Eigen)

这是因为Eigen虽然大部分是MPL2 licensed的，但是还有少部分代码是LGPL licensed的，如果修改了其代码，则必须开源。 这可能产生法律风险，而遭到法务部门的Complain

要注意alignment的问题

1
2
3
4
5
6
7
my_program: path/to/eigen/Eigen/src/Core/DenseStorage.h:44:
Eigen::internal::matrix_array<T, Size, MatrixOptions, Align>::internal::matrix_array()
[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
Assertion `(reinterpret_cast<size_t>(array) & (sizemask)) == 0 && "this assertion
is explained here: http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html
     READ THIS WEB PAGE !!! ****"' failed.
There are 4 known causes for this issue. Please read on to understand them and learn how to fix them.

例如下面的代码都是有问题的，可能导致整个程序Crash。

  • 结构体含有Eigen类型的成员变量
1
2
3
4
5
6
7
class Foo {
  //...
  Eigen::Vector2d v; // 这个类型只是一个例子，很多类型都有问题
  //...
};
//...
Foo *foo = new Foo;
  • Eigen类型的变量被放到STL容器中
1
2
3
4
// Eigen::Matrix2f这个类型只是一个例子，很多类型都有问题
std::vector<Eigen::Matrix2f> my_vector;
struct my_class { ... Eigen::Matrix2f m; ... }; 
std::map<int, my_class> my_map;
  • Eigen类型的变量被按值传入一个函数
1
2
// Eigen::Vector4d只是一个例子，很多类型都有这个问题
void func(Eigen::Vector4d v);
  • 在栈上定义Eigen类型的变量(只有GCC4.4及以下版本 on Windows被发现有这个问题，例如MinGW or TDM-GCC)
1
2
3
void foo() {
  Eigen::Quaternionf q; // Eigen::Quaternionf只是一个例子，很多类型都有这个问题
}