处理 Eigen 中两种存储布局的固定大小矩阵的连续向量

Question

外部库为我提供了 double 的原始指针，我想将其映射到 Eigen 类型。原始数组在逻辑上是一个大的有序集合，由大小相同的小密集固定大小矩阵组成。主要问题是小密集矩阵可能按行优先或列优先排序，我想同时容纳它们。

我目前的做法如下。请注意，一个小的固定大小块（在块数组中）的所有条目在内存中需要是连续的。

template<int bs, class Mattype>
void block_operation(double *const vals, const int numblocks)
{
    Eigen::Map<Mattype> mappedvals(vals, 
                Mattype::IsRowMajor ? numblocks*bs : bs,
                Mattype::IsRowMajor ? bs : numblocks*bs
            );

    for(int i = 0; i < numblocks; i++)
      if(Mattype::isRowMajor)
          mappedvals.template block<bs,bs>(i*bs,0) = block_operation_rowmajor(mappedvals);
      else
          mappedvals.template block<bs,bs>(0,i*bs) = block_operation_colmajor(mappedvals);
}

调用函数首先计算出 Mattype（在 2 个选项中），然后使用正确的模板参数调用上述函数。

因此我所有的算法都需要写两次，我的代码中穿插着这些布局检查。有没有办法以布局不可知的方式做到这一点？请记住，此代码需要尽可能快。

理想情况下，我会 Map 数据一次并将其用于所有需要的操作。但是，我能想出的唯一解决方案是每当我需要访问该块时，为每个小块调用一次 Map 构造函数。

template<int bs, StorageOptions layout>
inline Map<Matrix<double,bs,bs,layout>> extractBlock(double *const vals, 
                                                     const int bindex)
{
  return Map<Matrix<double,bs,bs,layout>>(vals+bindex*bs*bs);
}

这个函数是否会被优化到什么都没有（通过像 GCC 7.3 或 Intel 2017 这样的现代编译器在 -std=c++14 -O3 下），或者我每次调用这个函数时都会付出一点代价（每个函数一次）块，还有很多小块）？有更好的方法吗？

Answer 1

你的 extractBlock 很好，一个更简单但有点丑陋的解决方案是在 block_operation 的开头使用重新解释转换：

using BlockType = Matrix<double,bs,bs,layout|DontAlign>;
BlockType* blocks = reinterpret_cast<BlockType*>(vals);

for(int i...)
  block[i] = ...;

这仅适用于固定大小的矩阵。另请注意 DontAlign 这很重要，除非您可以保证 vals 根据 AVX 和 bs 的存在在 16 甚至 32 字节上对齐...所以只需使用 DontAlign！

处理 Eigen 中两种存储布局的固定大小矩阵的连续向量

Dealing with a contiguous vector of fixed-size matrices for both storage layouts in Eigen

c++

eigen

eigen3