|
|
不知道为何上篇没阅读量，可能是写的简单了点。但其实是照顾大家，先从简单的理解，然后再来看本篇或uBLAS或论文中的表达式模板实现才心里有数。
如前所述，在上篇中我们已经实现了对加法的表达式模板计算。
然而缺点：极不直观，需要显示转化为最后的值类型；编码实现丑陋，每种运算符代码大量重复。
我们的基本标准依然是Eigen，同为CPP实现，没道理它能做的功能我做不出来。而在Eigen中，惰性求值表达式是隐藏于普通表达式中的，完美合并在一体。因此这也需要我们对编码表达式的模板类下更多功夫。
标量包装
我们的惰性求值实质上是一种对矩阵中相同下标元素的运算，因此matlab中的逐元操作应该能迁移至惰性求值中。但看过上篇实现的都知道，求值最后会调用[]操作符，而标量就是个数，没法弄。
所以我们需要对标量进行包装，让它支持惰性求值中的[]操作：
template <typename T>
struct expr_scalar
{
private:
T const &s;
public:
constexpr expr_scalar(T const &v)
: s(v)
{
}
constexpr T const &operator[](std::size_t) const
{
return s;
}
constexpr std::size_t size() const
{
return 0;
}
};虽然我们尽量让所有功能在编译器完成，但遗憾的是包装标量会有一个引用的额外开销。
有了包装标量，它让表达式中的算术类型也能参与模板运算。
求值表达式设计
此处的设计融合了STL中的valarray以及uBLAS中的部分模板技术。
再已经有了上篇文章中MatrixSum类体感以后，我们考虑一种更周全、更抽象的设计。它是一种能在模板参数中储存多元运算符的类，它的主要操作就说求值操作符[]。而具体的矩阵加减和逐元都能从这个基类中继承而来：
template <typename T>
class expr
{
public:
const T &self() const { return static_cast<const T &>(*this); }
T &self() { return static_cast<T &>(*this); }
protected:
explicit expr(){};
constexpr size_t size() { return self().size_impl(); }
auto operator[](size_t idx) const { return self().at_impl(idx); }
auto operator()() const { return self()(); };
};
这里明显不同于前篇文章的MatrixSum，我们实现了CPP的隐式类型转换操作符operator T()，这样能让我们在模板与值类型之间隐式的转换。这也是实现uBLAS或者Eigen那种惰性求值融合进普通表达式的关键。当然，如此一来，只要发生了隐式类型转换就会默认求值，到也不如前篇中显示指定那么明显了。
这里我还用到了CRTP避免虚继承调用开销。将操作符或说算符抽象出求值表达式基类好处在于：在对多种算符作惰性求值时，大规模作简单化代码编写。
二元操作符
二元操作符继承于求值表达式基类。具体内容和前篇MatrixSum差不多：
template <typename Ops, typename lExpr, typename rExpr>
class biops : public expr<biops<Ops, lExpr, rExpr>>
{
public:
using base_type = expr<biops<Ops, lExpr, rExpr>>;
using base_type::size;
using base_type::operator[];
friend base_type;
explicit biops(const Ops &ops, const lExpr &lxpr, const rExpr &rxpr)
: m_ops(ops), m_lxpr(lxpr), m_rxpr(rxpr){};
constexpr size_t size_impl() { return gl_get_more(m_lxpr.size(), m_rxpr.size()); };
auto at_impl(size_t idx) const { return m_ops(m_lxpr[idx], m_rxpr[idx]); };
template <typename T>
operator T()
{
T res{};
for (size_t idx = 0; idx < res.size(); ++idx)
{
res[idx] = (*this)[idx];
}
return res;
}
template <typename T, disable_arith_type_t<T> * = nullptr>
auto operator+(const T &rhs)
{
return biops<expr_plus_t, biops<Ops, lExpr, rExpr>, T>(expr_plus, *this, rhs);
}
template <typename T, enable_arith_type_t<T> * = nullptr>
auto operator+(const T &rhs)
{
return biops<expr_plus_t, biops<Ops, lExpr, rExpr>, expr_scalar<T>>(expr_plus, *this, rhs);
}
template <typename T, disable_arith_type_t<T> * = nullptr>
auto operator-(const T &rhs)
{
return biops<expr_minus_t, biops<Ops, lExpr, rExpr>, T>(expr_minus, *this, rhs);
}
template <typename T, enable_arith_type_t<T> * = nullptr>
auto operator-(const T &rhs)
{
return biops<expr_minus_t, biops<Ops, lExpr, rExpr>, expr_scalar<T>>(expr_minus, *this, rhs);
}
template <typename T, enable_arith_type_t<T> * = nullptr>
auto operator*(const T &rhs)
{
return biops<expr_mul_t, biops<Ops, lExpr, rExpr>, expr_scalar<T>>(expr_mul, *this, rhs);
}
template <typename T, enable_arith_type_t<T> * = nullptr>
auto operator/(const T &rhs)
{
return biops<expr_div_t, biops<Ops, lExpr, rExpr>, expr_scalar<T>>(expr_div, *this, rhs);
}
private:
Ops m_ops;
lExpr m_lxpr;
rExpr m_rxpr;
};值得注意的是，我把加减乘除这些算子具体化成了仿函数：
struct expr_plus_t
{
constexpr explicit expr_plus_t() = default;
template <typename LType, typename RType>
auto operator()(const LType &lhs, const RType &rhs) const
{
return lhs + rhs;
}
};
constexpr expr_plus_t expr_plus{};这么做的好处在于，我不用为每标量包装类/求值表达式类专门写几次加减乘除。
Matrix接口
前面的部分已经完成了求值表达式模板的功能，现在需要将它接入Matrix类。这里有两种方法：一是让它继承expr，顶掉自己以前的加减和逐元，让矩阵本身成为表达式模板的零阶运算数，而biops则是高阶的运算数；二是再次继承expr，生成一个零阶运算数，并直接在Matrix的加减和逐元中返回由biops包装的零阶对象启动求值。
我用第二种方式，因为不想让Matrix的继承关系太复杂。包装的零阶计算对象为：
template <typename T>
class expr_result : expr<expr_result<T>>
{
public:
using base_type = expr<expr_result<T>>;
using base_type::size;
using base_type::operator[];
friend base_type;
explicit expr_result(const T &val) : value(val) {}
size_t size_impl() const { return value.size(); };
auto at_impl(size_t idx) const { return value[idx]; };
decltype(auto) operator()() const { return (value); }
private:
const T &value;
};
它仅仅是数据对象的包裹。实际上你不创建这个类型已有的一样能计算，但进行二元运算时，会将数据对象拷贝一份，这与我们做模板计算求值目标正好相反，所以这种对数据成员的间接引用是必须的。
接下来一个麻烦在于二元运算符的实现。当然，你可以使用：
template<typename lExpr, typename rExpr>
auto operator+(const lExpr &lhs, const rExpr &rhs)
{ return binary_ops<vec_plus_t,lExpr,rExpr>(vec_plus,lhs,rhs); }这样的全局重载为每一个计算式模板生成运算函数，但这样毫无疑问污染了全局空间。当引入更多第三方库时，影响是无法估计的。
因此，我们只能将运算符实现为expr子类以及Matrix的成员函数。以加法为例，我们需要分别考虑：矩阵相加，矩阵标量相加，矩阵模板表达式相加，标量模板表达式相加等四种情况。
对于矩阵与矩阵和标量运算，在矩阵内部即可处理：
template <typename T, enable_arith_type_t<T> * = nullptr>
auto operator+(const T &other) const
{
using result_t = details::expr_result<Matrix>;
using result_s = details::expr_result<details::expr_scalar<T>>;
return details::biops<details::expr_plus_t, result_t, result_s>(details::expr_plus, result_t(*this), result_s(other));
}
template <typename T, disable_arith_type_t<T> * = nullptr>
auto operator+(const T &other) const
{
using result_t = details::expr_result<Matrix>;
return details::biops<details::expr_plus_t, result_t, T>(details::expr_plus, result_t(*this), other);
}
auto operator+(const Matrix &other) const
{
using result_t = details::expr_result<Matrix>;
return details::biops<details::expr_plus_t, result_t, result_t>(details::expr_plus, result_t(*this), result_t(other));
}使用SFINAE将输入T划分为三类处理，同理其它运算符。在expr类中再次操作：
template <typename T, disable_arith_type_t<T> * = nullptr, typename T::expr_type>
auto operator+(const T &rhs)
{
return biops<expr_plus_t, biops<Ops, lExpr, rExpr>, T>(expr_plus, *this, rhs);
}
template <typename T, disable_arith_type_t<T> * = nullptr>
auto operator+(const T &rhs)
{
using result_t = expr_result<T>;
return biops<expr_plus_t, biops<Ops, lExpr, rExpr>, result_t>(expr_plus, *this, result_t(rhs));
}
template <typename T, enable_arith_type_t<T> * = nullptr>
auto operator+(const T &rhs)
{
using result_s = expr_result<expr_scalar<T>>;
return biops<expr_plus_t, biops<Ops, lExpr, rExpr>, result_s>(expr_plus, *this, result_s(rhs));
}
我在expr的基类中派发了一个标签expr_type，方便我辨识类型做SFINAE。这一部分实现是最为繁琐：为了保证符号限制在矩阵库内部，需要显示的，精细的处理每次运算的操作数。
BUG修复
上面的惰性求值模板当中略有BUG：这个模板类在标量参与运算时，会创建右值的标量包装类，然后进行二元操作后，求值表达式会储存这个模板右值的引用，在进行延时求值之后，会导致索引不存在内存然后Coredump。
解决问题的办法是通过类型萃取让求值表达式储存表达式类型的引用而对标量才用值拷贝。这里我实现一个简单的类型萃取：
template <typename T>
struct expr_traits
{
using ExprRef = T const &;
};
template <typename T>
struct expr_traits<expr_scalar<T>>
{
using ExprRef = expr_scalar<T>;
};主模板中声明一个类型标签，偏特化中更改类型标签的名称。这种技术在CPP模板编程中使用的极多。
有了类型萃取模板，在零阶递归的求值表达式中可以这样声明变量：
private:
typename expr_traits<T>::ExprRef value;这样就能根据输入参数自动在值类型与引用类型之间转换。
Benchmark
按上文改造接口以后，我们的惰性求值功能基本和主流的CPP线代库看起来一样了：
Matrix<4, 4> x = {1, 2, 3, 4, 5, 6, 7, 8};
Matrix<4, 4> y = x.T();
Matrix<4, 4> z = x * 2 + y * 2 + 3.3 + x * y;
PRINT_SINGLE_ELEMENTS(z, &#34;2*(x + x^T) + 3.3 +x*x^T = &#34;);
2*(x + x^T) + 3.3 +x*x^T = Matrix<4,4>:
33.3 49.3 47.3 55.3
49.3 67.3 65.3 75.3
47.3 65.3 61.3 71.3
55.3 75.3 71.3 83.3在dtor中加入打印，观察是否有惰性求值过程：
2*(x + x^T) + 3.3 +x*x^T = Matrix<4,4>:
33.3 49.3 47.3 55.3
49.3 67.3 65.3 75.3
47.3 65.3 61.3 71.3
55.3 75.3 71.3 83.3
dtor!
dtor!
dtor!可以发现，仅仅发生三次析构。没有任何多余对象。
依旧以Eigen做基准，对算式 2x + 2y + Mat(3.3) + xy+x+y+x+y 进行计算：
2022-11-19T02:33:44+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1609.46 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x6)
L1 Instruction 32 KiB (x6)
L2 Unified 256 KiB (x6)
L3 Unified 12288 KiB (x1)
-------------------------------------------------------------
Benchmark Time CPU Iterations
-------------------------------------------------------------
BM_MatrixMul 3629 ns 3498 ns 263529
BM_MatrixMulEigen 1556 ns 1535 ns 448000
BM_MatrixInv 2225 ns 2093 ns 298667
BM_MatrixEigenInv 3050 ns 3069 ns 224000
BM_MatrixExpr 681 ns 672 ns 1000000
BM_MatrixEigenExpr 645 ns 656 ns 1120000再对比没有加入惰性求值之前：
2022-11-19T02:31:48+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1617.86 MHz CPU s)
CPU Caches:
L1 Data 32 KiB (x6)
L1 Instruction 32 KiB (x6)
L2 Unified 256 KiB (x6)
L3 Unified 12288 KiB (x1)
-------------------------------------------------------------
Benchmark Time CPU Iterations
-------------------------------------------------------------
BM_MatrixMul 3671 ns 3599 ns 186667
BM_MatrixMulEigen 1241 ns 1256 ns 560000
BM_MatrixInv 2141 ns 2148 ns 320000
BM_MatrixEigenInv 3068 ns 3115 ns 235789
BM_MatrixExpr 1587 ns 1604 ns 448000
BM_MatrixEigenExpr 636 ns 628 ns 1120000加入惰性求值以后提升了大约一半的时间，基本和Eigen没有区别。并且这个时间随着运算对象增多以及加大都更有意义。
最新代码可在github获取 |
|
发表于 2022-12-3 18:57:05