c++ - online - ide para lenguaje c

¿Por qué los compiladores de C++ no hacen mejor el plegado constante?

Question

Estoy investigando formas de acelerar una gran sección del código C ++, que tiene derivados automáticos para la computación de jacobianos. Esto implica hacer una cierta cantidad de trabajo en los residuos reales, pero la mayoría del trabajo (basado en el tiempo de ejecución perfilado) está en el cálculo de los jacobianos.

Esto me sorprendió, ya que la mayoría de los jacobianos se propagan hacia adelante desde 0 y 1, por lo que la cantidad de trabajo debe ser 2-4x la función, no 10-12x. Para modelar cómo es una gran cantidad de trabajo jacobiano, hice un ejemplo súper mínimo con solo un producto de puntos (en lugar de sin, cos, sqrt y más que estaría en una situación real) que el compilador debería poder para optimizar a un solo valor de retorno:

#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);
}

Que debería ser lo mismo que

double testReturnFirst(const Array12d& b)
{
    return b(0);
}

Me decepcionó descubrir que, sin la función de matemática rápida habilitada, ni GCC 8.2, Clang 6 ni MSVC 19 pudieron realizar optimizaciones en absoluto sobre el producto de puntos ingenuo con una matriz llena de 0s. Incluso con matemáticas rápidas ( https://godbolt.org/z/GvPXFy ) las optimizaciones son muy deficientes en GCC y Clang (aún implican multiplicaciones y adiciones), y MSVC no hace ninguna optimización en absoluto.

No tengo experiencia en compiladores, pero ¿hay alguna razón para esto? Estoy bastante seguro de que, en una gran proporción de los cálculos científicos, ser capaz de realizar una mejor propagación / plegado constante haría más aparentes las optimizaciones, incluso si el plegado constante en sí no tuviera como resultado una aceleración.

Si bien me interesan las explicaciones de por qué esto no se hace en el lado del compilador, también me interesa lo que puedo hacer en un lado práctico para hacer que mi propio código sea más rápido cuando me enfrento a este tipo de patrones.

Answer 1

Me decepcionó descubrir que, sin la función de matemática rápida habilitada, ni GCC 8.2, Clang 6 ni MSVC 19 pudieron realizar optimizaciones en absoluto sobre el producto de puntos ingenuo con una matriz llena de 0s.

Desafortunadamente, no tienen otra opción. Dado que los flotadores IEEE han firmado ceros, agregar 0.0 no es una operación de identidad:

-0.0 + 0.0 = 0.0 // Not -0.0!

Del mismo modo, multiplicar por cero no siempre da cero:

0.0 * Infinity = NaN // Not 0.0!

Por lo tanto, los compiladores simplemente no pueden realizar estos pliegues constantes en el producto de puntos mientras retienen el cumplimiento de flotación IEEE; por lo que saben, su entrada puede contener ceros y / o infinitos firmados.

Tendrá que usar -ffast-math para obtener estos pliegues, pero eso puede tener consecuencias no deseadas. Puede obtener un control más preciso con indicadores específicos (de http://gcc.gnu.org/wiki/FloatingPointMath ). De acuerdo con la explicación anterior, agregar las siguientes dos banderas debe permitir el plegado constante:
-ffinite-math-only , -fno-signed-zeros

De hecho, obtienes el mismo ensamblaje que con -ffast-math esta manera: https://godbolt.org/z/vGULLA . Solo abandonas los ceros firmados (probablemente irrelevantes), los NaN y los infinitos. Presumiblemente, si siguiera generándolos en su código, obtendría un comportamiento indefinido, así que sopese sus opciones.

En cuanto a por qué su ejemplo no está optimizado mejor incluso con -ffast-math : Eso está en Eigen. Es de suponer que tienen vectorización en sus operaciones de matriz, que son mucho más difíciles de ver para los compiladores. Un bucle simple está correctamente optimizado con estas opciones: https://godbolt.org/z/OppEhY

Answer 2

Esto se debe a que Eigen vectoriza explícitamente su código como 3 vmulpd, 2 vaddpd y 1 reducción horizontal dentro de los 4 registros de componentes restantes (esto supone AVX, con SSE solo obtendrá 6 mulpd y 5 addpd). Con -ffast-math GCC y clang pueden eliminar los últimos 2 vmulpd y vaddpd (y esto es lo que hacen) pero realmente no pueden reemplazar el vmulpd restante y la reducción horizontal que Eigen ha generado explícitamente.

Entonces, ¿qué EIGEN_DONT_VECTORIZE si deshabilita la vectorización explícita de Eigen definiendo EIGEN_DONT_VECTORIZE ? Luego obtiene lo que esperaba ( https://godbolt.org/z/UQsoeH ), pero otras piezas de código podrían volverse mucho más lentas.

Si desea deshabilitar localmente la vectorización explícita y no tiene miedo de meterse con el interno de Eigen, puede introducir una opción DontVectorize en Matrix y deshabilitar la vectorización especializando los traits<> para este tipo de 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>;

El ejemplo completo allí: https://godbolt.org/z/bOEyzv

Answer 3

Una forma de forzar a un compilador a optimizar las multiplicaciones por 0 y 1` es desenrollar manualmente el bucle. Por simplicidad usemos

#include <array>
#include <cstddef>
constexpr std::size_t n = 12;
using Array = std::array<double, n>;

Luego podemos implementar una función de dot simple usando expresiones de plegado (o recursión si no están disponibles):

<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>{});
}

Ahora echemos un vistazo a tu función

double test(const Array& b)
{
    const Array a{1};    // = {1, 0, ...}
    return dot(a, b);
}

Con -ffast-math gcc 8.2 produces :

test(std::array<double, 12ul> const&):
  movsd xmm0, QWORD PTR [rdi]
  ret

clang 6.0.0 va en la misma línea:

test(std::array<double, 12ul> const&): # @test(std::array<double, 12ul> const&)
  movsd xmm0, qword ptr [rdi] # xmm0 = mem[0],zero
  ret

Por ejemplo, para

double test(const Array& b)
{
    const Array a{1, 1};    // = {1, 1, 0...}
    return dot(a, b);
}

obtenemos

test(std::array<double, 12ul> const&):
  movsd xmm0, QWORD PTR [rdi]
  addsd xmm0, QWORD PTR [rdi+8]
  ret

Adición. Clang desenrolla un ciclo for (std::size_t i = 0; i < n; ++i) ... sin todos estos trucos de expresiones de plegado, gcc no necesita y necesita ayuda.