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

文档总目录

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

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

u = v + w;

这里调用了运算符 =，向量 u 是 VectorXf 类的对象，即 Matrix。 src/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());
}

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

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

这里，Derived是VectorXf（因为u是一个VectorXf），而OtherDerived是CwiseBinaryOp。具体来说，如前面所述，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宏，CwiseBinaryOp的Flags成员从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个参数自动确定。

这两个参数Vectorization和Unrolling由一个辅助类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个，packetSize是4。因此，我们将不得不逐个复制最后2个系数，而不是通过数据包。在这里，alignedEnd为48。

接下来是实际的循环。

首先是向量化的部分：前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);
}

在这里，StoreMode是Aligned，表示我们正在进行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::pstore和internal::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.h中internal::ploadt和src/Core/arch/SSE/PacketMath.h中的internal::pload的定义。它们与internal::pstore非常相似。

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

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

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

internal::scalar_sum_op在 src/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.h中internal::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++模板元编程是编写编译器脚本的过程 - 已经证明这种脚本语言是图灵完备的。详情请参见维基百科。