¿Por qué el hilo openmp 32 es mucho más lento que 1 hilo?

Question

Estoy tratando de escribir una aplicación que calcule la norma l2 de 2 matrices. Tengo que hacer un paralelo con mi cálculo.

Aquí está el código que he paralelizado:

  double time_start_openmp = omp_get_wtime();
  #pragma omp parallel for
  for (i = 0; i < n; i++)
  {
       numberOfThreads = omp_get_num_threads();
       double local_diff = x[i] - xseq[i];
       diff_vector[i] = local_diff;
       l2_norm += (local_diff * local_diff);
  }

   time_end_openmp = omp_get_wtime();

   l2_norm = sqrt(l2_norm);

   openmp_exec_time = time_end_openmp - time_start_openmp;
   printf("OPENMP: %d %ld %f %.12e\n", n, numberOfThreads, openmp_exec_time, l2_norm);

Compilo el código como:

gcc -fopenmp -g -ggdb -Wall -lm -o test test.c 

Estoy ejecutando este código con 1 subproceso y 32 subprocesos. El resultado es exactamente lo contrario de lo esperado. Aquí hay un resultado de ejemplo:

[hayri@hayri-durmaz MatrixMultipication_MPI]$ export OMP_NUM_THREADS=32 [hayri@hayri-durmaz MatrixMultipication_MPI]$ ./test 10000
OPENMP: 10000 32 0.001084 0.000000000000e+00
[hayri@hayri-durmaz MatrixMultipication_MPI]$ export OMP_NUM_THREADS=1 [hayri@hayri-durmaz MatrixMultipication_MPI]$ ./test 10000
OPENMP: 10000 1 0.000106 0.000000000000e+00

¿Estoy viendo mal o usar 32 subprocesos es 10 veces más lento que 1 subproceso? Entonces, ¿qué estoy haciendo mal aquí?

Aquí está mi código completo:

#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
#include <math.h>

#define MATSIZE 2000

static size_t totalMemUsage = 0;

size_t vectors_dot_prod(double *x, double *y, size_t n)
{
    double res = 0.0;
    size_t i;
    for (i = 0; i < n; i++)
    {
        res += x[i] * y[i];
    }
    return res;
}

size_t vectors_dot_prod2(double *x, double *y, size_t n)
{
    size_t res = 0.0;
    size_t i = 0;
    for (; i <= n - 4; i += 4)
    {
        res += (x[i] * y[i] +
                x[i + 1] * y[i + 1] +
                x[i + 2] * y[i + 2] +
                x[i + 3] * y[i + 3]);
    }
    for (; i < n; i++)
    {
        res += x[i] * y[i];
    }
    return res;
}

void matrix_vector_mult(double **mat, double *vec, double *result, size_t rows, size_t cols)
{ // in matrix form: result = mat * vec;
    size_t i;
    for (i = 0; i < rows; i++)
    {
        result[i] = vectors_dot_prod2(mat[i], vec, cols);
    }
}

double get_random()
{

    double range = 1000;
    double div = RAND_MAX / range;
    double randomNumber = (rand() / div);
    // printf("%d\n", randomNumber);
    return randomNumber;
}

void print_2d_arr(double *arr, size_t row, size_t col)
{
    size_t i, j, index;

    for (i = 0; i < row; i++)
    {
        for (j = 0; j < col; j++)
        {
            index = i * col + j;
            printf("%3f ", arr[index]);
        }
        printf("\n");
    }
}
void print_1d_arr(double *arr, size_t row)
{
    size_t i;
    for (i = 0; i < row; i++)
    {
        printf("%f, ", arr[i]);
    }
    printf("\n");
}

size_t **fullfillArrayWithRandomNumbers(double *arr, size_t n)
{
    /*
    * Fulfilling the array with random numbers 
    * */
    size_t i;
    for (i = 0; i < n; i++)
    {
        arr[i] = get_random();
    }
    return 0;
}

double *allocarray1D(size_t size)
{
    double *array = calloc(size, sizeof(double));
    totalMemUsage = totalMemUsage + size * sizeof(double);
    return array;
}

size_t ParallelRowMatrixVectorMultiply(size_t n, double *a, double *b, double *x, MPI_Comm comm)
{
    size_t i, j;
    size_t nlocal;
    double *fb;
    int npes, myrank;
    MPI_Comm_size(comm, &npes);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    fb = (double *)malloc(n * sizeof(double));
    nlocal = n / npes;
    MPI_Allgather(b, nlocal, MPI_DOUBLE, fb, nlocal, MPI_DOUBLE, comm);
    for (i = 0; i < nlocal; i++)
    {
        x[i] = 0.0;
        for (j = 0; j < n; j++)
        {
            size_t index = i * n + j;
            x[i] += a[index] * fb[j];
        }
    }
    free(fb);
    return 0;
}

size_t ParallelRowMatrixVectorMultiply_WithoutAllgather(size_t n, double *a, double *b, double *x_partial, double *x, MPI_Comm comm)
{

    // Process 0 sends b to everyone
    MPI_Bcast(b, n, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    size_t i, j;
    size_t nlocal;
    // double *fb;
    int npes, myrank;
    MPI_Comm_size(comm, &npes);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
    // fb = (double *)malloc(n * sizeof(double));
    nlocal = n / npes;
    // MPI_Allgather(b, nlocal, MPI_DOUBLE, fb, nlocal, MPI_DOUBLE, comm);
    for (i = 0; i < nlocal; i++)
    {
        x_partial[i] = 0.0;
        for (j = 0; j < n; j++)
        {
            size_t index = i * n + j;
            // printf("%f x %f\n", a[index], b[j]);
            x_partial[i] += a[index] * b[j];
        }
    }
    // free(b);

    // Process 0 gathers x_partials to create x
    MPI_Gather(x_partial, nlocal, MPI_DOUBLE, x, nlocal, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    return 0;
}

size_t SequentialMatrixMultiply(size_t n, double *a, double *b, double *x)
{
    size_t i, j;
    for (i = 0; i < n; i++)
    {
        x[i] = 0.0;
        for (j = 0; j < n; j++)
        {
            size_t index = i * n + j;
            // printf("%f x %f\n", a[index], b[j]);
            x[i] += a[index] * b[j];
        }
    }
    return 0;
}

int main(int argc, char *argv[])
{
    // Global declerations
    size_t i;
    // MPI_Status status;

    // Initialize the MPI environment
    MPI_Init(&argc, &argv);

    // Get the number of processes
    int world_size;
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);

    // Get the rank of the process
    int taskid;
    MPI_Comm_rank(MPI_COMM_WORLD, &taskid);

    // Get the name of the processor
    char processor_name[MPI_MAX_PROCESSOR_NAME];
    int name_len;
    MPI_Get_processor_name(processor_name, &name_len);

    if (argc != 2)
    {
        if (taskid == 0)
            printf("Usage: %s <N>\n", argv[0]);
        MPI_Finalize();
        return 0;
    }
    srand(time(NULL) + taskid);
    size_t n = atoi(argv[1]);
    size_t nOverK = n / world_size;

    double *a = allocarray1D(n * n);
    double *b = allocarray1D(n);
    double *x = allocarray1D(n);
    double *x_partial = allocarray1D(nOverK);
    double *xseq = allocarray1D(n);

    double *a_partial = allocarray1D(n * nOverK);

    if (a == NULL || b == NULL || x == NULL || xseq == NULL || x_partial == NULL)
    {
        if (taskid == 0)
            printf("Allocation failed\n");
        MPI_Finalize();
        return 0;
    }
    // Process 0 creates A matrix.
    if (taskid == 0)
    {
        fullfillArrayWithRandomNumbers(a, n * n);
        // Process 0 produces the b
        fullfillArrayWithRandomNumbers(b, n);
    }

    // Process 0 sends a_partial to everyone
    if (!(world_size == 1 && n == 64000))
    {
        MPI_Scatter(a, n * nOverK, MPI_DOUBLE, a_partial, n * nOverK, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    }

    MPI_Barrier(MPI_COMM_WORLD);
    double time_start = MPI_Wtime();
    ParallelRowMatrixVectorMultiply_WithoutAllgather(n, a_partial, b, x_partial, x, MPI_COMM_WORLD);
    double time_end = MPI_Wtime();
    double parallel_exec_time = time_end - time_start;

    double *exec_times = allocarray1D(world_size);
    // Process 0 gathers x_partials to create x
    MPI_Gather(&parallel_exec_time, 1, MPI_DOUBLE, exec_times, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    // print_1d_arr(x, n);

    if (taskid == 0)
    {
        SequentialMatrixMultiply(n, a, b, xseq);
        // check difference between x and xseq using OpenMP
        //print_1d_arr(exec_times, world_size);
        // print_1d_arr(xseq, n);
        double max_exec, min_exec, avg_exec;
        min_exec = 1000;
        for (i = 0; i < world_size; i++)
        {
            if (max_exec < exec_times[i])
            {
                max_exec = exec_times[i];
            }
            if (min_exec > exec_times[i])
            {
                min_exec = exec_times[i];
            }
            avg_exec += exec_times[i];
        }
        avg_exec = avg_exec / world_size;

        long double time_start_openmp = omp_get_wtime();
        long double time_end_openmp, openmp_exec_time, min_exec_time, max_exec_time, avg_exec_time;
        max_exec_time = 0;
        max_exec_time = 1000;
        long double l2_norm = 0;
        size_t numberOfThreads = 0;
        size_t r = 0;
        double *diff_vector = allocarray1D(n);
        size_t nrepeat = 10000;

        if (world_size == 1)
        {
            #pragma omp parallel
            {
                numberOfThreads = omp_get_num_threads();
                #pragma omp parallel for private(i)
                for (i = 0; i < n; i++)
                {
                    double local_diff = x[i] - xseq[i];
                    diff_vector[i] = local_diff;
                    l2_norm += (local_diff * local_diff);
                }
            }
        }
        else
        {
            #pragma omp parallel
            {
                numberOfThreads = omp_get_num_threads();
                #pragma omp parallel for private(i)
                for (i = 0; i < n; i++)
                {
                    double local_diff = x[i] - xseq[i];
                    diff_vector[i] = local_diff;
                    l2_norm += (local_diff * local_diff);
                }
            }
        }
        l2_norm = sqrt(l2_norm);
        time_end_openmp = omp_get_wtime();
        openmp_exec_time = time_end_openmp - time_start_openmp;
        // print matrix size, number of processors, number of threads, time, time_openmp, L2 norm of difference of x and xseq (use %.12e while printing norm)
        if (world_size == 1)
        {
            printf("OPENMP: %d %ld %Lf %.12e\n", n, numberOfThreads, openmp_exec_time, openmp_exec_time, l2_norm);
            printf("NEW_OPENMP: %d %ld %f %.12e\n", n, numberOfThreads, openmp_exec_time, l2_norm);
        }
        printf("MIN_AVG_MAX: %d %d %f %f %f\n", n, world_size, min_exec, max_exec, avg_exec);
        printf("MPI: %d %d %f %.12Lf %.12e\n", n, world_size, max_exec, l2_norm, l2_norm);
        totalMemUsage = totalMemUsage / (1024 * 1024 * 1024);
        printf("TOTALMEMUSAGE: %zu\n", totalMemUsage);

        //printf("process: %d %d %d %f %.12e\n", taskid, n, world_size, parallel_exec_time, l2_norm);
        //printf("%d %ld %f %.12e\n", n, numberOfThreads, openmp_exec_time, l2_norm);
    }
    MPI_Finalize();
    return 0;
}

Aquí está la salida;

cn009
36
mpicc -fopenmp -g -ggdb  -lm -o rowmv rowmv.c 


OPENMP: 32000 1 0.000299 2.991110086441e-04
MIN_AVG_MAX: 32000 1 3.112523 3.112523 3.112523
MPI: 32000 1 3.112523 0.000000000000 9.532824124368e-130
TOTALMEMUSAGE: 15


OPENMP: 32000 2 0.000535 5.350699648261e-04
MIN_AVG_MAX: 32000 1 3.125519 3.125519 3.125519
MPI: 32000 1 3.125519 0.000000000000 9.532824124368e-130
TOTALMEMUSAGE: 15


OPENMP: 32000 4 0.000434 4.341900348663e-04
MIN_AVG_MAX: 32000 1 3.170650 3.170650 3.170650
MPI: 32000 1 3.170650 0.000000000000 9.532824124368e-130
TOTALMEMUSAGE: 15


OPENMP: 32000 8 0.000454 4.542167298496e-04
MIN_AVG_MAX: 32000 1 3.168685 3.168685 3.168685
MPI: 32000 1 3.168685 0.000000000000 9.532824124368e-130
TOTALMEMUSAGE: 15


OPENMP: 32000 16 0.000507 5.065393634140e-04
MIN_AVG_MAX: 32000 1 3.158761 3.158761 3.158761
MPI: 32000 1 3.158761 0.000000000000 9.532824124368e-130
TOTALMEMUSAGE: 15


OPENMP: 32000 32 0.000875 8.752988651395e-04
MIN_AVG_MAX: 32000 1 3.166051 3.166051 3.166051
MPI: 32000 1 3.166051 0.000000000000 9.532824124368e-130
TOTALMEMUSAGE: 15

Answer 1

¿Estoy viendo mal o usar 32 subprocesos es 10 veces más lento que 1 subproceso? Entonces, ¿qué estoy haciendo mal aquí?

En la parte del código que se perfila y se paraleliza con OpenMP:

 #pragma omp parallel
 {
    numberOfThreads = omp_get_num_threads();
    #pragma omp parallel for private(i)
    for (i = 0; i < n; i++)
    {
        double local_diff = x[i] - xseq[i];
        diff_vector[i] = local_diff;
        l2_norm += (local_diff * local_diff);
    }
 }

hay una condición de carrera, es decir, el acceso a la variable l2_norm. Además, puede eliminar el private(i), ya que la variable de índice ( es decir, i ) en el bucle en paralelo se establecerá implícitamente como privada por OpenMP. La condición de carrera se puede arreglar con la reducción de OpenMP . Además, su bucle no está distribuyendo las iteraciones entre los subprocesos como deseaba. Debido a que agregó nuevamente la cláusula paralela a eso #pragma omp for, y asumiendo que ha desactivado el paralelismo anidado, que de manera predeterminada es, cada uno de los subprocesos creados en el exterior parallel regionejecutará "secuencialmente" el código dentro de esa región, a saber:

    #pragma omp parallel for private(i)
    for (i = 0; i < n; i++)
    {
        double local_diff = x[i] - xseq[i];
        diff_vector[i] = local_diff;
        l2_norm += (local_diff * local_diff);
    }

Por lo tanto, cada subproceso ejecutará todas las Niteraciones del bucle que pretendía que se paralelizara. En consecuencia, eliminar el paralelismo y agregar gastos generales adicionales ( por ejemplo, creación de subprocesos) al código secuencial. Para solucionar esos problemas ( es decir, condición de carrera y región paralela "anidada" ) cambie este código a:

 #pragma omp parallel
 {
    numberOfThreads = omp_get_num_threads();
    #pragma omp for reduction(+:l2_norm)
    for (i = 0; i < n; i++)
    {
        double local_diff = x[i] - xseq[i];
        diff_vector[i] = local_diff;
        l2_norm += (local_diff * local_diff);
    }
 }

Ahora, habiendo solucionado esos problemas, todavía le queda otro problema (en cuanto al rendimiento), a saber, que el bucle paralelo se está realizando en el contexto de una paralelización híbrida de OpenMP + MPI, y no ha vinculado explícitamente los OpenMPsubprocesos (dentro de los MPIprocesos) a los núcleos correspondientes. Sin esa vinculación explícita, no se puede estar seguro de en qué núcleos terminarán esos hilos. Naturalmente, la mayoría de las veces, tener varios subprocesos ejecutándose en el mismo núcleo lógico aumentará la ejecución general de la aplicación que se paraleliza.

Si su aplicación usa subprocesos, entonces probablemente desee asegurarse de que no está vinculado en absoluto (especificando --bind-to none), o vinculado a múltiples núcleos usando un nivel de vinculación apropiado o una cantidad específica de elementos de procesamiento por aplicación proceso. Puede resolver este problema de una de las siguientes maneras:

1. deshabilitar el enlace con la bandera MPI --bind-to none, para permitir que los subprocesos se asignen a diferentes núcleos;
2. o realizar el encuadernado de hilos, en consecuencia. Consulte este hilo SO sobre cómo asignar los hilos a los núcleos en paralelizaciones híbridas como MPI + OpenMP.

Al establecer explícitamente el número de subprocesos por proceso en consecuencia, puede evitar que varios subprocesos terminen en el mismo núcleo y, en consecuencia, evitar que los subprocesos dentro del mismo núcleo luchen por los mismos recursos.

Consejo:

En mi opinión, primero debe probar el rendimiento del OpenMPsolo, sin ningún proceso MPI. En este contexto, probar la escalabilidad de código mediante la medición de la versión secuencial contra 2las discusiones, a continuación 4, 8y así sucesivamente, aumentando gradualmente el número de hilos. Eventualmente, habrá una serie de subprocesos para los que el código simplemente dejará de escalar. Naturalmente, la cantidad de trabajo paralelo que realizan los subprocesos tiene que ser lo suficientemente grande para superar la sobrecarga del paralelismo. Por lo tanto, también debe probar con entradas cada vez más grandes.

Después de haber perfilado, probado y mejorado su OpenMPversión, puede extender esa paralelización de memoria compartida con múltiples procesos usando MPI.

Answer 2

Además de la condición de carrera al actualizar una variable compartida como se indica en la respuesta de @ dreamcrash, su código no distribuye el trabajo correctamente.

#pragma omp parallel
{
    numberOfThreads = omp_get_num_threads();
    #pragma omp parallel for private(i)
                ~~~~~~~~
    for (i = 0; i < n; i++)
    {
        double local_diff = x[i] - xseq[i];
        diff_vector[i] = local_diff;
        l2_norm += (local_diff * local_diff);
    }
}

La parallelconstrucción en el bucle interno lo convierte en una forconstrucción paralela combinada anidada . Significa que cada subproceso del equipo que ejecuta el bucle paralelo exterior genera una nueva región paralela y distribuye el ibucle entre los subprocesos que contiene. No hay distribución en la región paralela exterior y terminas con N subprocesos que repiten exactamente el mismo trabajo. De forma predeterminada, el paralelismo anidado está deshabilitado, por lo que la región paralela anidada se ejecuta secuencialmente y su código está haciendo esto de manera efectiva:

#pragma omp parallel
{
    numberOfThreads = omp_get_num_threads();
    for (i = 0; i < n; i++)
    {
        double local_diff = x[i] - xseq[i];
        diff_vector[i] = local_diff;
        l2_norm += (local_diff * local_diff);
    }
}

No hay distribución de trabajo y todos los subprocesos escriben en las mismas ubicaciones en la diff_vector[]matriz.

Por un lado, este código en general está limitado a la memoria, ya que la cantidad de cálculo por byte de datos es baja: las CPU modernas pueden hacer muchas multiplicaciones y restas por ciclo, mientras que obtener datos de la memoria y escribir los resultados allí lleva muchos ciclos. Los problemas vinculados a la memoria no se vuelven más rápidos con más subprocesos, ya que el factor limitante es el ancho de banda de la memoria. Esto no es un gran problema en su caso porque las entradas de la matriz de 32K ocupan 256 KB de memoria y eso cabe en la mayoría de los cachés de CPU, y el caché L3 es increíblemente rápido, pero aún es más grande que el caché L1 más rápido de un Núcleo de CPU. Por otro lado, escribir en las mismas áreas de memoria desde múltiples subprocesos da como resultado un intercambio verdadero y falso, con la invalidación de caché entre subprocesos asociada, que generalmente da como resultado que el código paralelo se ejecute mucho más lento que la versión secuencial.

Existen herramientas que pueden ayudarlo a analizar el rendimiento de su código y detectar problemas. Como ya escribí en un comentario, Intel VTune es uno de ellos y está disponible gratuitamente como parte del kit de herramientas oneAPI. Intel Inspector es otro (de nuevo, gratuito y parte del kit de herramientas oneAPI) y encuentra problemas como carreras de datos. Las dos herramientas funcionan muy bien juntas y no podría recomendarlas lo suficiente a ningún aspirante a programador paralelo.

También hay una condición de carrera menor en la que se escribe numberOfThreads, pero dado que todos los valores escritos son iguales, no es un gran problema lógico. La versión correcta del código en cuestión debería ser:

#pragma omp parallel
{
    #pragma omp master
    numberOfThreads = omp_get_num_threads();

    #pragma omp parallel reduction(+:l2_norm)
    for (i = 0; i < n; i++)
    {
        double local_diff = x[i] - xseq[i];
        diff_vector[i] = local_diff;
        l2_norm += (local_diff * local_diff);
    }
}