Блог

У меня есть облако точек в a std::vector<double> в шаблоне x, y, z и a std::vector<int> индексов, где каждая тройка последовательных целых чисел является связностью грани. В основном простая структура данных треугольной сетки.

Мне нужно вычислить площади всех граней, и я сравниваю несколько методов:

Я могу обернуть фрагменты данных следующим Eigen::Map<const Eigen::Vector3d> образом:

 static void face_areas_eigenmap(const std::vector<double>amp; V,
                                const std::vector<int>amp; F,
                                std::vector<double>amp; 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{amp;V[v0 * 3]};
    Eigen::Map<const Eigen::Vector3d> x1{amp;V[v1 * 3]};
    Eigen::Map<const Eigen::Vector3d> x2{amp;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>amp; V,
                                const std::vector<int>amp; F,
                                std::vector<double>amp; 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{amp;V[v0 * 3]};
    Eigen::Vector3d x1{amp;V[v1 * 3]};
    Eigen::Vector3d x2{amp;V[v2 * 3]};

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

Наконец, я также рассматриваю жестко запрограммированную версию с явным перекрестным произведением и нормой:

 static void face_areas_ptr(const std::vector<double>amp; V,
                           const std::vector<int>amp; F, std::vector<double>amp; FA) {
  for (auto f = 0; f < F.size() / 3;   f) {
    const auto* x0 = amp;V[F[f * 3] * 3];
    const auto* x1 = amp;V[F[f * 3   1] * 3];
    const auto* x2 = amp;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 , я не ожидал никаких изменений в производительности, поскольку map — это в основном указатель.

 -----------------------------------------------------------------
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::Map собственными выражениями, о которых нужно знать?

Глядя на вывод компилятора, кажется, что вторая версия заставляет компилятор выделять меньше нагрузок на память, объединяя некоторые из них в векторные нагрузки. https://godbolt.org/z/qs38P41eh

Собственный код для cross не содержит никакой явной векторизации. Это зависит от того, хорошо ли компилятор справляется с этим. И поскольку вы вызываете cross для выражения (вычитания), компилятор сдается слишком рано. По сути, это ошибка компилятора в том, что он не нашел ту же оптимизацию.

Ваш третий код работает так же, как и второй, потому что компилятор распознает вычитание (создание s0 и s1) как то, что он может сделать векторизованным, что приводит к эквивалентному коду. Вы можете добиться того же с помощью Eigen, если сделаете это следующим образом:

     Eigen::Map<const Eigen::Vector3d> x0{amp;V[v0 * 3]};
    Eigen::Map<const Eigen::Vector3d> x1{amp;V[v1 * 3]};
    Eigen::Map<const Eigen::Vector3d> x2{amp;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();