jkf*_*ing 59 c++ compiler-construction automatic-differentiation eigen ceres-solver
我正在研究加速大部分C++代码的方法,这些代码具有用于计算jacobians的自动衍生产品.这涉及在实际残差中做一些工作,但大部分工作(基于异形执行时间)是在计算雅可比人.
这使我感到惊讶,因为大多数雅各比人都是从0和1向前传播,因此工作量应该是功能的2-4倍,而不是10-12倍.为了模拟大量的jacobian工作是什么样的,我用一个点积(而不是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根本不做任何优化.
我没有编译器的背景,但是有这个原因吗？我相当肯定,即使恒定折叠本身不会导致加速,在很大比例的科学计算中能够做更好的恒定传播/折叠会使更多的优化变得明显.
虽然我有兴趣解释为什么在编译器方面没有这样做,但我也对在实际方面可以做什么感兴趣,以便在面对这些模式时更快地使我自己的代码.
gga*_*ael 73
这是因为Eigen明确地将您的代码矢量化为3 vmulpd,2 vaddpd和其余4个组件寄存器中的1个水平缩减(这假定为AVX,只有SSE,您将获得6个mulpd和5个addpd).随着-ffast-mathGCC和铛允许删除最后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
Max*_*hof 38
我很失望地发现,如果没有启用快速数学运算,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:在Eigen上.据推测,他们对矩阵运算进行了矢量化,这对编译器来说要难以理解.使用以下选项正确优化了一个简单的循环:https://godbolt.org/z/OppEhY
Evg*_*Evg 12
强制编译器优化乘法0和1的一种方法是手动展开循环.为简单起见,我们使用
#include <array>
#include <cstddef>
constexpr std::size_t n = 12;
using Array = std::array<double, n>;
然后我们可以dot使用fold表达式实现一个简单的函数(如果它们不可用则递归):
<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-mathgcc 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没有并且需要一些帮助.