Почему компиляторы C ++ не делают лучшее постоянное свертывание?
Я исследую способы ускорить большую часть кода C ++, который имеет автоматические производные для вычисления якобиан. Это включает выполнение некоторого объема работы с фактическими остатками, но большая часть работы (основанная на профилированном времени выполнения) заключается в вычислении якобиан.
Это удивило меня, так как большинство якобианов распространяются вперед от 0 до 1, поэтому объем работы должен быть в 2-4 раза больше функции, а не в 10-12 раз. Чтобы смоделировать, на что похожа большая часть якобианской работы, я сделал супер минимальный пример с простым точечным произведением (вместо sin, cos, sqrt и других, которые были бы в реальной ситуации), что компилятор должен иметь возможность оптимизировать до единого возвращаемого значения:
#include <Eigen/Core>
#include <Eigen/Geometry>
using Array12d = Eigen::Matrix<double,12,1>;
double testReturnFirstDot(const Array12d& b)
{
Array12d a;
a.array() = 0.;
a(0) = 1.;
return a.dot(b);
}
Который должен быть таким же, как
double testReturnFirst(const Array12d& b)
{
return b(0);
}
Я был разочарован, обнаружив, что без включенной быстрой математики ни GCC 8.2, ни Clang 6, ни MSVC 19 не смогли вообще выполнить какую-либо оптимизацию для наивного точечного продукта с матрицей, полной 0. Даже с фаст-математикой (https://godbolt.org/z/GvPXFy) оптимизации в GCC и Clang очень плохие (все еще включают умножения и добавления), а MSVC вообще не выполняет никаких оптимизаций.
У меня нет опыта работы с компиляторами, но есть ли причина для этого? Я вполне уверен, что в значительной части научных вычислений возможность сделать лучшее постоянное распространение / свертывание сделает более очевидными оптимизацию, даже если само постоянное сгибание не приведет к ускорению.
Хотя меня интересуют объяснения того, почему это не делается на стороне компилятора, меня также интересует, что я могу сделать с практической стороны, чтобы ускорить создание собственного кода при работе с шаблонами такого типа.
Решение
Это потому, что Eigen явно векторизовал ваш код как 3 vmulpd, 2 vaddpd и 1 горизонтальное сокращение в оставшихся 4-х компонентных регистрах (это предполагает AVX, только с SSE вы получите 6 mulpd и 5 addpd). С -ffast-math GCC и clang разрешено удалять последние 2 vmulpd и vaddpd (и это то, что они делают), но они не могут реально заменить оставшиеся vmulpd и горизонтальное сокращение, которые были явно сгенерированы Eigen.
Так что, если вы отключите явную векторизацию Эйгена, определив EIGEN_DONT_VECTORIZE? Тогда вы получите то, что ожидали (https://godbolt.org/z/UQsoeH) но другие части кода могут стать намного медленнее.
Если вы хотите локально отключить явную векторизацию и не боитесь возиться с внутренним Eigen, вы можете ввести DontVectorize возможность Matrix и отключить векторизацию, специализируясь traits<> за это Matrix тип:
static const int DontVectorize = 0x80000000;
namespace Eigen {
namespace internal {
template<typename _Scalar, int _Rows, int _Cols, int _MaxRows, int _MaxCols>
struct traits<Matrix<_Scalar, _Rows, _Cols, DontVectorize, _MaxRows, _MaxCols> >
: traits<Matrix<_Scalar, _Rows, _Cols> >
{
typedef traits<Matrix<_Scalar, _Rows, _Cols> > Base;
enum {
EvaluatorFlags = Base::EvaluatorFlags & ~PacketAccessBit
};
};
}
}
using ArrayS12d = Eigen::Matrix<double,12,1,DontVectorize>;
Полный пример там: https://godbolt.org/z/bOEyzv
Другие решения
Я был разочарован, обнаружив, что без включенной быстрой математики ни GCC 8.2, ни Clang 6, ни MSVC 19 не смогли вообще выполнить какую-либо оптимизацию для наивного точечного продукта с матрицей, полной 0.
У них нет другого выбора, к сожалению. Поскольку IEEE-числа имеют подписанные нули, добавляя 0.0 не является идентификационной операцией:
-0.0 + 0.0 = 0.0 // Not -0.0!
Аналогично, умножение на ноль не всегда дает ноль:
0.0 * Infinity = NaN // Not 0.0!
Таким образом, компиляторы просто не могут выполнять эти постоянные сгибы в точечном произведении, сохраняя соответствие стандарту IEEE с плавающей запятой — насколько они знают, ваши входные данные могут содержать подписанные нули и / или бесконечности.
Вам придется использовать -ffast-math чтобы получить эти складки, но это может иметь нежелательные последствия. Вы можете получить более детальный контроль с определенными флагами (от http://gcc.gnu.org/wiki/FloatingPointMath). Согласно приведенному выше объяснению, добавление следующих двух флагов должно позволить постоянное свертывание:
-ffinite-math-only, -fno-signed-zeros
Действительно, вы получаете ту же сборку, что и с -ffast-math сюда: https://godbolt.org/z/vGULLA. Вы отказываетесь только от подписанных нулей (вероятно, не относящихся к делу), NaN и бесконечностей. Предположительно, если бы вы все еще производили их в своем коде, вы бы получили неопределенное поведение, поэтому взвесите ваши варианты.
Что касается того, почему ваш пример не оптимизирован лучше даже с -ffast-mathЭто на Эйгене. Предположительно, они имеют векторизацию в своих матричных операциях, которые компиляторам гораздо сложнее увидеть. Простой цикл должным образом оптимизирован с этими параметрами: https://godbolt.org/z/OppEhY
Один из способов заставить компилятор оптимизировать умножения на 0 и 1 — это вручную развернуть цикл. Для простоты давайте использовать
#include <array>
#include <cstddef>
constexpr std::size_t n = 12;
using Array = std::array<double, n>;
Тогда мы можем реализовать простой dot функция с использованием выражений складывания (или рекурсии, если они недоступны):
<utility>
template<std::size_t... is>
double dot(const Array& x, const Array& y, std::index_sequence<is...>)
{
return ((x[is] * y[is]) + ...);
}
double dot(const Array& x, const Array& y)
{
return dot(x, y, std::make_index_sequence<n>{});
}
Теперь давайте посмотрим на вашу функцию
double test(const Array& b)
{
const Array a{1}; // = {1, 0, ...}
return dot(a, b);
}
С -ffast-math GCC 8,2 производит:
test(std::array<double, 12ul> const&):
movsd xmm0, QWORD PTR [rdi]
ret
Clang 6.0.0 идет в том же духе:
test(std::array<double, 12ul> const&): # @test(std::array<double, 12ul> const&)
movsd xmm0, qword ptr [rdi] # xmm0 = mem[0],zero
ret
Например, для
double test(const Array& b)
{
const Array a{1, 1}; // = {1, 1, 0...}
return dot(a, b);
}
мы получаем
test(std::array<double, 12ul> const&):
movsd xmm0, QWORD PTR [rdi]
addsd xmm0, QWORD PTR [rdi+8]
ret
Дополнение. Clang разворачивает for (std::size_t i = 0; i < n; ++i) ... цикл без всех этих уловок выражений сгиба, gcc не делает и нуждается в некоторой помощи.