eigen（一）简介

最新推荐文章于 2022-10-14 16:25:22 发布

过得很好梵蒂冈

最新推荐文章于 2022-10-14 16:25:22 发布

阅读量1.2k

点赞数

本文链接：https://blog.csdn.net/qq_44944990/article/details/89429048

版权

一、概况

Eigen是有关线性代数（矩阵、向量等）的c++模板库。支持SSE2/3/4, ARM NEON (32-bit and 64-bit), PowerPC AltiVec/VSX (32-bit and 64-bit) instruction sets, S390x SIMD (ZVector)

（一）安装

所有的源码在头文件，无需编译

（二）优化

使用vectorization
x86使用-msse2编译选项
x86-64默认支持SSE2
32-bit ARM NEON使用-mfpu=neon -mfloat-abi=softfp|hard
64-bit ARM SIMD默认支持
小矩阵（维度为2-4）使用fixed-size
注意temporary objects的创建
矩阵表达式的写法

二、使用流程

1.包含头文件
头文件模块，简单的矩阵乘法使用#include

#include <iostream>#include <Eigen/Dense>using Eigen::MatrixXd;int main(){  MatrixXd m(2,2);  m(0,0) = 3;  m(1,0) = 2.5;  m(0,1) = -1;  m(1,1) = m(1,0) + m(0,1);  std::cout << m << std::endl;}
  
  1
2
3
4
5
6
7
8
9
10
11
12

2.编译指定eigen在本机的安装路径

g++ -I /path/to/eigen/ my_program.cpp -o my_program
  
  1

三、数据类型

数据类型包括两类Matrix和Array。

typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, Options> MyMatrixType;typedef Array<Scalar, RowsAtCompileTime, ColsAtCompileTime, Options> MyArrayType;
  
  1
2

存储顺序默认使用列存储，可以使用data()函数获得存储数据的首地址.

1. Matrix

Matrix

//使用caseMatrix<double, 6, Dynamic>                  // Dynamic number of columns (heap allocation)Matrix<double, Dynamic, 2>                  // Dynamic number of rows (heap allocation)Matrix<double, Dynamic, Dynamic, RowMajor>  // Fully dynamic, row major (heap allocation)Matrix<double, 13, 3>                       // Fully fixed (usually allocated on stack)//typedef简化Matrix<float,Dynamic,Dynamic>   <=>   MatrixXfMatrix<double,Dynamic,1>        <=>   VectorXdMatrix<int,1,Dynamic>           <=>   RowVectorXiMatrix<float,3,3>               <=>   Matrix3fMatrix<float,4,1>               <=>   Vector4f
  
  1
2
3
4
5
6
7
8
9
10
11

2.Array

定义类似于Matrix，支持coefficient-wise运算，比如每个元素加一个常数；两个矩阵对应的元素相乘。

//typdef 简化Array<float,Dynamic,Dynamic>    <=>   ArrayXXfArray<double,Dynamic,1>         <=>   ArrayXdArray<int,1,Dynamic>            <=>   RowArrayXiArray<float,3,3>                <=>   Array33fArray<float,4,1>                <=>   Array4f
  
  1
2
3
4
5
6

3.Matrix和Array之间相互转化

matrix=>array: matrix的.array()函数
array=>matrix: array的.matrix()函数
eigen不允许在表达式中混合使用matrix和array，但是允许使用=进行隐式转化，比如：

MatrixXf m(2,2);MatrixXf n(2,2);MatrixXf result(2,2);result = m.array() * n.array();
  
  1
2
3
4

matrix自身有成员函数.cwiseProduct()可以执行coefficient-wise product。

四、初始化

1.构造函数

Vector4d  v4;Vector2f  v1(x, y);Array3i   v2(x, y, z);Vector4d  v3(x, y, z, w);MatrixXf  m5; // empty objectMatrixXf  m6(nb_rows, nb_columns);
  
  1
2
3
4
5
6

2.comma

Vector3f  v1;     v1 << x, y, z;ArrayXf   v2(4);  v2 << 1, 2, 3, 4;Matrix3f  m1;   m1 << 1, 2, 3,                      4, 5, 6,                      7, 8, 9;
  
  1
2
3
4
5

3. map

对于原始的数据，可以使用Map复用原来的数据地址。

//pf为float*Map<MatrixXf> mf(pf,rows,columns);//pi为int*Map<const Vector4i> mi(pi);
  
  1
2
3
4

使用map

float data[] = {1,2,3,4};Map<Vector3f> v1(data);       // uses v1 as a Vector3f objectMap<ArrayXf>  v2(data,3);     // uses v2 as a ArrayXf objectMap<Array22f> m1(data);       // uses m1 as a Array22f objectMap<MatrixXf> m2(data,2,2);   // uses m2 as a MatrixXf objectfloat data[] = {1,2,3,4,5,6,7,8,9};Map<VectorXf,0,InnerStride<2> >  v1(data,3);                      // = [1,3,5]Map<VectorXf,0,InnerStride<> >   v2(data,3,InnerStride<>(3));     // = [1,4,7]Map<MatrixXf,0,OuterStride<3> >  m2(data,2,3);                    // both lines     |1,4,7|Map<MatrixXf,0,OuterStride<> >   m1(data,2,3,OuterStride<>(3));   // are equal to:  |2,5,8
  
  1
2
3
4
5
6
7
8
9
10
11

map和原始数组之间的相互转化

//array=>Matrixdouble *x;MatrixXd m = Map<MatrixXd> (x, rows, cols);//Matrix=>arraydouble *y;MatrixXd n;Map<MatrixXd>(y, n.rows(), n.cols())=n;
  
  1
2
3
4
5
6
7

五、运算

1.aliasing

aliasing：如果等号的两边出现同一个matrix或者array，就会引起结果异常。
1.对于标量乘法、加减法，没有影响

mat = 2 * mat;      //no aliasingarr = arr.square(); //no aliasing
  
  1
2

2.只有matrix乘法，eigen认定肯定会出现aliasing
默认引入临时matrix来存放中间结果

//1.aliasingmatA=matA*matA//eigen默认将上式分解成temp=matA*matAmatA=temp//2.no aliasing//使用noalias()来避免上面的分解matB=matA*matA //仍然会分解，虽然没有必要matB.noalias()=matA*matA //不会分解
  
  1
2
3
4
5
6
7
8
9
10

3.除了矩阵乘法以外的运算如果会出现aliasing问题，可以使用eval()函数或者xxxInPlace()函数避免

Matrix2i a; a << 1, 2, 3, 4;a = a.transpose(); //error: aliasinga = a.transpose().eval(); //oka.transposeInPlace(); //ok
  
  1
2
3
4

2.尽量使用复杂表达式，eigen可以更好地优化

六、多线程

1.支持多线程

eigen的部分算法（矩阵乘法、PartialPivLU）支持多线程，编译的时候增加-fopenmp，可以通过以下方式设置使用的线程数：

OMP_NUM_THREADS=n ./my_programomp_set_num_threads(n);Eigen::setNbThreads(n);
  
  1
2
3

也可以通过EIGEN_DONT_PARALLELIZE宏定义关闭eigen的多线程。
如果程序本身使用了多线程，需要使用以下方式初始化eigen

#include <Eigen/Core>int main(int argc, char** argv){  Eigen::initParallel();  ...}
  
  1
2
3
4
5
6
7

2.使用mkl库

Intel mkl(Math Kernal Library)可以支持更高程度优化的多线程数学计算。eigen使用MKL流程：
- define the EIGEN_USE_MKL_ALL macro before including any Eigen’s header
- link your program to MKL libraries (see the MKL linking advisor)
- on a 64bits system, you must use the LP64 interface (not the ILP64 one)
MKL的优化适用于Dynamic的matrix，数据类型仅限于float/double/complex/complex，其他类型的数据或者real和complex混用的数据不会使用mkl优化，而是使用内置的算法优化。