添加新的表达式类型
英文原文(Adding a new expression type)
本页面针对非常高级的用户，他们不害怕处理一些Eigen的内部细节。在大多数情况下，可以通过使用自定义一元或二元函数避免使用自定义表达式，而极其复杂的矩阵操作可以通过零元函数(nullary-expressions)来实现，如前一页所述。
本页面通过示例介绍了如何在Eigen中实现新的轻量级表达式类型。它由三个部分组成：表达式类型本身、包含有关表达式编译时信息的特性类和评估器类，用于将表达式评估为矩阵。
设定
循环矩阵是每列与左边的列相同的矩阵，只不过它循环向下移位。例如，这是一个 4×4 循环矩阵：
$$ \begin{bmatrix} 1&8&4&2 \\ 2&1&8&4 \\ 4&2&1&8 \\ 8&4&2&1 \end{bmatrix} $$
循环矩阵由其第一列唯一确定。我们希望编写一个函数 makeCirculant，在给定第一列的情况下，它返回一个表示循环矩阵的表达式。
为简单起见，我们将 makeCirculant 函数限制为稠密矩阵。也允许数组或稀疏矩阵可能是有意义的，但我们不会在这里这样做。我们也不希望支持矢量化。
开始
我们将分部分介绍实现makeCirculant函数的文件。首先，我们包含适当的头文件并向前声明我们将称之为Circulant的表达式类。makeCirculant函数将返回此类型的对象。实际上，类Circulant是一个类模板；模板参数ArgType是传递给makeCirculant函数的向量的类型。
#include <Eigen/Core>
#include <iostream>
template <class ArgType> class Circulant;
Traits类
对于每个表达式类 X，在 Eigen::internal 命名空间中应该有一个特征类 Traits<X>，其中包含有关 X 的信息（称为编译时）。
如“设定”中所解释的，我们设计了Circulant表达式类来引用密集矩阵。循环矩阵的条目与传递给makeCirculant函数的向量的条目具有相同的类型。用于索引条目的类型也相同。为简单起见，我们仅返回列主矩阵。最后，循环矩阵是一个方阵（行数等于列数），行数等于传递给makeCirculant函数的列向量的行数。如果这是一个动态大小向量，则循环矩阵的大小在编译时是未知的。
namespace Eigen {
namespace internal {
template <class ArgType>
struct traits<Circulant<ArgType> >
{
typedef Eigen::Dense StorageKind;
typedef Eigen::MatrixXpr XprKind;
typedef typename ArgType::StorageIndex StorageIndex;
typedef typename ArgType::Scalar Scalar;
enum {
Flags = Eigen::ColMajor,
RowsAtCompileTime = ArgType::RowsAtCompileTime,
ColsAtCompileTime = ArgType::RowsAtCompileTime,
MaxRowsAtCompileTime = ArgType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = ArgType::MaxRowsAtCompileTime
};
};
}
}
表达式类
接下来的步骤是定义表达式类本身。要继承自MatrixBase，以公开密集矩阵的接口。在构造函数中，我们检查是否传递了列向量（参见断言），并将用于构建循环矩阵的向量存储在成员变量m_arg中。最后，表达式类应计算相应循环矩阵的大小。如上所述，这是一个方阵，其列数与用于构造矩阵的向量一样多。
template <class ArgType>
class Circulant : public Eigen::MatrixBase<Circulant<ArgType> >
{
public:
Circulant(const ArgType& arg)
: m_arg(arg)
{
EIGEN_STATIC_ASSERT(ArgType::ColsAtCompileTime == 1,
YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
}
typedef typename Eigen::internal::ref_selector<Circulant>::type Nested;
typedef Eigen::Index Index;
Index rows() const {
return m_arg.rows(); }
Index cols() const {
return m_arg.rows(); }
typedef typename Eigen::internal::ref_selector<ArgType>::type ArgTypeNested;
ArgTypeNested m_arg;
};
评估器
最后一个大片段实现了Circulant表达式的评估器。评估器计算循环矩阵的条目。这是在.coeff()成员函数中完成的。通过找到构成循环矩阵的向量的相应条目来计算这些条目。当循环矩阵是从给定复杂表达式的向量构建时，获取此条目可能实际上并不容易，因此我们使用对应于该向量的评估器。
CoeffReadCost 常量记录计算循环矩阵一项的成本；我们忽略索引计算，并认为这与计算构建循环矩阵的向量条目的成本相同。
在构造函数中，我们保存定义循环矩阵的列向量的求值器。还保存该向量的大小；请记住，我们可以查询表达式对象来查找大小，但不能查询求值器。
namespace Eigen {
namespace internal {
template<typename ArgType>
struct evaluator<Circulant<ArgType> >
: evaluator_base<Circulant<ArgType> >
{
typedef Circulant<ArgType> XprType;
typedef typename nested_eval<ArgType, XprType::ColsAtCompileTime>::type ArgTypeNested;
typedef remove_all_t<ArgTypeNested> ArgTypeNestedCleaned;
typedef typename XprType::CoeffReturnType CoeffReturnType;
enum {
CoeffReadCost = evaluator<ArgTypeNestedCleaned>::CoeffReadCost,
Flags = Eigen::ColMajor
};
evaluator(const XprType& xpr)
: m_argImpl(xpr.m_arg), m_rows(xpr.rows())
{
}
CoeffReturnType coeff(Index row, Index col) const
{
Index index = row - col;
if (index < 0) index += m_rows;
return m_argImpl.coeff(index);
}
evaluator<ArgTypeNestedCleaned> m_argImpl;
const Index m_rows;
};
}
}
入口
template <class ArgType>
Circulant<ArgType> makeCirculant(const Eigen::MatrixBase<ArgType>& arg)
{
return Circulant<ArgType>(arg.derived());
}
完成这些之后，makeCirculant 函数就非常简单了。它只是创建一个表达式对象并返回它。
一个简单用于测试的 main 函数
最后，一个简短的主函数展示了如何调用 makeCirculant 函数。
int main()
{
Eigen::VectorXd vec(4);
vec << 1, 2, 4, 8;
Eigen::MatrixXd mat;
mat = makeCirculant(vec);
std::cout << mat << std::endl;
}
如果组合所有代码片段，则会产生以下输出，表明程序按预期工作：
1 8 4 2
2 1 8 4
4 2 1 8
8 4 2 1
