Copia fácil de caché de una matriz con reajuste por índice conocido, recopilar, dispersar
algorithm performance (7)
Supongamos que tenemos una matriz de datos y otra matriz con índices.
data = [1, 2, 3, 4, 5, 7]
index = [5, 1, 4, 0, 2, 3]
Queremos crear una nueva matriz a partir de elementos de data en la posición del index . El resultado debe ser
[4, 2, 5, 7, 3, 1]
El algoritmo ingenuo funciona para O (N) pero realiza un acceso aleatorio a la memoria.
¿Puede sugerir un algoritmo amigable de caché de CPU con la misma complejidad?
PD: En mi caso concreto, todos los elementos de la matriz de datos son enteros.
Las matrices PPS pueden contener millones de elementos.
PPPS Estoy de acuerdo con SSE / AVX o cualquier otra optimización específica de x64
data = [1, 2, 3, 4, 5, 7]
index = [5, 1, 4, 0, 2, 3]Queremos crear una nueva matriz a partir de elementos de datos en la posición del índice. El resultado debe ser
result -> [4, 2, 5, 7, 3, 1]
Un solo hilo, una pasada
Creo que, para unos pocos millones de elementos y en un solo hilo , el enfoque ingenuo podría ser el mejor aquí.
Tanto los data como el index se acceden (leen) secuencialmente, lo que ya es óptimo para la memoria caché de la CPU. Eso deja la escritura aleatoria, pero escribir en la memoria no es tan fácil como leer en la memoria caché de todos modos.
Esto solo necesitaría un pase secuencial a través de los datos y el índice. Y es probable que algunas (a veces muchas) de las escrituras también sean compatibles con el caché .
Usando múltiples bloques para el result - múltiples hilos
Podríamos asignar o usar bloques de tamaño compatibles con la memoria caché para el resultado (los bloques son regiones en la result array ), y hacer un bucle a través del index y los data varias veces (mientras permanecen en la memoria caché).
En cada bucle, solo escribimos los elementos para obtener un result que se ajuste al bloque de resultados actual . Esto sería ''compatible con la memoria caché'' para las escrituras también, pero necesita múltiples bucles (el número de bucles podría incluso ser bastante alto, es decir, el size of data / size of result-block ).
Lo anterior podría ser una opción cuando se usan varios subprocesos : los data y el index , de solo lectura, serían compartidos por todos los núcleos en algún nivel de la memoria caché (dependiendo de la arquitectura de la memoria caché). Los bloques de result en cada subproceso serían totalmente independientes (un núcleo nunca tiene que esperar el resultado de otro núcleo, o una escritura en la misma región). Por ejemplo: 10 millones de elementos: cada hilo podría estar trabajando en un bloque de resultados independiente de, por ejemplo, 500.000 elementos (el número debería ser una potencia de 2).
Combinar los valores como un par y ordenarlos primero: esto ya llevaría mucho más tiempo que la opción ingenua (y tampoco sería tan amigable con el caché).
Además, si solo hay unos pocos millones de elementos (enteros), no habrá mucha diferencia. Si estuviéramos hablando de miles de millones o datos que no caben en la memoria, otras estrategias podrían ser preferibles (como, por ejemplo, la asignación de memoria del conjunto de resultados si no cabe en la memoria).
Combine el índice y los datos en una sola matriz. Luego use algún algoritmo de clasificación que sea fácil de almacenar en caché para ordenar estos pares (por índice). Luego deshacerse de los índices. (Puede combinar la fusión / eliminación de índices con la primera / última pasada del algoritmo de clasificación para optimizar esto un poco).
Para una clasificación O (N) fácil de almacenar en caché, use la ordenación de radix con un radix suficientemente pequeño (como máximo la mitad de las líneas de caché en la caché de la CPU).
Aquí está la implementación en C del algoritmo similar a la ordenación de radix:
void reorder2(const unsigned size)
{
const unsigned min_bucket = size / kRadix;
const unsigned large_buckets = size % kRadix;
g_counters[0] = 0;
for (unsigned i = 1; i <= large_buckets; ++i)
g_counters[i] = g_counters[i - 1] + min_bucket + 1;
for (unsigned i = large_buckets + 1; i < kRadix; ++i)
g_counters[i] = g_counters[i - 1] + min_bucket;
for (unsigned i = 0; i < size; ++i)
{
const unsigned dst = g_counters[g_index[i] % kRadix]++;
g_sort[dst].index = g_index[i] / kRadix;
g_sort[dst].value = g_input[i];
__builtin_prefetch(&g_sort[dst + 1].value, 1);
}
g_counters[0] = 0;
for (unsigned i = 1; i < (size + kRadix - 1) / kRadix; ++i)
g_counters[i] = g_counters[i - 1] + kRadix;
for (unsigned i = 0; i < size; ++i)
{
const unsigned dst = g_counters[g_sort[i].index]++;
g_output[dst] = g_sort[i].value;
__builtin_prefetch(&g_output[dst + 1], 1);
}
}
Se diferencia del tipo de radix en dos aspectos: (1) no hace pases de conteo porque todos los contadores son conocidos de antemano; (2) evita el uso de valores de potencia de 2 para radix.
Este código C ++ se usó para la evaluación comparativa (si desea ejecutarlo en un sistema de 32 bits, disminuya ligeramente la constante kMaxSize ).
Aquí están los resultados de la prueba (en la CPU Haswell con 6Mb de caché):
Es fácil ver que las matrices pequeñas (por debajo de ~ 2 000 000 elementos) son compatibles con la memoria caché, incluso para algoritmos ingenuos. También es posible que note que el enfoque de clasificación comienza a ser hostil para el caché en el último punto del diagrama (con size/radix cerca de 0,75 líneas de caché en caché L3). Entre estos límites, el enfoque de clasificación es más eficiente que el algoritmo ingenuo.
En teoría (si comparamos solo el ancho de banda de memoria necesario para estos algoritmos con líneas de caché de 64 bytes y valores de 4 bytes), el algoritmo de clasificación debería ser 3 veces más rápido. En la práctica tenemos una diferencia mucho menor, alrededor del 20%. Esto podría mejorarse si usamos valores de 16 bits más pequeños para data matriz de data (en este caso, el algoritmo de clasificación es aproximadamente 1.5 veces más rápido).
Un problema más con el enfoque de clasificación es su comportamiento en el peor de los casos cuando el size/radix está cerca de alguna potencia de 2. Esto puede ser ignorado (porque no hay tantos tamaños "malos") o arreglado haciendo este algoritmo un poco más complicado.
Si aumentamos el número de pases a 3, los 3 pases usan principalmente caché L1, pero el ancho de banda de la memoria aumenta en un 60%. Utilicé este código para obtener resultados experimentales: TL; DR . Después de determinar (experimentalmente) el mejor valor de la raíz, obtuve resultados algo mejores para tamaños superiores a 4 000 000 (donde el algoritmo de 2 pases usa el caché L3 para una pasada) pero resultados algo peores para arreglos más pequeños (donde el algoritmo de 2 pasos usa L2 caché para ambas pasadas). Como puede esperarse, el rendimiento es mejor para datos de 16 bits.
Conclusión: la diferencia de rendimiento es mucho menor que la diferencia en la complejidad de los algoritmos, por lo que el enfoque ingenuo es casi siempre mejor; Si el rendimiento es muy importante y solo se utilizan valores de 2 o 4 bytes, es preferible el método de clasificación.
Me preocupa que esto no sea un patrón ganador.
Teníamos un código que funcionaba bien y lo optimizamos eliminando una copia.
El resultado fue que tuvo un bajo rendimiento (debido a problemas de caché). No puedo ver cómo se puede producir un algoritmo de una sola pasada que resuelva el problema. El uso de OpenMP puede permitir que las paradas se compartan entre varios subprocesos.
Noté que su índice cubre completamente el dominio pero está en orden aleatorio.
Si tuviera que ordenar el índice pero también aplicar las mismas operaciones a la matriz de índice a la matriz de datos, la matriz de datos se convertiría en el resultado que está buscando.
Hay una gran variedad de algoritmos para seleccionar, todos satisfacerían sus criterios de caché. Pero su complejidad varía. Consideraría o bien quicksort o mergesort.
Si estás interesado en esta respuesta puedo elaborar con pseudo código.
Si su problema tiene que ver con muchos más datos de los que muestra aquí, la forma más rápida, y probablemente la más compatible con la memoria caché, sería realizar una operación de clasificación de combinación grande y amplia.
Así que dividiría los datos de entrada en partes razonables, y tendría un hilo separado operando en cada parte. El resultado de esta operación sería dos matrices muy similares a la entrada (un índice de datos y un índice de destino), sin embargo, los índices se clasificarían. Luego tendría un subproceso final que realizaría una operación de combinación de los datos en la matriz de salida final.
Siempre que los segmentos se elijan bien, este debe ser un algoritmo muy fácil de almacenar en caché. Por sabiamente, quiero decir, de modo que los datos utilizados por los distintos subprocesos se asignen a diferentes líneas de caché (del procesador elegido) para evitar el apaleamiento del caché.
Si tiene muchos datos y ese es el cuello de botella, tendrá que usar un algoritmo basado en bloques donde lea y escriba lo más posible en los mismos bloques. Tomará hasta 2 pases sobre los datos para asegurarse de que la nueva matriz esté completamente llena y que el tamaño del bloque deba establecerse de manera apropiada. El pseudocódigo está abajo.
def populate(index,data,newArray,cache)
blockSize = 1000
for i = 0; i < size(index); i++
//We cached this value earlier
if i in cache
newArray[i] = cache[i]
remove(cache,i)
else
newIndex = index[i]
newValue = data[i]
//Check if this index is in our block
if i%blockSize != newIndex%blockSize
//This index is not in our current block, cache it
cache[newIndex] = newValue
else
//This value is in our current block
newArray[newIndex] = newValue
cache = {}
newArray = []
populate(index,data,newArray,cache)
populate(index,data,newArray,cache)
Análisis
La solución ingenua accede al índice y la matriz de datos en orden, pero se accede a la nueva matriz en orden aleatorio. Dado que se accede al nuevo array de manera aleatoria, esencialmente terminas con O (N ^ 2) donde N es el número de bloques en el array.
La solución basada en bloque no salta de bloque a bloque. Lee el índice, los datos y la nueva matriz, todos en secuencia para leer y escribir en los mismos bloques. Si un índice estará en otro bloque, se almacenará en caché y se recuperará cuando aparezca el bloque al que pertenece o si ya se ha pasado el bloque, se recuperará en la segunda pasada. Un segundo pase no dolerá en absoluto. Esto es O (N).
La única advertencia está en tratar con el caché. Aquí hay muchas oportunidades para ser creativo, pero en general, si muchas de las lecturas y escrituras terminan en bloques diferentes, la memoria caché crecerá y esto no es óptimo. Depende de la composición de sus datos, de la frecuencia con que esto ocurra y de la implementación de su caché.
Imaginemos que toda la información dentro del caché existe en un bloque y cabe en la memoria. Y digamos que el caché tiene elementos y. El enfoque ingenuo habría accedido al azar al menos y veces. El enfoque basado en bloques obtendrá a aquellos en la segunda pasada.
Supongo que el reordenamiento ocurre solo una vez de la misma manera. Si sucede varias veces, entonces crear una mejor estrategia de antemano (con un algoritmo de clasificación apropiado) mejorará el rendimiento
Escribí el siguiente programa para probar si una simple división del objetivo en N bloques ayuda, y mi descubrimiento fue:
a) incluso para los peores casos, no fue posible que el rendimiento de un solo hilo (el uso de escrituras segmentadas) no supere la estrategia ingenua, y suele ser peor por al menos un factor de 2
b) Sin embargo, el rendimiento se aproxima a la unidad para algunas subdivisiones (probablemente depende del procesador) y los tamaños de matriz, lo que indica que en realidad mejoraría el rendimiento de varios núcleos
La consecuencia de esto es: Sí, es más "fácil de almacenar en caché" que no subdividir, pero para un solo hilo (y solo un reordenamiento) esto no te ayudará un poco.
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
void main(char **ARGS,int ARGC) {
int N=1<<26;
double* source = malloc(N*sizeof(double));
double* target = malloc(N*sizeof(double));
int* idx = malloc(N*sizeof(double));
int i;
for(i=0;i<N;i++) {
source[i]=i;
target[i]=0;
idx[i] = rand() % N ;
};
struct timeval now,then;
gettimeofday(&now,NULL);
for(i=0;i<N;i++) {
target[idx[i]]=source[i];
};
gettimeofday(&then,NULL);
printf("%f/n",(0.0+then.tv_sec*1e6+then.tv_usec-now.tv_sec*1e6-now.tv_usec)/N);
gettimeofday(&now,NULL);
int j;
int targetblocks;
int M = 24;
int targetblocksize = 1<<M;
targetblocks = (N/targetblocksize);
for(i=0;i<N;i++) {
for(j=0;j<targetblocks;j++) {
int k = idx[i];
if ((k>>M) == j) {
target[k]=source[i];
};
};
};
gettimeofday(&then,NULL);
printf("%d,%f/n",targetblocks,(0.0+then.tv_sec*1e6+then.tv_usec-now.tv_sec*1e6-now.tv_usec)/N);
};