C++矩阵寄算（二）：Fast than light！栈上运算

7 个月前 · 来自专栏 CPP11魔法技术

幽州狼

公无渡河，公竟渡河。坠河而死，其奈公何？

昨天摸完基本实现后，一直想能不能把小矩阵对象放栈上。因为栈上对象在Cache里，访存比堆上快得多。也幻想过STL的vector和string一样，有SSO小对象优化，但测了一番终究是没有。

（栈上对象访问一般都是在一级缓存，编译器直接栈指针做偏移就得到值了，是最快的访问方式。堆上变量需要访问栈上指针再通过这个指针寻址。并且堆上变量经常被换出缓存，容易触发缺页中断导致性能大幅下降。）

然后改用array做容器，结果测试随便爆栈。

虽然各种爆栈，但benchmark发现栈上对象运算速度大概是堆上10倍，比Eigen还快一倍以上。

这下没法忽视这个优化点了。首先思考能否通过Allocator，达到小对象栈上，大对象堆上效果。写好结果发现内存池下移动语义没法处理。

中间想过很多方法，CRTP特化派生类数据成员，虚基类做方法接口派生只有成员等等。。。

最后返璞归真：通过类模板偏特化，做一个小矩阵array容器，大矩阵vector容器数据成员基类；在派生类继承数据成员，实现各类线代方法，最后把派生类包装成标准STL容器。

偏特化实现堆栈内存区分

我们需要根据矩阵的行列数特化出两个类型，特化子类型只有一个数据成员，就是array或者vector。于是我们至少需要三个模板参数：行，列，是否小矩阵。

基类定义如下：

    constexpr size_t gl_sm(size_t A, size_t B)
        return A * B < 280 ? 1 : 0;
    template <size_t M, size_t N, std::size_t A = gl_sm(M, N)>
    class MatrixBase;

gl_sm是编译期表达式，用以判断大小。因为CPP里面大于符号和模板符号冲突。第三个参数用0表示大矩阵1表示小矩阵，然后根据矩阵维数赋予默认值。

然后做偏特化：

    template <std::size_t M, std::size_t N>
    class MatrixBase<M, N, 1>
    protected:
        MatrixBase() : m_data{} {}
        ~MatrixBase() = default;
        MatrixBase(const MatrixBase &) = default;
        MatrixBase(MatrixBase &&) = default;
        MatrixBase &operator=(const MatrixBase &) = default;
        MatrixBase &operator=(MatrixBase &&) = default;
    protected:
        std::array<double, M * N> m_data;
    template <std::size_t M, std::size_t N>
    class MatrixBase<M, N, 0>
        template <typename T, typename RT = void>
        using enable_arith_type_t = typename std::enable_if<std::is_arithmetic<T>::value, RT>::type;
    protected:
        MatrixBase() : m_data(M * N, 0.0) {}
        ~MatrixBase() = default;
        MatrixBase(const MatrixBase &) = default;
        MatrixBase(MatrixBase &&) = default;
        MatrixBase &operator=(const MatrixBase &) = default;
        MatrixBase &operator=(MatrixBase &&) = default;
    protected:
        std::vector<double> m_data;
    };

主要是一些构造函数，我略去了一部分从各类数组构造的部分。要注意CPP中array是标准容器中唯一默认初始化不为零的容器。vector需要初始化原因是提前分配空间，避免不断复制自己。

有这个基类以后，我们再用面向用户的矩阵类Matrix继承它。继承时候因为基类模板的第三参数已经自动填写，所以会根据矩阵大小自动实例化成相应vector或者array。我们也将基类的数据成员声明成一个名称，方便操作。对于vector和array有不同的接口时，就在基类中封装，保证Matrix调用到的始终是一个名字。这样就既实现根据矩阵大小在堆或栈上分配内存，又隐藏了细节。

benchmark

还是昨天的例子，依旧是和Eigen比较。测试项目是随机稠密矩阵 A*A^{T}

测试一：维数 A[500,15]*A^{T}[15,500] ：

这个测试中数组长度500*15超过预设280，矩阵A和其转置以及结果都在堆上存在，速度应该和昨天一样，约Eigen一半：

2022-11-14T21:38:42+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1613.97 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         2279471 ns      2264493 ns          345
BM_MatrixMulEigen    1331496 ns      1325335 ns          448
BM_MatrixInv            2386 ns         2372 ns       263529
BM_MatrixEigenInv       3281 ns         3296 ns       213333

测试二：维数 A[15,15]*A^{T}[15,15] ：

这个测试中数组长度15*15未超过预设280，矩阵A和其转置以及结果都在栈上存在，速度应该非常快：

2022-11-14T21:43:21+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1611.89 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             222 ns          214 ns      2986667
BM_MatrixMulEigen        609 ns          600 ns      1120000
BM_MatrixInv            2030 ns         1995 ns       344615
BM_MatrixEigenInv       3459 ns         3488 ns       179200

测试三：维数 A[20,12]*A^{T}[12,20] ：

两个栈上一个堆上，速度差不多以堆上运算为准：

2022-11-14T21:49:28+08:00
Running C:\Users\Desktop\matrix\build\Release\matrix_bench.exe
Run on (12 X 1629.7 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            3578 ns         3453 ns       203636
BM_MatrixMulEigen       1170 ns         1200 ns       560000
BM_MatrixInv            2062 ns         2086 ns       344615
BM_MatrixEigenInv       2860 ns         2849 ns       235789

测试可以发现，Eigen的运算速度随着矩阵大小波动小一些，合理猜测它实现了类似Allocator一样的内存池策略，所以不像我们快的时候比Eigen快一倍，慢的1/3情况。

堆栈分配的界限是根据栈大小来的，我以15*15为界限。在SE3群中，最多的就是4维6维运算。不担心爆栈崩溃该参数可以随意调整。

包装迭代器

这样改变类布局后留了个尾巴。以前类的唯一数据成员通过flat可以直接访问；现在则是根据类的实例化有vector和array两种情况。对外的迭代器接口不统一了。因此需要按STL要求自己写一个迭代器。

struct iterator
        using iterator_category = std::bidirectional_iterator_tag;
        using difference_type = std::ptrdiff_t;
        using value_type = double;
        using pointer = double *;
        using reference = double &;
        iterator(pointer ptr) : m_ptr(ptr) {}
        reference operator*() const
            return *m_ptr;
        pointer operator->()
            return m_ptr;
        iterator &operator++()
            m_ptr++;
            return *this;
        iterator operator++(int)
            iterator tmp = *this;
            ++(*this);
            return tmp;
        iterator &operator--()
            m_ptr--;
            return *this;
        iterator operator--(int)
            iterator tmp = *this;
            --(*this);
            return tmp;