Eigen：对于这个模板表达式，为什么 Map 比 Vector3d 慢？

Question

我在 x、y、z 模式的 std::vector<double> 中有一个点云，以及一个 std::vector<int> 索引，其中连续整数的每个三元组是一张脸的连通性。基本上是一个简单的三角形网格数据结构。

我必须计算所有面孔的面积，我正在对几种方法进行基准测试：

我可以像这样将数据块包装在 Eigen::Map<const Eigen::Vector3d> 中：

static void face_areas_eigenmap(const std::vector<double>& V,
                                const std::vector<int>& F,
                                std::vector<double>& FA) {
  // Number of faces is size / 3.
  for (auto f = 0; f < F.size() / 3; ++f) {
    // Get vertex indices of face f.
    auto v0 = F[f * 3];
    auto v1 = F[f * 3 + 1];
    auto v2 = F[f * 3 + 2];
    
    // View memory at each vertex position as a vector.
    Eigen::Map<const Eigen::Vector3d> x0{&V[v0 * 3]};
    Eigen::Map<const Eigen::Vector3d> x1{&V[v1 * 3]};
    Eigen::Map<const Eigen::Vector3d> x2{&V[v2 * 3]};
    
    // Compute and store face area.
    FA[f] = 0.5 * (x1 - x0).cross(x2 - x0).norm();
  }
}

或者我可以选择这样创建 Eigen::Vector3d：

static void face_areas_eigenvec(const std::vector<double>& V,
                                const std::vector<int>& F,
                                std::vector<double>& FA) {
  for (auto f = 0; f < F.size() / 3; ++f) {
    auto v0 = F[f * 3];
    auto v1 = F[f * 3 + 1];
    auto v2 = F[f * 3 + 2];
    
    // This is the only change, swap Map for Vector3d.
    Eigen::Vector3d x0{&V[v0 * 3]};
    Eigen::Vector3d x1{&V[v1 * 3]};
    Eigen::Vector3d x2{&V[v2 * 3]};

    FA[f] = 0.5 * (x1 - x0).cross(x2 - x0).norm();
  }
}

最后我也在考虑带有显式叉积和范数的硬编码版本：

static void face_areas_ptr(const std::vector<double>& V,
                           const std::vector<int>& F, std::vector<double>& FA) {
  for (auto f = 0; f < F.size() / 3; ++f) {
    const auto* x0 = &V[F[f * 3] * 3];
    const auto* x1 = &V[F[f * 3 + 1] * 3];
    const auto* x2 = &V[F[f * 3 + 2] * 3];

    std::array<double, 3> s0{x1[0] - x0[0], x1[1] - x0[1], x1[2] - x0[2]};
    std::array<double, 3> s1{x2[0] - x0[0], x2[1] - x0[1], x2[2] - x0[2]};

    std::array<double, 3> c{s0[1] * s1[2] - s0[2] * s1[1],
                            s0[2] * s1[0] - s0[0] * s1[2],
                            s0[0] * s1[1] - s0[1] * s1[0]};

    FA[f] = 0.5 * std::sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
  }
}

我已经对这些方法进行了基准测试，使用 Eigen::Map 的版本总是最慢的，尽管它与使用 Eigen::Vector3d 的版本做同样的事情，但我期望性能没有变化，因为地图是基本上是一个指针。

-----------------------------------------------------------------
Benchmark                       Time             CPU   Iterations
-----------------------------------------------------------------
BM_face_areas_eigenvec   59757936 ns     59758018 ns           11
BM_face_areas_ptr        58305018 ns     58304436 ns           11
BM_face_areas_eigenmap   62356850 ns     62354710 ns           10

我试过使用与指针版本相同的代码在地图版本中切换特征模板表达式：

std::array<double, 3> s0{x1[0] - x0[0], x1[1] - x0[1], x1[2] - x0[2]};
std::array<double, 3> s1{x2[0] - x0[0], x2[1] - x0[1], x2[2] - x0[2]};

std::array<double, 3> c{s0[1] * s1[2] - s0[2] * s1[1],
                        s0[2] * s1[0] - s0[0] * s1[2],
                        s0[0] * s1[1] - s0[1] * s1[0]};

FA[f] = 0.5 * std::sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);

神奇的是时间是可比较的：

-----------------------------------------------------------------
Benchmark                       Time             CPU   Iterations
-----------------------------------------------------------------
BM_face_areas_array      58967864 ns     58967891 ns           11
BM_face_areas_ptr        60034545 ns     60034682 ns           11
BM_face_areas_eigenmap   60382482 ns     60382027 ns           11

Eigen 表达式中 Eigen::Map 有什么问题需要注意吗？

Answer 1

查看编译器输出，似乎第二个版本通过将其中一些内存负载聚合到矢量负载中，使编译器发出更少的内存负载。 https://godbolt.org/z/qs38P41eh

cross 的 Eigen 代码不包含任何显式矢量化。这取决于编译器能否很好地处理它。而且由于您在表达式（减法）上调用 cross，编译器会过早放弃。基本上是编译器没有找到相同优化的错。

您的第三个代码与第二个代码的工作方式相同，因为编译器将减法（创建 s0 和 s1）识别为它可以执行矢量化的操作，从而生成等效代码。如果你这样做，你可以用 Eigen 达到同样的效果：

    Eigen::Map<const Eigen::Vector3d> x0{&V[v0 * 3]};
    Eigen::Map<const Eigen::Vector3d> x1{&V[v1 * 3]};
    Eigen::Map<const Eigen::Vector3d> x2{&V[v2 * 3]};
    
    Eigen::Vector3d s0 = x1 - x0;
    Eigen::Vector3d s1 = x2 - x0;

    // Compute and store face area.
    FA[f] = 0.5 * s0.cross(s1).norm();

Eigen：对于这个模板表达式，为什么 Map 比 Vector3d 慢？

Eigen: Why is Map slower than Vector3d for this template expression?

c++

linear-algebra

eigen

c++20

eigen3