c++ - combinaciones entre arreglos y registros

Uso de vectores anidados frente a una envoltura de vector aplanado, comportamiento extraño

Question

El problema

Durante mucho tiempo tuve la impresión de que usar un std::vector<std::vector...> anidado para simular una matriz N-dimensional es, en general, malo, ya que la memoria no garantiza que sea contigua, y uno puede tiene errores de caché. Pensé que es mejor usar un vector plano y un mapa de múltiples dimensiones a 1D y viceversa. Así que decidí probarlo (código listado al final). Es bastante sencillo, programé la lectura / escritura en un vector 3D anidado frente a mi propia envoltura 3D de un vector 1D. Compilé el código con g++ y clang++ , con la optimización -O3 activada. Para cada ejecución cambié las dimensiones, por lo que puedo obtener una buena idea sobre el comportamiento. Para mi sorpresa, estos son los resultados que obtuve en mi máquina MacBook Pro (Retina, 13 pulgadas, finales de 2012), 2.5GHz i5, 8GB RAM, OS X 10.10.5:

g ++ 5.2

dimensions nested flat
X Y Z (ms) (ms)
100 100 100 -> 16 24
150 150 150 -> 58 98
200 200 200 -> 136 308
250 250 250 -> 264 746
300 300 300 -> 440 1537

clang ++ (LLVM 7.0.0)

dimensions nested flat
X Y Z (ms) (ms)
100 100 100 -> 16 18
150 150 150 -> 53 61
200 200 200 -> 135 137
250 250 250 -> 255 271
300 300 300 -> 423 477

Como puede ver, la envoltura "aplanar" nunca está superando a la versión anidada. Además, la implementación libstdc ++ de g ++ se desempeña bastante mal en comparación con la implementación libc ++, por ejemplo, para 300 x 300 x 300 la versión plana es casi 4 veces más lenta que la versión anidada. libc ++ parece tener el mismo rendimiento.

Mis preguntas:

1. ¿Por qué la versión plana no es más rápida? ¿No debería ser? ¿Me falta algo en el código de prueba?
2. Además, ¿por qué libstdc ++ de g ++ funciona tan mal al usar vectores planos? De nuevo, ¿no debería funcionar mejor?

El código que utilicé:

#include <chrono>
#include <cstddef>
#include <iostream>
#include <memory>
#include <random>
#include <vector>
// Thin wrapper around flatten vector
template<typename T>
class Array3D
{
std::size_t _X, _Y, _Z;
std::vector<T> _vec;
public:
Array3D(std::size_t X, std::size_t Y, std::size_t Z):
_X(X), _Y(Y), _Z(Z), _vec(_X * _Y * _Z) {}
T& operator()(std::size_t x, std::size_t y, std::size_t z)
{
return _vec[z * (_X * _Y) + y * _X + x];
}
const T& operator()(std::size_t x, std::size_t y, std::size_t z) const
{
return _vec[z * (_X * _Y) + y * _X + x];
}
};
int main(int argc, char** argv)
{
std::random_device rd{};
std::mt19937 rng{rd()};
std::uniform_real_distribution<double> urd(-1, 1);
const std::size_t X = std::stol(argv[1]);
const std::size_t Y = std::stol(argv[2]);
const std::size_t Z = std::stol(argv[3]);
// Standard library nested vector
std::vector<std::vector<std::vector<double>>>
vec3D(X, std::vector<std::vector<double>>(Y, std::vector<double>(Z)));
// 3D wrapper around a 1D flat vector
Array3D<double> vec1D(X, Y, Z);
// TIMING nested vectors
std::cout << "Timing nested vectors.../n";
auto start = std::chrono::steady_clock::now();
volatile double tmp1 = 0;
for (std::size_t x = 0 ; x < X; ++x)
{
for (std::size_t y = 0 ; y < Y; ++y)
{
for (std::size_t z = 0 ; z < Z; ++z)
{
vec3D[x][y][z] = urd(rng);
tmp1 += vec3D[x][y][z];
}
}
}
std::cout << "/tSum: " << tmp1 << std::endl; // we make sure the loops are not optimized out
auto end = std::chrono::steady_clock::now();
std::cout << "Took: ";
auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << ms << " milliseconds/n";
// TIMING flatten vector
std::cout << "Timing flatten vector.../n";
start = std::chrono::steady_clock::now();
volatile double tmp2 = 0;
for (std::size_t x = 0 ; x < X; ++x)
{
for (std::size_t y = 0 ; y < Y; ++y)
{
for (std::size_t z = 0 ; z < Z; ++z)
{
vec1D(x, y, z) = urd(rng);
tmp2 += vec1D(x, y, z);
}
}
}
std::cout << "/tSum: " << tmp2 << std::endl; // we make sure the loops are not optimized out
end = std::chrono::steady_clock::now();
std::cout << "Took: ";
ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << ms << " milliseconds/n";
}

EDITAR

Cambiando el Array3D<T>::operator() regresa a

return _vec[(x * _Y + y) * _Z + z];

De acuerdo con la sugerencia de @ 1201ProgramAlarm, de hecho se elimina el comportamiento "extraño" de g ++, en el sentido de que las versiones planas y anidadas se llevan aproximadamente al mismo tiempo. Sin embargo, sigue siendo intrigante. Pensé que el anidado sería mucho peor debido a problemas de caché. ¿Puedo tener suerte y tener toda la memoria asignada de forma contigua?

Answer 1

¿Puedo tener suerte y tener toda la memoria asignada de forma contigua?

Probablemente si. He modificado un poco su muestra, así que tenemos un punto de referencia que se concentra más en las diferencias entre los dos enfoques:

• el relleno de la matriz se realiza en una pasada separada, por lo que se excluye la velocidad del generador aleatorio
• eliminado volátil
• arreglado un pequeño error ( tmp1fue impreso en lugar de tmp2)
• se agregó una #if 1parte, que se puede usar para aleatorizar la vec3Dubicación en la memoria. Al final resultó que, tiene una gran diferencia en mi máquina.

Sin aleatorización (he usado parámetros: 300 300 300):

Timing nested vectors...
Sum: -131835
Took: 2122 milliseconds
Timing flatten vector...
Sum: -131835
Took: 2085 milliseconds

Así que hay una pequeña diferencia, la versión plana es un poco más rápido. (He hecho varias pruebas, y he puesto el tiempo mínimo aquí).

Con aleatorización:

Timing nested vectors...
Sum: -117685
Took: 3014 milliseconds
Timing flatten vector...
Sum: -117685
Took: 2085 milliseconds

Entonces, el efecto de caché se puede ver aquí: la versión anidada es ~ 50% más lenta . Esto se debe a que la CPU no puede predecir qué área de memoria se utilizará, por lo que su prefetcher de caché no es eficiente.

Aquí está el código modificado:

#include <chrono>
#include <cstddef>
#include <iostream>
#include <memory>
#include <random>
#include <vector>
template<typename T>
class Array3D
{
std::size_t _X, _Y, _Z;
std::vector<T> _vec;
public:
Array3D(std::size_t X, std::size_t Y, std::size_t Z):
_X(X), _Y(Y), _Z(Z), _vec(_X * _Y * _Z) {}
T& operator()(std::size_t x, std::size_t y, std::size_t z)
{
return _vec[(x * _Y + y) * _Z + z];
}
const T& operator()(std::size_t x, std::size_t y, std::size_t z) const
{
return _vec[(x * _Y + y) * _Z + z];
}
};
double nested(std::vector<std::vector<std::vector<double>>> &vec3D, std::size_t X, std::size_t Y, std::size_t Z) {
double tmp1 = 0;
for (int iter=0; iter<100; iter++)
for (std::size_t x = 0 ; x < X; ++x)
{
for (std::size_t y = 0 ; y < Y; ++y)
{
for (std::size_t z = 0 ; z < Z; ++z)
{
tmp1 += vec3D[x][y][z];
}
}
}
return tmp1;
}
double flatten(Array3D<double> &vec1D, std::size_t X, std::size_t Y, std::size_t Z) {
double tmp2 = 0;
for (int iter=0; iter<100; iter++)
for (std::size_t x = 0 ; x < X; ++x)
{
for (std::size_t y = 0 ; y < Y; ++y)
{
for (std::size_t z = 0 ; z < Z; ++z)
{
tmp2 += vec1D(x, y, z);
}
}
}
return tmp2;
}
int main(int argc, char** argv)
{
std::random_device rd{};
std::mt19937 rng{rd()};
std::uniform_real_distribution<double> urd(-1, 1);
const std::size_t X = std::stol(argv[1]);
const std::size_t Y = std::stol(argv[2]);
const std::size_t Z = std::stol(argv[3]);
std::vector<std::vector<std::vector<double>>>
vec3D(X, std::vector<std::vector<double>>(Y, std::vector<double>(Z)));
#if 1
for (std::size_t i = 0 ; i < X*Y; i++)
{
std::size_t xa = rand()%X;
std::size_t ya = rand()%Y;
std::size_t xb = rand()%X;
std::size_t yb = rand()%Y;
std::swap(vec3D[xa][ya], vec3D[xb][yb]);
}
#endif
// 3D wrapper around a 1D flat vector
Array3D<double> vec1D(X, Y, Z);
for (std::size_t x = 0 ; x < X; ++x)
{
for (std::size_t y = 0 ; y < Y; ++y)
{
for (std::size_t z = 0 ; z < Z; ++z)
{
vec3D[x][y][z] = vec1D(x, y, z) = urd(rng);
}
}
}
std::cout << "Timing nested vectors.../n";
auto start = std::chrono::steady_clock::now();
double tmp1 = nested(vec3D, X, Y, Z);
auto end = std::chrono::steady_clock::now();
std::cout << "/tSum: " << tmp1 << std::endl; // we make sure the loops are not optimized out
std::cout << "Took: ";
auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << ms << " milliseconds/n";
std::cout << "Timing flatten vector.../n";
start = std::chrono::steady_clock::now();
double tmp2 = flatten(vec1D, X, Y, Z);
end = std::chrono::steady_clock::now();
std::cout << "/tSum: " << tmp2 << std::endl; // we make sure the loops are not optimized out
std::cout << "Took: ";
ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
std::cout << ms << " milliseconds/n";
}

Answer 2

Al leer las otras respuestas, no estoy realmente satisfecho con la precisión y el nivel de detalle de las respuestas, así que intentaré una explicación yo mismo:

El tema del hombre aquí no es indirecto sino una cuestión de localidad espacial :

Básicamente, hay dos cosas que hacen que el almacenamiento en caché sea especialmente efectivo:

• Localidad temporal , lo que significa que probablemente se volverá a acceder a una palabra de memoria que se ha accedido recientemente en un futuro próximo. Esto podría ocurrir, por ejemplo, en nodos cerca de la raíz de un árbol de búsqueda binario al que se accede con frecuencia.

• Localidad espacial , lo que significa que si se ha accedido a una palabra de memoria, es probable que también se acceda pronto a las palabras de memoria anteriores o posteriores a esta palabra. Esto sucede en nuestro caso, para matrices anidadas y aplanadas.

Para evaluar el impacto que los efectos de indirección y caché podrían tener en este problema, supongamos que tenemos X = Y = Z = 1024

A juzgar por esta pregunta , una sola línea de caché (L1, L2 o L3) tiene 64 bytes de longitud, lo que significa 8 valores dobles. Asumamos que el caché L1 tiene 32 kB (4096 dobles), el caché L2 tiene 256 kB (32k dobles) y el caché L3 tiene 8 MB (1M dobles).

Esto significa que, suponiendo que el caché no contenga ningún otro dato (lo que es una suposición audaz, lo sé), en el caso aplanado solo cada cuarto valor de y conduce a una falta de caché L1 (la latencia del caché L2 es probablemente de alrededor de 10 20 ciclos), solo cada 32º valor de y conduce a una falla de caché L2 (la latencia de la caché L3 es un valor inferior a 100 ciclos) y solo en caso de una falla de caché L3 tenemos que acceder a la memoria principal. No quiero abrir todo el cálculo aquí, ya que tener en cuenta toda la jerarquía de caché lo hace un poco más difícil, pero digamos que casi todos los accesos a la memoria pueden almacenarse en caché en el caso aplanado.

En la formulación original de esta pregunta, el índice aplanado se calculó de manera diferente ( z * (_X * _Y) + y * _X + x ), un aumento del valor que cambia en el bucle más interno (z) siempre significa un salto de _X * _Y * 64 bit , lo que lleva a un diseño de memoria mucho más local que aumenta las fallas de caché en una gran cantidad.

En el caso anidado, la respuesta depende mucho del valor de Z:

• Si Z es bastante grande, la mayoría de los accesos a la memoria son contiguos, ya que se refieren a las entradas de un solo vector en el vector<vector<vector>>> anidado vector<vector<vector>>> , que están distribuidos de manera contigua. Solo cuando se incrementa el valor de y o x, necesitamos usar realmente el direccionamiento indirecto para recuperar el puntero de inicio del siguiente vector más interno.
• Si Z es bastante pequeño, necesitamos cambiar nuestro acceso a la memoria de "puntero base" con bastante frecuencia, lo que podría provocar errores de caché si las áreas de almacenamiento de los vectores más internos se colocan de forma bastante aleatoria en la memoria. Sin embargo, si son más o menos contiguos, podríamos observar diferencias mínimas o nulas en el rendimiento del caché.

Dado que había una pregunta sobre el resultado del ensamblaje, permítanme darles una breve descripción:

Si compara la salida de ensamblaje de la matriz anidada y aplanada, observará muchas similitudes: hay tres bucles anidados equivalentes y las variables de conteo x, y, z se almacenan en los registros. La única diferencia real, aparte del hecho de que la versión anidada usa dos contadores para cada índice externo para evitar la multiplicación por 24 en cada cálculo de la dirección, y la versión aplanada hace lo mismo para el bucle más interno y se multiplica por 8 : el bucle y donde, en lugar de simplemente incrementar y calcular el índice aplanado, necesitamos hacer tres cargas de memoria interdependientes para determinar el puntero base para nuestro bucle interno:

mov rax, QWORD PTR [rdi]
mov rax, QWORD PTR [rax+r11]
mov rax, QWORD PTR [rax+r10]

Pero como esto solo ocurre cada Zth vez y los punteros para el ''vector medio'' son probablemente almacenados en caché, la diferencia de tiempo es insignificante.

Answer 3

Es debido a cómo está ordenando sus índices en la clase 3D. Ya que tu bucle más interno está cambiando z, esa es la parte más grande de tu índice, por lo que obtienes muchas fallas de caché. Reorganizar su indexación a

_vec[(x * _Y + y) * _Z + z]

y deberías ver un mejor rendimiento.

Answer 4

(Esto realmente no responde a la pregunta. Creo que lo leí al revés inicialmente, asumiendo que el OP acaba de encontrar lo que esperaba, que los vectores anidados son más lentos que planos).

Debería esperar que la versión de vector anidado sea más lenta para cualquier otra cosa que no sea el acceso secuencial. Después de corregir el orden de indexación principal de fila / columna para su versión plana, debería ser más rápido para muchos usos, especialmente porque es más fácil para un compilador auto-vectorizar con SIMD sobre una gran variedad plana que sobre muchos cortos std::vector<>.

Una línea de caché es sólo 64B. Eso es 8 doubles. La ubicación en un nivel de página es importante debido a las entradas de TLB limitadas, y la captación previa requiere accesos secuenciales, pero de todos modos obtendrá eso (lo suficientemente cerca) con vectores anidados que se asignan todos a la vez con la mayoría de las implementaciones de malloc. (Esta es una marca de microbios trivial que no hace nada antes de asignar sus vectors. En un programa real que asigna y libera algo de memoria antes de hacer muchas asignaciones pequeñas, algunas de ellas pueden estar dispersas por más).

Además de la localidad, los niveles adicionales de direccionamiento indirecto son potencialmente problemáticos.

Una referencia / puntero a un std :: vector solo apunta al bloque de tamaño fijo que contiene el tamaño actual, el espacio asignado y el puntero al búfer. IDK si alguna implementación coloca el búfer justo después de los datos de control como parte del mismo bloque malloced, pero probablemente eso sea imposible porque sizeof(std::vector<int>)tiene que ser constante para que pueda tener un vector de vectores. Echa un vistazo a la asm en godbolt : una función que solo devuelve v[10]una carga con una matriz arg, pero dos cargas con un std :: vector arg.

En la implementación de vector anidado, la carga v[x][y][z]requiere 4 pasos (suponiendo que un puntero o una referencia a vya está en un registro).

• cargar v.buffer_pointeroo v.bpcomo sea que la implementacion lo llame. (Un puntero a una matriz de std::vector<std::vector<double>>)
• carga v.bp[x].bp(un puntero a una matriz de std::vector<double>)
• carga v.bp[x].bp[y].bp(un puntero a una matriz de double)
• cargar v.bp[x].bp[y].bp[z](lo doubleque queramos)

Una matriz 3D adecuada, simulada con una sola std::vector, simplemente hace:

• carga v.bp(un puntero a una matriz de double)
• cargar v.bp[(x*w + y)*h + z](lo doubleque queramos)

Los accesos múltiples al mismo arreglo 3D simulado con diferentes xey requieren el cálculo de un nuevo índice, pero v.bppermanecerán en un registro. Así que en lugar de 3 fallos de caché, solo obtenemos uno.

Al atravesar la matriz 3D en orden se oculta la penalización de la implementación del vector anidado, porque hay un bucle sobre todos los valores en el vector más interno que oculta la sobrecarga de cambiar x e y. La captación previa de los punteros contiguos en los vectores externos ayuda aquí, y Zes lo suficientemente pequeño en sus pruebas de que el bucle sobre un vector más interno no expulsará el puntero para el siguiente yvalor.

Lo que todo programador debe saber sobre la memoria se está quedando un poco desactualizado, pero cubre los detalles del almacenamiento en caché y la localidad. La captación previa de software no es tan importante como lo fue en P4, así que no le prestes demasiada atención a esa parte de la guía.

Answer 5

¿Por qué los vectores anidados tienen aproximadamente la misma velocidad que los planos en su microbenchmark, después de corregir el orden de indexación ? Usted esperaría que la matriz plana fuera más rápida (vea la respuesta de Tobias sobre posibles problemas de localidad , y mi otra respuesta para saber por qué los vectores anidados chupan en general, pero no tan mal para el acceso secuencial). Pero su prueba específica está haciendo tantas cosas que permiten que la ejecución fuera de orden oculte la sobrecarga de usar vectores anidados y / o que solo ralentice las cosas tanto que la sobrecarga adicional se pierda en el ruido de medición.

Puse su código fuente de corrección de errores de rendimiento en Godbolt para que podamos ver el asm del bucle interno compilado por g ++ 5.2, con -O3 . (La bifurcación de Apple puede ser similar a clang3.7, pero solo veré la versión de gcc). Hay un montón de código de las funciones de C ++, pero puede hacer clic derecho en una línea de origen para desplazarse por las ventanas de ASM a El código para esa línea. Además, pase el mouse sobre una línea de origen para poner en negrita el asm que implementa esa línea o viceversa.

Los dos bucles internos de gcc para la versión anidada son los siguientes (con algunos comentarios agregados a mano):

## outer-most loop not shown
.L213: ## middle loop (over `y`)
test rbp, rbp # Z
je .L127 # inner loop runs zero times if Z==0
mov rax, QWORD PTR [rsp+80] # MEM[(struct vector * *)&vec3D], MEM[(struct vector * *)&vec3D]
xor r15d, r15d # z = 0
mov rax, QWORD PTR [rax+r12] # MEM[(struct vector * *)_195], MEM[(struct vector * *)_195]
mov rdx, QWORD PTR [rax+rbx] # D.103857, MEM[(double * *)_38]
## Top of inner-most loop.
.L128:
lea rdi, [rsp+5328] # tmp511, ## function arg: pointer to the RNG object, which is a local on the stack.
lea r14, [rdx+r15*8] # D.103851, ## r14 = &(vec3D[x][y][z])
call double std::generate_canonical<double, 53ul, std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >(std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul>&) #
addsd xmm0, xmm0 # D.103853, D.103853 ## return val *= 2.0: [0.0, 2.0]
mov rdx, QWORD PTR [rsp+80] # MEM[(struct vector * *)&vec3D], MEM[(struct vector * *)&vec3D] ## redo the pointer-chasing from vec3D.data()
mov rdx, QWORD PTR [rdx+r12] # MEM[(struct vector * *)_150], MEM[(struct vector * *)_150]
subsd xmm0, QWORD PTR .LC6[rip] # D.103859, ## and subtract 1.0: [-1.0, 1.0]
mov rdx, QWORD PTR [rdx+rbx] # D.103857, MEM[(double * *)_27]
movsd QWORD PTR [r14], xmm0 # *_155, D.103859 # store into vec3D[x][y][z]
movsd xmm0, QWORD PTR [rsp+64] # D.103853, tmp1 # reload volatile tmp1
addsd xmm0, QWORD PTR [rdx+r15*8] # D.103853, *_62 # add the value just stored into the array (r14 = rdx+r15*8 because nothing else modifies the pointers in the outer vectors)
add r15, 1 # z,
cmp rbp, r15 # Z, z
movsd QWORD PTR [rsp+64], xmm0 # tmp1, D.103853 # spill tmp1
jne .L128 #,
#End of inner-most loop
.L127: ## middle-loop
add r13, 1 # y,
add rbx, 24 # sizeof(std::vector<> == 24) == the size of 3 pointers.
cmp QWORD PTR [rsp+8], r13 # %sfp, y
jne .L213 #,
## outer loop not shown.

Y para el bucle plano:

 ## outer not shown.
.L214:
test rbp, rbp # Z
je .L135 #,
mov rax, QWORD PTR [rsp+280] # D.103849, vec1D._Y
mov rdi, QWORD PTR [rsp+288] # D.103849, vec1D._Z
xor r15d, r15d # z
mov rsi, QWORD PTR [rsp+296] # D.103857, MEM[(double * *)&vec1D + 24B]
.L136: ## inner-most loop
imul rax, r12 # D.103849, x
lea rax, [rax+rbx] # D.103849,
imul rax, rdi # D.103849, D.103849
lea rdi, [rsp+5328] # tmp520,
add rax, r15 # D.103849, z
lea r14, [rsi+rax*8] # D.103851, # &vec1D(x,y,z)
call double std::generate_canonical<double, 53ul, std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul> >(std::mersenne_twister_engine<unsigned long, 32ul, 624ul, 397ul, 31ul, 2567483615ul, 11ul, 4294967295ul, 7ul, 2636928640ul, 15ul, 4022730752ul, 18ul, 1812433253ul>&) #
mov rax, QWORD PTR [rsp+280] # D.103849, vec1D._Y
addsd xmm0, xmm0 # D.103853, D.103853
mov rdi, QWORD PTR [rsp+288] # D.103849, vec1D._Z
mov rsi, QWORD PTR [rsp+296] # D.103857, MEM[(double * *)&vec1D + 24B]
mov rdx, rax # D.103849, D.103849
imul rdx, r12 # D.103849, x # redo address calculation a 2nd time per iteration
subsd xmm0, QWORD PTR .LC6[rip] # D.103859,
add rdx, rbx # D.103849, y
imul rdx, rdi # D.103849, D.103849
movsd QWORD PTR [r14], xmm0 # MEM[(double &)_181], D.103859 # store into the address calculated earlier
movsd xmm0, QWORD PTR [rsp+72] # D.103853, tmp2
add rdx, r15 # tmp374, z
add r15, 1 # z,
addsd xmm0, QWORD PTR [rsi+rdx*8] # D.103853, MEM[(double &)_170] # tmp2 += vec1D(x,y,z). rsi+rdx*8 == r14, so this is a reload of the store this iteration.
cmp rbp, r15 # Z, z
movsd QWORD PTR [rsp+72], xmm0 # tmp2, D.103853
jne .L136 #,
.L135: ## middle loop: increment y
add rbx, 1 # y,
cmp r13, rbx # Y, y
jne .L214 #,
## outer loop not shown.

Su MacBook Pro (finales de 2012) tiene una CPU Intel IvyBridge , por lo que estoy usando los números para esa microarquitectura de las tablas de instrucciones de Agner Fog y la guía de microarca . Las cosas deberían ser casi iguales en otras CPU Intel / AMD.

El único IvB i5 móvil de 2.5GHz es el i5-3210M, por lo que su CPU tiene 3MiB de caché L3. Esto significa que incluso su caso de prueba más pequeño (100 ^ 3 * 8B por double ~ = 7.63MiB) es más grande que su caché de último nivel, por lo que ninguno de sus casos de prueba cabe en el caché. Probablemente sea algo bueno, ya que asigna e inicializa por defecto tanto anidado como plano antes de probar cualquiera de ellos. Sin embargo, realiza la prueba en el mismo orden que asignó, por lo que si la matriz anidada sigue siendo caché después de poner a cero la matriz plana, la matriz plana aún puede estar activa en la caché L3 después del bucle de tiempo sobre la matriz anidada.

Si hubiera utilizado un ciclo de repetición para recorrer la misma matriz varias veces, podría haber tenido el tiempo suficiente para medir tamaños de matriz más pequeños.

Estás haciendo varias cosas que son súper raras y lo hacen tan lento que la ejecución fuera de orden puede ocultar la latencia adicional de cambiar y , incluso si tus vectores z internos no son perfectamente contiguos.

1. Se ejecuta un PRNG lento dentro del bucle temporizado. std::uniform_real_distribution<double> urd(-1, 1); es una sobrecarga adicional sobre std::mt19937 rng{rd()}; , que ya es lento en comparación con la latencia de FP-add (3 ciclos), o en comparación con el rendimiento de carga de caché L1D de 2 por ciclo. Todo este tiempo adicional de ejecución del PRNG brinda a la ejecución fuera de orden la oportunidad de ejecutar las instrucciones de indexación de matrices para que la dirección final esté lista en el momento en que están los datos. A menos que tenga muchos errores de caché, en su mayoría solo está midiendo la velocidad de PRNG , ya que produce resultados mucho más lentos que 1 por ciclo de reloj.

   g ++ 5.2 no incluye completamente el código urd(rng) , y la convención de llamadas del sistema V x86-64 no tiene registros de XMM con preservación de llamadas. Por tmp1 tanto, tmp1 / tmp2 se debe derramar / recargar para cada elemento, incluso si no fueran volatile .

   También pierde su lugar en el vector Z, y tiene que rehacer los 2 niveles externos de direccionamiento indirecto antes de acceder al siguiente elemento z . Esto se debe a que no conoce la parte interna de la función a la que está llamando y asume que podría tener un puntero a la memoria del vector<> externo vector<> . (La versión plana hace dos multiplicaciones para la indexación en el bucle interno, en lugar de un simple puntero-agregar).

   clang (con libc ++) hace que el PRNG esté completamente en línea, por lo que pasar a la siguiente z es simplemente add reg, 8 para incrementar un puntero en las versiones plana y anidada. Puede obtener el mismo comportamiento de gcc si obtiene un iterador fuera del bucle interno, o si obtiene una referencia al vector interno, en lugar de rehacer el operator[] y esperando que el compilador lo levante por usted.

   El rendimiento / latencia de Intel / AMD FP add / sub / mul no depende de los datos, excepto para denormales. ( x87 también disminuye la velocidad para NaN y tal vez infinito , pero SSE no. El código de 64 bits usa SSE incluso para escalar / float / escalar). Por lo tanto, podría haber inicializado su matriz con ceros, o con un PRNG fuera del ciclo de tiempo. . (O déjelos en cero, ya que el vector<double> constructor lo hace por usted, y en realidad se necesita un código adicional para no hacerlo en los casos en los que va a escribir otra cosa). El rendimiento de la división y el sqrt depende de los datos en algunas CPU, y mucho más lento que agregar / sub / mul.

2. Escribe cada elemento justo antes de leerlo, dentro del bucle interno . En la fuente, esto parece una tienda / recarga. Desafortunadamente, eso es lo que hace gcc en realidad, pero clang con libc ++ (que integra el PRNG) transforma el cuerpo del bucle:

    // original
   vec3D[x][y][z] = urd(rng);
   tmp1 += vec3D[x][y][z];
   // what clang''s asm really does
   double xmm7 = urd(rng);
   vec3D[x][y][z] = xmm7;
   tmp1 += xmm7;

   En el asm de Clang:

    # do { ...
   addsd xmm7, xmm4 # last instruction of the PRNG
   movsd qword ptr [r8], xmm7 # store it into the Z vector
   addsd xmm7, qword ptr [rsp + 88]
   add r8, 8 # pointer-increment to walk along the Z vector
   dec r13 # i--
   movsd qword ptr [rsp + 88], xmm7
   jne .LBB0_74 # }while(i != 0);

   Se permite hacer esto porque vec3D no es volatile o atomic<> , por lo que sería un comportamiento indefinido para cualquier otro hilo que esté escribiendo esta memoria al mismo tiempo. Esto significa que puede optimizar el almacenamiento / recarga de objetos en memoria en solo un almacén (y simplemente usar el valor que almacenó, sin recargar). O bien, optimice la tienda por completo si puede demostrar que es una tienda muerta (una tienda que nada puede leer, por ejemplo, a una variable static no utilizada).

   En la versión de gcc , realiza la indexación para la tienda antes de la llamada PRNG y la indexación para la recarga después. Así que creo que gcc no está seguro de que la llamada a la función no modifique un puntero, porque los punteros a los vectores externos han escapado a la función. (Y el PRNG no está en línea).

   Sin embargo, incluso un almacenamiento / recarga real en el asm es aún menos sensible a los fallos de caché que a una simple carga.

   Almacenar-> el reenvío de carga todavía funciona incluso si la tienda pierde en el caché. Por lo tanto, una falla de caché en un vector Z no retrasa directamente la ruta crítica. Solo lo ralentiza si la ejecución fuera de orden no puede ocultar la latencia de la falta de caché. (Una tienda puede retirarse tan pronto como los datos se escriben en el almacén-búfer (y todas las instrucciones anteriores se han retirado). No estoy seguro de si una carga puede retirarse antes de que la línea de caché llegue a L1D, si se ha puede ser capaz de hacerlo, ya que x86 permite la reordenación de StoreLoad (se permite que las tiendas se vuelvan visibles globalmente después de las cargas). En ese caso, una tienda / recarga solo agrega 6 ciclos de latencia para el resultado de PRNG ( fuera de la ruta crítica de un estado de PRNG al siguiente estado de PRNG). Y es solo el rendimiento de un cuello de botella si la memoria caché se pierde tanto que el búfer de la tienda se llena y evita que se ejecuten nuevos uops del almacén, lo que a su vez evita que los uops nuevos emitirse en el núcleo fuera de orden cuando la Estación de Reservaciones o ROB se llena con uops no ejecutados o no retirados (respectivamente).

   Con la indexación invertida (versión original del código plano), probablemente el principal cuello de botella fueron las tiendas dispersas. IDK porque clang lo hizo mucho mejor que gcc allí. Tal vez clang logró invertir un bucle y atravesar la memoria en orden secuencial después de todo. (Dado que incorporaba completamente el PRNG, no había llamadas de función que requirieran que el estado de la memoria coincidiera con el orden del programa).

   Atravesar cada vector Z en orden significa que las fallas de caché están relativamente alejadas (incluso si cada vector Z no es contiguo al anterior), lo que da mucho tiempo para que se ejecuten las tiendas. O incluso si una carga reenviada desde la tienda no puede retirarse hasta que la caché L1D sea la propietaria de la línea de caché (en el estado Modificado del protocolo MESI), la ejecución especulativa tiene los datos correctos y no tuvo que esperar la latencia de la caché-señorita. La ventana de instrucciones fuera de orden es probablemente lo suficientemente grande como para evitar que la ruta crítica se detenga antes de que la carga pueda retirarse. (Las cargas de fallas de caché normalmente son muy malas, porque las instrucciones dependientes no pueden ejecutarse sin datos para que puedan operar. Por lo tanto, crean burbujas en la tubería. Con una falta de caché completa de DRAM que tiene una latencia de más de 300 ciclos, y la ventana fuera de orden es 168 uops en IvB, no puede ocultar toda la latencia para que el código se ejecute a 1 uop (aproximadamente 1 instrucción) por reloj.) Para almacenes puros, fuera de orden La ventana de orden se extiende más allá del tamaño del ROB, ya que no necesitan comprometerse con L1D para retirarse. De hecho, no pueden comprometerse hasta después de su jubilación, porque ese es el punto en el que se sabe que no son especulativos. (Por lo tanto, hacerlos visibles globalmente antes de eso evitaría un retroceso en la detección de una excepción o especulación errónea).

   No tengo libc++ instalado en mi escritorio, por lo que no puedo comparar esa versión con g ++. Con g ++ 5.4, encuentro Nested: 225 milisegundos y Flat: 239 milisegundos. Sospecho que los multiplicadores adicionales de indexación de matrices son un problema, y ​​compito con las instrucciones ALU que usa el PRNG. En contraste, la versión anidada que rehace un montón de persecución de punteros que se encuentran en la memoria caché L1D puede suceder en paralelo. Mi escritorio es un Skylake i7-6700k a 4.4GHz. SKL tiene un tamaño de ROB (Reordenar Buffer) de 224 uops, y RS de 97 uops, por lo que la ventana fuera de orden es muy grande . También tiene una latencia FP-add de 4 ciclos (a diferencia de los uarches anteriores donde era 3).

3. volatile double tmp1 = 0; Su acumulador es volatile , lo que obliga al compilador a almacenarlo / recargarlo en cada iteración del bucle interno. La latencia total de la cadena de dependencia transmitida por el bucle en el bucle interno es de 9 ciclos: 3 para addsd y 6 para el reenvío desde movsd al almacen hasta la recarga de movsd . (Clang pliega la recarga en un operando de memoria con addsd xmm7, qword ptr [rsp + 88] , pero la misma diferencia. ( [rsp+88] está en la pila, donde se almacenan las variables con almacenamiento automático, si es necesario derramarlas de registros.)

   Como se indicó anteriormente, la llamada a la función no en línea para gcc también forzará un derrame / recarga en la convención de llamadas del sistema V x86-64 (utilizada por todo menos en Windows). Pero un compilador inteligente podría haber hecho 4 llamadas PRNG, por ejemplo, y luego 4 almacenes de arreglos. (Si hubieras usado un iterador para asegurarte de que gcc sabía que los vectores que contenían otros vectores no estaban cambiando).

   El uso de -ffast-math hubiera permitido que el compilador se vectorizara automáticamente (si no fuera por PRNG y volatile ). Esto le permitiría ejecutar los arreglos lo suficientemente rápido como para que la falta de ubicación entre diferentes vectores Z pueda ser un problema real. También permitiría a los compiladores desenrollarse con múltiples acumuladores, para ocultar la latencia de FP-add. por ejemplo, podrían (y clang) hacer asm equivalente a:

   float t0=0, t1=0, t2=0, t3=0;
   for () {
   t0 += a[i + 0];
   t1 += a[i + 1];
   t2 += a[i + 2];
   t3 += a[i + 3];
   }
   t0 = (t0 + t1) + (t2 + t3);

   Eso tiene 4 cadenas de dependencia separadas, por lo que puede mantener 4 FP agregados en vuelo. Como IvB tiene una latencia de 3 ciclos, una por rendimiento de reloj para addsd , solo necesitamos mantener 4 en vuelo para saturar su rendimiento. (Skylake tiene 4c de latencia, 2 por rendimiento de reloj, lo mismo que mul o FMA, por lo que necesita 8 acumuladores para evitar cuellos de botella de latencia. En realidad, aún más es mejor . Como lo hizo la persona que hizo la pregunta, Haswell lo hizo mejor con más acumuladores cuando se acerque a maximizar el rendimiento de carga)

   Algo así sería una prueba mucho mejor de cuán eficiente es hacer un bucle en un Array3D. Si desea evitar que el bucle se optimice por completo, solo use el resultado . Pruebe su marca microbiológica para asegurarse de que al aumentar el tamaño del problema el tiempo cambie; si no, entonces algo se optimizó, o no estás probando lo que crees que estás probando. ¡No hagas que el bucle interno sea volatile !

Escribir microbenchmarks no es fácil. Debe comprender lo suficiente como para escribir uno que pruebe lo que cree que está probando. : P Este es un buen ejemplo de lo fácil que es equivocarse.

¿Puedo tener suerte y tener toda la memoria asignada de forma contigua?

Sí, eso probablemente suceda con muchas asignaciones pequeñas hechas en orden, cuando no ha asignado y liberado nada antes de hacer esto. Si fueran lo suficientemente grandes (generalmente una página de 4 kB o más), glibc malloc cambiaría a usar mmap(MAP_ANONYMOUS) y luego el núcleo elegiría direcciones virtuales aleatorias ( ASLR ). Entonces, con una Z más grande, se puede esperar que la localidad empeore. Pero, por otro lado, los vectores Z más grandes significan que pasas más tiempo en bucle en un vector contiguo, por lo que una falla de caché al cambiar y ( x ) se vuelve relativamente menos importante.

Aparentemente, hacer un bucle secuencial sobre sus datos con su información no lo expone, porque el puntero adicional accede a la memoria caché, por lo que la persecución del puntero tiene una latencia lo suficientemente baja para que la ejecución del OOO lo oculte con su ciclo lento.

Prefetch tiene un tiempo muy fácil mantenerse al día aquí.

Diferentes compiladores / bibliotecas pueden hacer una gran diferencia con esta prueba extraña. En mi sistema (Arch Linux, i7-6700k Skylake con 4.4GHz turbo máx.), El mejor de 4 ejecuciones a 300 300 300 para g ++ 5.4 -O3 fue:

Timing nested vectors...
Sum: 579.78
Took: 225 milliseconds
Timing flatten vector...
Sum: 579.78
Took: 239 milliseconds
Performance counter stats for ''./array3D-gcc54 300 300 300'':
532.066374 task-clock (msec) # 1.000 CPUs utilized
2 context-switches # 0.004 K/sec
0 cpu-migrations # 0.000 K/sec
54,523 page-faults # 0.102 M/sec
2,330,334,633 cycles # 4.380 GHz
7,162,855,480 instructions # 3.07 insn per cycle
632,509,527 branches # 1188.779 M/sec
756,486 branch-misses # 0.12% of all branches
0.532233632 seconds time elapsed

vs. g ++ 7.1 -O3 (que aparentemente decidió dividirse en algo que g ++ 5.4 no hizo)

Timing nested vectors...
Sum: 932.159
Took: 363 milliseconds
Timing flatten vector...
Sum: 932.159
Took: 378 milliseconds
Performance counter stats for ''./array3D-gcc71 300 300 300'':
810.911200 task-clock (msec) # 1.000 CPUs utilized
0 context-switches # 0.000 K/sec
0 cpu-migrations # 0.000 K/sec
54,523 page-faults # 0.067 M/sec
3,546,467,563 cycles # 4.373 GHz
7,107,511,057 instructions # 2.00 insn per cycle
794,124,850 branches # 979.299 M/sec
55,074,134 branch-misses # 6.94% of all branches
0.811067686 seconds time elapsed

vs. clang4.0 -O3 (con libstdc ++ de gcc, no libc ++)

perf stat ./array3D-clang40-libstdc++ 300 300 300
Timing nested vectors...
Sum: -349.786
Took: 1657 milliseconds
Timing flatten vector...
Sum: -349.786
Took: 1631 milliseconds
Performance counter stats for ''./array3D-clang40-libstdc++ 300 300 300'':
3358.297093 task-clock (msec) # 1.000 CPUs utilized
9 context-switches # 0.003 K/sec
0 cpu-migrations # 0.000 K/sec
54,521 page-faults # 0.016 M/sec
14,679,919,916 cycles # 4.371 GHz
12,917,363,173 instructions # 0.88 insn per cycle
1,658,618,144 branches # 493.887 M/sec
916,195 branch-misses # 0.06% of all branches
3.358518335 seconds time elapsed

No me interesé en qué fue lo que Clang hizo mal, o probé con -ffast-math y / or -march=native . (Sin embargo, esos no harán mucho a menos que elimines la volatile ).

perf stat -d no muestra más errores de caché (L1 o último nivel) para clang que gcc. Pero sí muestra que Clang está haciendo más del doble de cargas de L1D.

Lo intenté con una matriz no cuadrada. Es casi exactamente al mismo tiempo mantener el número total de elementos igual, pero cambiando la dimensión final a 5 o 6.

Incluso un cambio menor a la C ayuda, y se "aplana" más rápido que el anidado con gcc (desde 240 ms hasta 220 ms para 300 ^ 3, pero sin apenas hacer una diferencia para anidado):

 // vec1D(x, y, z) = urd(rng);
double res = urd(rng);
vec1D(x, y, z) = res; // indexing calculation only done once, after the function call
tmp2 += vec1D(x, y, z);
// using iterators would still avoid redoing it at all.