Reducciones en paralelo en tiempo logarítmico

algorithm parallel-processing

Question

Dadas n sumas parciales, es posible sumar todas las sumas parciales en pasos paralelos log2. Por ejemplo, suponga que hay ocho hilos con ocho sumas parciales: s0, s1, s2, s3, s4, s5, s6, s7 . Esto podría reducirse en log2(8) = 3 pasos secuenciales como este;

thread0     thread1    thread2    thread4
s0 += s1    s2 += s3   s4 += s5   s6 +=s7
s0 += s2    s4 += s6
s0 += s4

Me gustaría hacer esto con OpenMP pero no quiero usar la cláusula de reduction de OpenMP. Se me ocurrió una solución, pero creo que se puede encontrar una solución mejor tal vez utilizando la cláusula de task de OpenMP.

Esto es más general que la adición escalar. Permítanme elegir un caso más útil: una reducción de matriz (consulte here , here y here para obtener más información sobre las reducciones de matriz).

Digamos que quiero hacer una reducción de matriz en una matriz a . Aquí hay un código que llena las matrices privadas en paralelo para cada subproceso.

int bins = 20;
int a[bins];
int **at;  // array of pointers to arrays
for(int i = 0; i<bins; i++) a[i] = 0;
#pragma omp parallel
{
    #pragma omp single   
    at = (int**)malloc(sizeof *at * omp_get_num_threads());        
    at[omp_get_thread_num()] = (int*)malloc(sizeof **at * bins);
    int a_private[bins];
    //arbitrary function to fill the arrays for each thread
    for(int i = 0; i<bins; i++) at[omp_get_thread_num()][i] = i + omp_get_thread_num();
}

En este punto, tengo una matriz de punteros a matrices para cada subproceso. Ahora quiero agregar todas estas matrices juntas y escribir la suma final en a. Aquí está la solución que se me ocurrió.

#pragma omp parallel
{
    int n = omp_get_num_threads();
    for(int m=1; n>1; m*=2) {
        int c = n%2;
        n/=2;
        #pragma omp for
        for(int i = 0; i<n; i++) {
            int *p1 = at[2*i*m], *p2 = at[2*i*m+m];
            for(int j = 0; j<bins; j++) p1[j] += p2[j];
        }
        n+=c;
    }
    #pragma omp single
    memcpy(a, at[0], sizeof *a*bins);
    free(at[omp_get_thread_num()]);
    #pragma omp single
    free(at);
}

Déjame intentar y explicar lo que hace este código. Supongamos que hay ocho hilos. Definamos el operador += para significar sumar sobre la matriz. por ejemplo, s0 += s1 es

for(int i=0; i<bins; i++) s0[i] += s1[i]

entonces este código haría

n   thread0     thread1    thread2    thread4
4   s0 += s1    s2 += s3   s4 += s5   s6 +=s7
2   s0 += s2    s4 += s6
1   s0 += s4

Pero este código no es ideal como me gustaría.

Un problema es que hay algunas barreras implícitas que requieren que todos los hilos se sincronicen. Estas barreras no deberían ser necesarias. La primera barrera es entre llenar los arreglos y hacer la reducción. La segunda barrera está en la #pragma omp for declaración en la reducción. Pero no puedo usar la cláusula nowait con este método para eliminar la barrera.

Otro problema es que hay varios hilos que no necesitan ser utilizados. Por ejemplo con ocho hilos. El primer paso en la reducción solo necesita cuatro subprocesos, el segundo paso dos subprocesos y el último paso solo un subproceso. Sin embargo, este método implicaría los ocho hilos en la reducción. Sin embargo, los otros hilos no hacen mucho de todos modos y deberían ir directamente a la barrera y esperar, por lo que probablemente no sea un gran problema.

Mi instinto es que se puede encontrar un mejor método utilizando la cláusula de task omp. Desafortunadamente, tengo poca experiencia con la cláusula de task y todos mis esfuerzos hasta ahora con ella hacen una reducción mejor de lo que ahora he fallado.

¿Alguien puede sugerir una mejor solución para reducir el tiempo logarítmico utilizando, por ejemplo, la cláusula de task de OpenMP?

Encontré un método que resuelve el problema de la barrera. Esto se reduce de forma asincrónica. El único problema que queda es que todavía pone hilos que no participan en la reducción en un bucle ocupado. Este método utiliza algo así como una pila para empujar punteros a la pila (pero nunca los muestra) en secciones críticas (esta fue una de las claves ya que las secciones críticas no tienen barreras implícitas . La pila se opera en serie pero la reducción en paralelo .

Aquí hay un ejemplo de trabajo.

#include <stdio.h>
#include <omp.h>
#include <stdlib.h>
#include <string.h>

void foo6() {
    int nthreads = 13;
    omp_set_num_threads(nthreads);
    int bins= 21;
    int a[bins];
    int **at;
    int m = 0;
    int nsums = 0;
    for(int i = 0; i<bins; i++) a[i] = 0;
    #pragma omp parallel
    {
        int n = omp_get_num_threads();
        int ithread = omp_get_thread_num();
        #pragma omp single
        at = (int**)malloc(sizeof *at * n * 2);
        int* a_private = (int*)malloc(sizeof *a_private * bins);

        //arbitrary fill function
        for(int i = 0; i<bins; i++) a_private[i] = i + omp_get_thread_num();

        #pragma omp critical (stack_section)
        at[nsums++] = a_private;

        while(nsums<2*n-2) {
            int *p1, *p2;
            char pop = 0;
            #pragma omp critical (stack_section)
            if((nsums-m)>1) p1 = at[m], p2 = at[m+1], m +=2, pop = 1;
            if(pop) {
                for(int i = 0; i<bins; i++) p1[i] += p2[i];
                #pragma omp critical (stack_section)
                at[nsums++] = p1;
            }
        }

        #pragma omp barrier
        #pragma omp single
        memcpy(a, at[2*n-2], sizeof **at *bins);
        free(a_private);
        #pragma omp single
        free(at);
    }
    for(int i = 0; i<bins; i++) printf("%d ", a[i]); puts("");
    for(int i = 0; i<bins; i++) printf("%d ", (nthreads-1)*nthreads/2 +nthreads*i); puts("");
}

int main(void) {
    foo6();
}

Siento que se puede encontrar un mejor método usando tareas que no coloquen los hilos que no se usan en un bucle ocupado.

Answer 1

En realidad, es bastante simple de implementar de manera limpia con tareas utilizando un enfoque recursivo de divide y vencerás. Este es casi el código de un textbook .

void operation(int* p1, int* p2, size_t bins)
{
    for (int i = 0; i < bins; i++)
        p1[i] += p2[i];
}

void reduce(int** arrs, size_t bins, int begin, int end)
{
    assert(begin < end);
    if (end - begin == 1) {
        return;
    }
    int pivot = (begin + end) / 2;
    /* Moving the termination condition here will avoid very short tasks,
     * but make the code less nice. */
#pragma omp task
    reduce(arrs, bins, begin, pivot);
#pragma omp task
    reduce(arrs, bins, pivot, end);
#pragma omp taskwait
    /* now begin and pivot contain the partial sums. */
    operation(arrs[begin], arrs[pivot], bins);
}

/* call this within a parallel region */
#pragma omp single
reduce(at, bins, 0, n);

Por lo que puedo decir, no hay sincronizaciones innecesarias y no hay encuestas raras en secciones críticas. También funciona naturalmente con un tamaño de datos diferente al número de rangos. Lo encuentro muy limpio y fácil de entender. De hecho, creo que esto es mejor que sus dos soluciones.

Pero veamos cómo funciona en la práctica *. Para eso podemos usar Score-p y Vampir :

* _{bins=10000 por lo que la reducción en realidad lleva un poco de tiempo. Ejecutado en un sistema Haswell de 24 núcleos sin turbo. gcc 4.8.4, -O3. Agregué un poco de búfer alrededor de la ejecución real para ocultar la inicialización / postprocesamiento}

La imagen revela lo que está sucediendo en cualquier hilo dentro de la aplicación en un eje de tiempo horizontal. Las implementaciones de árbol de arriba a abajo:

1. omp for bucle
2. omp critical tipo omp critical de tareas.
3. omp task

Esto muestra muy bien cómo se ejecutan las implementaciones específicas. Ahora parece que el ciclo for es en realidad el más rápido, a pesar de las sincronizaciones innecesarias. Pero todavía hay una serie de fallas en este análisis de rendimiento. Por ejemplo, no fijé los hilos. En la práctica, NUMA (acceso de memoria no uniforme) es muy importante: ¿el núcleo tiene estos datos en su propia caché / memoria de su propio socket? Aquí es donde la solución de la tarea se vuelve no determinista. La variación muy significativa entre las repeticiones no se considera en la comparación simple.

Si la operación de reducción se vuelve variable en tiempo de ejecución, entonces la solución de la tarea será mejor que tu sincronizado para el bucle.

La solución critical tiene un aspecto interesante, los hilos pasivos no esperan continuamente, por lo que es más probable que consuman recursos de la CPU. Esto puede ser malo para el rendimiento, por ejemplo, en el caso del modo turbo.

Recuerde que la solución de la task tiene más potencial de optimización al evitar las tareas de generación que regresan de inmediato. El rendimiento de estas soluciones también depende en gran medida del tiempo de ejecución específico de OpenMP. El tiempo de ejecución de Intel parece ser mucho peor para las tareas.

Mi recomendación es:

• Implemente la solución más fácil de mantener con una complejidad algorítmica óptima.
• Mida qué partes del código realmente importan para el tiempo de ejecución
• Analice en base a mediciones reales cuál es el cuello de botella. En mi experiencia, se trata más de NUMA y la programación en lugar de una barrera innecesaria.
• Realice la microoptimización en función de sus mediciones reales

Solución lineal

Aquí está la línea de tiempo para el lineal proccess_data_v1 de esta pregunta .

Reducción de OpenMP 4

Entonces pensé en la reducción de OpenMP. La parte difícil parece ser obtener los datos de la matriz at dentro del bucle sin una copia. Inicializo la matriz de trabajo con NULL y simplemente muevo el puntero la primera vez:

void meta_op(int** pp1, int* p2, size_t bins)
{
    if (*pp1 == NULL) {
        *pp1 = p2;
        return;
    }
    operation(*pp1, p2, bins);
}

// ...

// declare before parallel region as global
int* awork = NULL;

#pragma omp declare reduction(merge : int* : meta_op(&omp_out, omp_in, 100000)) initializer (omp_priv=NULL)

#pragma omp for reduction(merge : awork)
        for (int t = 0; t < n; t++) {
            meta_op(&awork, at[t], bins);
        }

Sorprendentemente, esto no se ve muy bien:

_{la parte superior es icc 16.0.2 , la parte inferior es gcc 5.3.0 , ambas con -O3 .}

Ambos parecen implementar la reducción serializada. Traté de buscar en gcc / libgomp , pero no es inmediatamente evidente para mí lo que está sucediendo. Desde el código / desensamblaje intermedio, parecen estar envolviendo la fusión final en un GOMP_atomic_start / end , y eso parece ser un mutex global. De manera similar, icc envuelve la llamada a la operation en un kmpc_critical . Supongo que no hubo mucha optimización en las costosas operaciones de reducción personalizada. Una reducción tradicional se puede hacer con una operación atómica compatible con hardware.

Observe cómo cada operation es más rápida porque la entrada se almacena en caché localmente, pero debido a la serialización, en general es más lenta. Nuevamente, esta no es una comparación perfecta debido a las altas variaciones, y las capturas de pantalla anteriores se realizaron con diferentes versiones de gcc . Pero la tendencia es clara, y también tengo datos sobre los efectos de caché.