[Eigen中文文档] 深入了解 Eigen - Eigen内部发生了什么(三)

简介: 本页的目标是了解 Eigen 如何编译。

文档总目录

英文原文(What happens inside Eigen, on a simple example)

至此,表达式 v + w 已完成求值,接下来,在编译该行代码的过程中,输入运算符 =

u = v + w;

这里调用了运算符 =,向量 uVectorXf类的对象,即 Matrixsrc/Core/Matrix.h 中,在 Matrix 类的定义中,我们看到:

template<typename OtherDerived>
inline Matrix& operator=(const MatrixBase<OtherDerived>& other)
{
   
    eigen_assert(m_storage.data()!=0 && "you cannot use operator= with a non initialized matrix (instead use set()");
    return Base::operator=(other.derived());
}

这里,BaseMatrixBase<Matrix>typedef 类型。所以,所谓的 Base::operator= 就是 MatrixBase<Matrix>::operator= 。它的原型在 src/Core/MatrixBase.h 中:

template<typename OtherDerived>
Derived& operator=(const MatrixBase<OtherDerived>& other);

这里,DerivedVectorXf(因为u是一个VectorXf),而OtherDerivedCwiseBinaryOp。具体来说,如前面所述,OtherDerived是:

CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf>

所以被调用的operator=的完整原型是:

VectorXf& MatrixBase<VectorXf>::operator=(const MatrixBase<CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf> > & other);

operator=字面意思是将两个 VectorXf的总和复制到另一个 VectorXf中。这个operator=的实现在文件 src/Core/Assign.h 中:

template<typename Derived>
template<typename OtherDerived>
inline Derived& MatrixBase<Derived>
  ::operator=(const MatrixBase<OtherDerived>& other)
{
   
    return internal::assign_selector<Derived,OtherDerived>::run(derived(), other.derived());
}

接下来,我们了解一下 internal::assign_selector,这是它的声明(所有内容仍然在同一个文件 src/Core/Assign.h 中):

template<typename Derived, typename OtherDerived,
         bool EvalBeforeAssigning = int(OtherDerived::Flags) & EvalBeforeAssigningBit,
         bool NeedToTranspose = Derived::IsVectorAtCompileTime
            && OtherDerived::IsVectorAtCompileTime
            && int(Derived::RowsAtCompileTime) == int(OtherDerived::ColsAtCompileTime)
            && int(Derived::ColsAtCompileTime) == int(OtherDerived::RowsAtCompileTime)
            && int(Derived::SizeAtCompileTime) != 1>
struct internal::assign_selector;

internal::assign_selector有4个模板参数,但是后面的两个由前面的两个自动确定。

EvalBeforeAssigning用于强制执行EvalBeforeAssigningBit。某些表达式具有此标志,使它们在将其分配给另一个表达式之前自动计算为临时值。这是Product表达式的情况,为了避免在执行m = m * m;时出现奇怪的混叠现象。CwiseBinaryOp表达式在这里没有EvalBeforeAssigningBit:从一开始就说过不想在这里引入一个临时值。在 src/Core/CwiseBinaryOp.h 中,internal::traits<CwiseBinaryOp> 中的Flags不包括EvalBeforeAssigningBit。然后,通过EIGEN_GENERIC_PUBLIC_INTERFACE宏,CwiseBinaryOpFlags成员从internal::traits导入,在这里模板参数EvalBeforeAssigning的值为false

NeedToTranspose是为了用户想要将行向量复制到列向量的情况而设定的特殊例外。我们允许这种情况违反一般规则,即在分配中需要维度匹配。在这里左侧和右侧都是列向量,即ColsAtCompileTime等于1。因此,NeedToTranspose也是false

因此,在部分特化中:

internal::assign_selector<Derived, OtherDerived, false, false>

它的定义如下:

template<typename Derived, typename OtherDerived>
struct internal::assign_selector<Derived,OtherDerived,false,false> 
{
   
    static Derived& run(Derived& dst, const OtherDerived& other) 
    {
    return dst.lazyAssign(other.derived()); }
};

接下来,我们了解一下 lazyAssign 是如何工作的:

template<typename Derived>
template<typename OtherDerived>
inline Derived& MatrixBase<Derived>
  ::lazyAssign(const MatrixBase<OtherDerived>& other)
{
   
    EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived,OtherDerived)
    eigen_assert(rows() == other.rows() && cols() == other.cols());
    internal::assign_impl<Derived, OtherDerived>::run(derived(),other.derived());
    return derived();
}

其中,

internal::assign_impl<Derived, OtherDerived>::run(derived(),other.derived());

internal::assign_impl内部是怎么实现的?他的定义如下:

template<typename Derived1, typename Derived2,
     int Vectorization = internal::assign_traits<Derived1, Derived2>::Vectorization,
     int Unrolling = internal::assign_traits<Derived1, Derived2>::Unrolling>
struct internal::assign_impl;

同样,internal::assign_selector接受4个模板参数,但是后面的2个参数由前面的2个参数自动确定。

这两个参数VectorizationUnrolling由一个辅助类internal::assign_traits确定。它的作用是确定要使用的向量化策略(即Vectorization)和展开策略(即Unrolling)。

这里不会深入介绍如何选择这些策略(这在同一文件顶部的internal::assign_traits实现中)。只介绍这里Vectorization的值为LinearVectorization,而Unrolling的值为NoUnrolling(由于向量具有动态大小,因此没有办法在编译时展开循环,这一点是显而易见的)。

所以 internal::assign_impl 的部分特化是:

internal::assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>

他的定义如下:

template<typename Derived1, typename Derived2>
struct internal::assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
{
   
  static void run(Derived1 &dst, const Derived2 &src)
  {
   
    const int size = dst.size();
    const int packetSize = internal::packet_traits<typename Derived1::Scalar>::size;
    const int alignedStart = internal::assign_traits<Derived1,Derived2>::DstIsAligned ? 0
                           : internal::first_aligned(&dst.coeffRef(0), size);
    const int alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize;

    for(int index = 0; index < alignedStart; index++)
      dst.copyCoeff(index, src);

    for(int index = alignedStart; index < alignedEnd; index += packetSize)
    {
   
      dst.template copyPacket<Derived2, Aligned, internal::assign_traits<Derived1,Derived2>::SrcAlignment>(index, src);
    }

    for(int index = alignedEnd; index < size; index++)
      dst.copyCoeff(index, src);
  }
};

线性向量化意味着左侧和右侧表达式可以被线性访问,也就是你可以通过一个整数索引引用它们的系数,而不是必须通过两个整数行、列引用它们的系数。

正如我们在开始时所说,向量化是使用4个浮点数块工作的。在这里,PacketSize 为 4。

我们需要解决两个潜在的问题:

  • 首先,如果数据包是128位对齐的,向量化效果会更好。这对于写访问尤其重要。因此,在写入dst的系数时,我们希望通过4个数据包分组这些系数,以便每个数据包都是128位对齐的。一般来说,这需要跳过dst的前几个系数。这就是alignedStart的作用。然后,我们逐个复制这些前几个系数,而不是通过数据包。在我们的情况下,dst表达式是VectorXf,记住我们在向量的构造中分配了对齐的数组。Eigen的DstIsAligned已经处理了这一切,无需进行任何运行时检查,因此alignedStart为零。

  • 其次,要复制的系数数量通常不是packetSize的整数倍。在这里,要复制的系数有50个,packetSize4。因此,我们将不得不逐个复制最后2个系数,而不是通过数据包。在这里,alignedEnd48

接下来是实际的循环。

首先是向量化的部分:前50个系数中的前48个将通过4个数据包逐个复制:

for(int index = alignedStart; index < alignedEnd; index += packetSize)
{
   
    dst.template copyPacket<Derived2, Aligned, internal::assign_traits<Derived1,Derived2>::SrcAlignment>(index, src);
}

什么是 copyPacket? 它在 src/Core/Coeffs.h 中定义,如下:

template<typename Derived>
template<typename OtherDerived, int StoreMode, int LoadMode>
inline void MatrixBase<Derived>::copyPacket(int index, const MatrixBase<OtherDerived>& other)
{
   
    eigen_internal_assert(index >= 0 && index < size());
    derived().template writePacket<StoreMode>(index, other.derived().template packet<LoadMode>(index));
}

writePacket()packet() 是什么。

首先,这里的writePacket()是左侧VectorXf的一个方法。他的定义在 src/Core/Matrix.h 中:

template<int StoreMode>
inline void writePacket(int index, const PacketScalar& x)
{
   
    internal::pstoret<Scalar, PacketScalar, StoreMode>(m_storage.data() + index, x);
}

在这里,StoreModeAligned,表示我们正在进行128位对齐的写访问,PacketScalar是一个表示 4个浮点数的SSE数据包 的类型,internal::pstoret是一个函数,将这样的数据包写入内存。它们的定义是特定于架构的,我们可以在src/Core/arch/SSE/PacketMath.h中找到它们:

src/Core/arch/SSE/PacketMath.h中决定PacketScalar类型的代码行(通过Matrix.h中的typedef)是:

template<> struct internal::packet_traits<float>  
{
    typedef __m128  type; enum {
   size=4}; };

在这里,__m128是一个特定于SSE的类型。请注意,这里的enum size是用来定义上面的packetSize的。

下面是internal::pstoret的实现:

template<> inline void internal::pstore(float*  to, const __m128&  from) 
{
    _mm_store_ps(to, from); }

在这里,__mm_store_ps是一个特定于SSE的内部函数,表示单个SSE指令。internal::pstoreinternal::pstoret之间的区别在于internal::pstoret是一个调度程序,处理对齐和未对齐的情况,你可以在src/Core/GenericPacketMath.h中找到其定义:

template<typename Scalar, typename Packet, int LoadMode>
inline void internal::pstoret(Scalar* to, const Packet& from)
{
    if(LoadMode == Aligned)
        internal::pstore(to, from);
    else
        internal::pstoreu(to, from);
}

这解释了writePacket()的工作原理。现在让我们看看packet()调用:

derived().template writePacket<StoreMode>(index, other.derived().template packet<LoadMode>(index));

在这里,other是总和表达式v + w.derived()只是从MatrixBase转换为其子类,这里是CwiseBinaryOp。接下来查看 src/Core/CwiseBinaryOp.h

class Matrix
{
   
    // ...
    template<int LoadMode>
    inline PacketScalar packet(int index) const
    {
   
      return internal::ploadt<Scalar, LoadMode>(m_storage.data() + index);
    }
};

让我们来看一下GenericPacketMath.hinternal::ploadtsrc/Core/arch/SSE/PacketMath.h中的internal::pload的定义。它们与internal::pstore非常相似。

让我们回到CwiseBinaryOp::packet()。一旦来自向量vw的数据包被返回,这个函数会做什么?它在这些数据包上调用m_functor.packetOp()。那么m_functor是什么?在这里,我们必须记住我们正在处理什么样的CwiseBinaryOp的特定模板特化:

CwiseBinaryOp<internal::scalar_sum_op<float>, VectorXf, VectorXf>

m_functor是一个空类internal::scalar_sum_op<float>的对象。如上所述,不必担心为什么我们需要构造这个空类的对象——这是一个实现细节,重点是其他一些函数对象需要存储成员数据。

internal::scalar_sum_opsrc/Core/Functors.h 中定义:

template<typename Scalar> struct internal::scalar_sum_op EIGEN_EMPTY_STRUCT {
   
    inline const Scalar operator() (const Scalar& a, const Scalar& b) const 
    {
    return a + b; }
    template<typename PacketScalar>
    inline const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const
    {
    return internal::padd(a,b); }
};

正如你所看到的,packetOp()函数所做的就是调用internal::padd对这两个数据包进行求和。下面是src/Core/arch/SSE/PacketMath.hinternal::padd的定义:

template<> inline __m128  internal::padd(const __m128&  a, const __m128&  b) 
{
    return _mm_add_ps(a,b); }

在这里,_mm_add_ps是一个SSE特有的内部函数,表示单个SSE指令。

总结一下,如下循环:

for(int index = alignedStart; index < alignedEnd; index += packetSize)
{
   
    dst.template copyPacket<Derived2, Aligned, internal::assign_traits<Derived1,Derived2>::SrcAlignment>(index, src);
}

已经被编译成以下代码:对于 index 从0到11(= 48/4 - 1)的循环,使用两个__mm_load_ps SSE指令从向量v和向量w中读取第i个数据包(每个数据包包含4个浮点数),然后使用__mm_add_ps指令将它们相加,最后使用__mm_store_ps指令存储结果。

剩下的第二个循环处理最后几个(这里是最后2个)系数:

for(int index = alignedEnd; index < size; index++)
    dst.copyCoeff(index, src);

它的工作方式与我们刚刚解释的方式类似,只是更简单,因为这里没有涉及SSE矢量化。copyPacket() 变成了copyCoeff()packet()变成了coeff()writePacket()变成了coeffRef()。如果你看到这里,你可能可以自己理解这部分。

我们看到,在编译过程中,所有Eigen的C++抽象都消失了,我们确实可以精确控制所发出的汇编指令。这就是C++的美妙之处!由于我们对发出的汇编指令有如此精确的控制,但选择正确指令的逻辑如此复杂,因此我们可以说Eigen确实像一个优化编译器。如果你愿意,你也可以说Eigen的行为像一个编译器的脚本。从某种意义上说,C++模板元编程是编写编译器脚本的过程 - 已经证明这种脚本语言是图灵完备的。详情请参见维基百科

相关文章
|
存储 编译器
[Eigen中文文档] 深入了解 Eigen - 类层次结构
本页面介绍了Eigen类层次结构中 Core 类的设计及其相互关系。一般用户可能不需要关注这些细节,但对于高级用户和Eigen开发人员可能会有用。
291 0
|
C++
[Eigen中文文档] 按值将Eigen对象传递给函数
对于 Eigen,这一点更为重要:按值传递固定大小的可向量化 Eigen 对象不仅效率低下,而且可能是非法的或使程序崩溃! 原因是这些 Eigen 对象具有对齐修饰符,在按值传递时会不遵守这些修饰符。
165 0
|
存储 编译器
|
存储 C语言 C++
|
存储 算法 NoSQL
[Eigen中文文档] 稀疏矩阵操作
在许多应用中(例如,有限元方法),通常要处理非常大的矩阵,其中只有少数系数不为零。在这种情况下,可以通过使用仅存储非零系数的特殊表示来减少内存消耗并提高性能。这样的矩阵称为稀疏矩阵。
463 0
|
存储 缓存
[Eigen中文文档] 深入了解 Eigen - 惰性求值与混叠(Aliasing)
Eigen具有智能的编译时机制,可以实现惰性求值并在适当的情况下删除临时变量。它会自动处理大多数情况下的混叠问题,例如矩阵乘积。自动行为可以通过使用MatrixBase::eval()和MatrixBase::noalias()方法手动覆盖。
299 0
|
存储 索引
[Eigen中文文档] 扩展/自定义Eigen(三)
本页面针对非常高级的用户,他们不害怕处理一些Eigen的内部细节。在大多数情况下,可以通过使用自定义一元或二元函数避免使用自定义表达式,而极其复杂的矩阵操作可以通过零元函数(nullary-expressions)来实现,如前一页所述。 本页面通过示例介绍了如何在Eigen中实现新的轻量级表达式类型。它由三个部分组成:表达式类型本身、包含有关表达式编译时信息的特性类和评估器类,用于将表达式评估为矩阵。
135 1
|
存储 NoSQL API
[Eigen中文文档] Matrix类
在Eigen中,所有矩阵和向量都是Matrix模板类的对象。向量只是行数或者列数为1的特殊矩阵。
399 1
[Eigen中文文档] 在 CMake 项目中使用 Eigen
Eigen提供了CMake(CMake 3.0或更高版本)支持,使得该库可以轻松地在CMake项目中使用。
702 1
|
并行计算 算法 安全
[Eigen中文文档] Eigen 和多线程
某些 Eigen 算法可以利用硬件中存在的多个内核。
438 0