c - dollars - openmp ubuntu

¿Por qué este código no se escala linealmente?

Question

Escribí este código de solucionador SOR. No se preocupe demasiado por lo que hace este algoritmo, no es la preocupación aquí. Pero solo por el bien de la integridad: puede resolver un sistema lineal de ecuaciones, dependiendo de qué tan bien condicionado esté el sistema.

Lo ejecuto con una matriz de separación de 2097152 filas mal condicionada (que nunca converge), con un máximo de 7 columnas distintas de cero por fila.

Traducción: el bucle externo do-while while realizará 10000 iteraciones (el valor que paso como max_iters ), el medio for 2097152 iteraciones, dividido en trozos de work_line , dividido entre los subprocesos OpenMP. El bucle interno más interno tendrá 7 iteraciones, excepto en muy pocos casos (menos del 1%) donde puede ser menor.

Existe una dependencia de datos entre los hilos en los valores de sol array. Cada iteración del medio for actualizar un elemento, pero lee hasta otros 6 elementos de la matriz. Dado que SOR no es un algoritmo exacto, al leer, puede tener cualquiera de los valores anteriores o actuales en esa posición (si está familiarizado con los solucionadores, este es un Gauss-Siedel que tolera el comportamiento de Jacobi en algunos lugares por el bien de paralelismo).

typedef struct{
    size_t size;

    unsigned int *col_buffer;
    unsigned int *row_jumper;
    real *elements;
} Mat;

int work_line;

// Assumes there are no null elements on main diagonal
unsigned int solve(const Mat* matrix, const real *rhs, real *sol, real sor_omega, unsigned int max_iters, real tolerance)
{
    real *coefs = matrix->elements;
    unsigned int *cols = matrix->col_buffer;
    unsigned int *rows = matrix->row_jumper;
    int size = matrix->size;
    real compl_omega = 1.0 - sor_omega;
    unsigned int count = 0;
    bool done;

    do {
        done = true;
        #pragma omp parallel shared(done)
        {
            bool tdone = true;

            #pragma omp for nowait schedule(dynamic, work_line)
            for(int i = 0; i < size; ++i) {
                real new_val = rhs[i];
                real diagonal;
                real residual;
                unsigned int end = rows[i+1];

                for(int j = rows[i]; j < end; ++j) {
                    unsigned int col = cols[j];
                    if(col != i) {
                        real tmp;
                        #pragma omp atomic read
                        tmp = sol[col];

                        new_val -= coefs[j] * tmp;
                    } else {
                        diagonal = coefs[j];
                    }
                }

                residual = fabs(new_val - diagonal * sol[i]);
                if(residual > tolerance) {
                    tdone = false;
                }

                new_val = sor_omega * new_val / diagonal + compl_omega * sol[i];
                #pragma omp atomic write
                sol[i] = new_val;
            }

            #pragma omp atomic update
            done &= tdone;
        }
    } while(++count < max_iters && !done);

    return count;
}

Como puede ver, no hay bloqueo dentro de la región paralela, por lo que, como siempre nos enseñaron, es el tipo de problema paralelo del 100%. Eso no es lo que veo en la práctica.

Todas mis pruebas se ejecutaron en una CPU Intel (R) Xeon (R) E5-2670 v2 a 2.50GHz, 2 procesadores, 10 núcleos cada uno, hiper-hilo habilitado, sumando hasta 40 núcleos lógicos.

En mis primeras series de work_line , work_line se fijó en 2048, y el número de subprocesos varió de 1 a 40 (40 carreras en total). Este es el gráfico con el tiempo de ejecución de cada ejecución (segundos x número de subprocesos):

La sorpresa fue la curva logarítmica, así que pensé que como la línea de trabajo era tan grande, los cachés compartidos no se utilizaban muy bien, así que desenterré este archivo virtual /sys/devices/system/cpu/cpu0/cache/index0/coherency_line_size que me dijo que el caché L1 de este procesador sincroniza las actualizaciones en grupos de 64 bytes (8 se duplica en el sol matriz). Así que puse el work_line a 8:

Entonces pensé que 8 era demasiado bajo para evitar las paradas de NUMA y establecer work_line en 16:

Mientras ejecutaba lo anterior, pensé "¿Quién soy yo para predecir qué work_line es buena? Veamos ...", y programada para ejecutar todas las work_line de work_line de 8 a 2048, pasos de 8 (es decir, cada múltiplo de la línea de caché, desde 1 a 256). Los resultados para 20 y 40 subprocesos (segundos x tamaño de la división del medio for bucle, divididos entre los subprocesos):

Creo que los casos con una work_line baja sufren una mala sincronización de la caché, mientras que una work_line mayor no ofrece beneficios más allá de un cierto número de subprocesos (supongo que porque la ruta de la memoria es el cuello de botella). Es muy triste que un problema que parece 100% paralelo presenta un comportamiento tan malo en una máquina real. Entonces, antes de convencerme de que los sistemas multi-core son una mentira muy bien vendida, les pregunto aquí primero:

¿Cómo puedo hacer que este código se amplíe linealmente al número de núcleos? ¿Qué me estoy perdiendo? ¿Hay algo en el problema que hace que no sea tan bueno como parece al principio?

Actualizar

Siguiendo las sugerencias, probé tanto dynamic programación static como la dynamic , pero eliminando los atómicos de lectura / escritura en el sol matriz. Para referencia, las líneas azul y naranja son las mismas del gráfico anterior (solo hasta work_line = 248; ). Las líneas amarillas y verdes son las nuevas. Por lo que pude ver: la static hace una diferencia significativa para la línea de trabajo baja, pero después de 96 los beneficios de la dynamic superan sus costos generales, haciéndolo más rápido. Las operaciones atómicas no hacen ninguna diferencia.

Answer 1

Incluso si no tiene un bloqueo mutex explícito en su código, tiene un recurso compartido entre sus procesos: la memoria y su bus. No ve esto en su código porque es el hardware el que se encarga de manejar todas las diferentes solicitudes de las CPU, pero sin embargo, es un recurso compartido.

Por lo tanto, cada vez que uno de sus procesos escribe en la memoria, la ubicación de la memoria tendrá que volver a cargarse desde la memoria principal por todos los demás procesos que la utilizan, y todos tienen que usar el mismo bus de memoria para hacerlo. El bus de memoria se satura, y no tiene más ganancia de rendimiento de los núcleos de CPU adicionales que solo sirven para empeorar la situación.

Answer 2

Intente ejecutar el IPCM ( Intel Performance Counter Monitor ). Puede ver el ancho de banda de la memoria y ver si se maximiza con más núcleos. Mi intuición es que tienes un ancho de banda de memoria limitado.

Como parte posterior rápida del cálculo de la envoltura, encuentro que el ancho de banda de lectura sin caché es de aproximadamente 10 GB / s en un Xeon. Si su reloj es de 2,5 GHz, esa es una palabra de 32 bits por ciclo de reloj. Su bucle interno es básicamente una operación de adición múltiple cuyos ciclos puede contar con una mano, más algunos ciclos para la sobrecarga del bucle. No me sorprende que después de 10 subprocesos, no obtengas ninguna ganancia de rendimiento.

Answer 3

La multiplicación de vector de matriz dispersa está vinculada a la memoria (consulte here ) y se puede mostrar con un modelo de línea de techo simple. Los problemas relacionados con la memoria se benefician del mayor ancho de banda de la memoria de los sistemas NUMA de múltiples tomas, pero solo si la inicialización de los datos se realiza de tal manera que los datos se distribuyen entre los dos dominios NUMA. Tengo algunas razones para creer que está cargando la matriz en serie y, por lo tanto, toda su memoria está asignada en un solo nodo NUMA. En ese caso, no se beneficiará del doble ancho de banda de memoria disponible en un sistema de doble zócalo y realmente no importa si usa la schedule(dynamic) o la schedule(static) . Lo que podría hacer es habilitar la política NUMA de intercalación de memoria para que la asignación de memoria se distribuya entre los dos nodos NUMA. Por lo tanto, cada subproceso terminará con un 50% de acceso a la memoria local y un 50% de acceso a la memoria remota en lugar de tener todos los subprocesos de la segunda CPU afectados por el 100% de acceso a la memoria remota. La forma más fácil de habilitar la política es mediante numactl :

$ OMP_NUM_THREADS=... OMP_PROC_BIND=1 numactl --interleave=all ./program ...

OMP_PROC_BIND=1 habilita la fijación de hilos y debería mejorar un poco el rendimiento.

También me gustaría señalar que esto:

done = true;
#pragma omp parallel shared(done)
{
    bool tdone = true;

    // ...

    #pragma omp atomic update
    done &= tdone;
}

Es probablemente una reimplementación no muy eficiente de:

done = true;
#pragma omp parallel reduction(&:done)
{
    // ...
        if(residual > tolerance) {
            done = false;
        }
    // ...
}

No tendrá una diferencia de rendimiento notable entre las dos implementaciones debido a la cantidad de trabajo realizado en el bucle interno, pero aún así no es una buena idea volver a implementar las primitivas de OpenMP existentes en aras de la portabilidad y la legibilidad.

Answer 4

Sospecho que tienes problemas de almacenamiento en caché. Cuando un hilo actualiza un valor en la matriz sol , invalida los cachés en otras CPU que almacenan esa misma línea de caché. Esto obliga a que se actualicen los cachés, lo que a su vez conduce a la parada de la CPU.

Answer 5

Su bucle interno tiene una omp atomic read , y su bucle medio tiene una omp atomic write en una ubicación que podría ser la misma que leyó una de las lecturas. OpenMP está obligado a garantizar que las escrituras atómicas y las lecturas de la misma ubicación se serialicen, por lo que de hecho es probable que sea necesario introducir un bloqueo, aunque no haya ningún explícito.

Incluso podría necesitar bloquear toda la matriz de sol menos que de alguna manera pueda averiguar qué lecturas pueden entrar en conflicto con las escrituras y, en realidad, los procesadores OpenMP no son necesariamente tan inteligentes.

Ningún código se escala de forma absolutamente lineal, pero puede estar seguro de que hay muchos códigos que se escalan mucho más cerca que de forma lineal.