Pourquoi openmp 32 thread est beaucoup plus lent qu'un thread?

Dec 17 2020

J'essaie d'écrire une application calculant la norme l2 de 2 tableaux. Je dois mettre en parallèle mon calcul.

Voici le code que j'ai parallélisé:

  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);

Je compile le code comme:

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

J'exécute ce code avec 1 threads et 32 ​​threads. La sortie est exactement le contraire de ce qui est attendu. Voici un exemple de sortie:

[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

Est-ce que je vois mal ou utiliser 32 threads est 10 fois plus lent qu'un thread? Alors, qu'est-ce que je fais de mal ici?

Voici mon code complet:

#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;
}

Voici la sortie;


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

Réponses

2 dreamcrash Dec 17 2020 at 01:33

Est-ce que je vois mal ou utiliser 32 threads est 10 fois plus lent qu'un thread? Alors, qu'est-ce que je fais de mal ici?

Dans la partie du code qui est à la fois profilée et parallélisée avec 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);
    }
 }

il y a une condition de concurrence, à savoir l'accès à la variable l2_norm. De plus, vous pouvez supprimer le private(i), car la variable d' index ( c'est-à-dire i ) dans la boucle parallélisée sera définie implicitement comme privée par OpenMP. La condition de concurrence peut être corrigée avec la réduction OpenMP . De plus, votre boucle ne distribue pas réellement les itérations entre les threads comme vous le souhaitez. Parce que vous avez ajouté à nouveau la clause parallel à cela #pragma omp for, et en supposant que vous avez désactivé le parallélisme imbriqué, ce qui est le cas par défaut, chacun des threads créés à l'extérieur parallel regionexécutera "séquentiellement" le code dans cette région, à savoir:

    #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);
    }

Par conséquent, chaque thread exécutera toutes les Nitérations de la boucle que vous vouliez paralléliser. Par conséquent, en supprimant le parallélisme et en ajoutant une surcharge supplémentaire ( par exemple, la création de thread) au code séquentiel. Pour résoudre ces problèmes ( c.-à-d. Condition de concurrence et région parallèle "imbriquée" ), modifiez ce code en:

 #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);
    }
 }

Maintenant, après avoir résolu ces problèmes, vous vous retrouvez avec un autre problème (en termes de performances), à savoir que la boucle parallèle est effectuée dans le contexte d'une parallélisation hybride de OpenMP + MPI, et vous n'avez pas explicitement lié les OpenMPthreads (dans les MPIprocessus) à les noyaux correspondants. Sans cette liaison explicite, on ne peut pas être sûr dans quels cœurs ces threads aboutiront. Naturellement, le plus souvent, le fait d'avoir plusieurs threads s'exécutant dans le même cœur logique augmentera l'exécution globale de l'application parallélisée.

Si votre application utilise des threads, vous voulez probablement vous assurer que vous n'êtes pas du tout lié (en spécifiant --bind-to none), ou lié à plusieurs cœurs en utilisant un niveau de liaison approprié ou un nombre spécifique d'éléments de traitement par application traiter. Vous pouvez résoudre ce problème en:

  1. la désactivation de la liaison avec l'indicateur MPI --bind-to none, pour permettre aux threads d'être affectés à différents cœurs;
  2. ou effectuez la liaison des threads, en conséquence. Consultez ce thread SO pour savoir comment mapper les threads aux cœurs dans les parallélisations hybrides telles que MPI + OpenMP.

En définissant explicitement le nombre de threads par processus en conséquence, vous pouvez éviter que plusieurs threads se retrouvent dans le même noyau, et par conséquent, éviter que les threads du même noyau se battent pour les mêmes ressources.

Conseils:

OMI, vous devez d'abord tester les performances du OpenMPseul, sans aucun processus MPI. Dans ce contexte, tester l'évolutivité du code en mesurant la version séquentielle contre 2fils, puis 4, 8et ainsi de suite, en augmentant progressivement le nombre de threads. Finalement, il y aura un certain nombre de threads pour lesquels le code arrête simplement la mise à l'échelle. Naturellement, la quantité de travail parallèle effectuée par les filetages doit être suffisamment importante pour surmonter la surcharge du parallélisme. Par conséquent, vous devriez également tester avec des entrées de plus en plus grandes.

Après avoir profilé, testé et amélioré votre OpenMPversion, vous pouvez ensuite étendre cette parallélisation de la mémoire partagée avec plusieurs processus en utilisant MPI.

1 HristoIliev Dec 17 2020 at 19:49

Outre la condition de concurrence dans la mise à jour d'une variable partagée comme indiqué dans la réponse de @ dreamcrash, votre code ne distribue pas correctement le travail.

#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 parallelconstruction de la boucle interne en fait une forconstruction parallèle combinée imbriquée . Cela signifie que chaque thread de l'équipe exécutant la boucle parallèle externe génère une toute nouvelle région parallèle et distribue la iboucle sur les threads qu'elle contient. Il n'y a pas de distribution dans la région parallèle externe et vous vous retrouvez avec N threads qui répètent exactement le même travail. Par défaut, le parallélisme imbriqué est désactivé, de sorte que la région parallèle imbriquée s'exécute séquentiellement et votre code le fait effectivement:

#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);
    }
}

Il n'y a pas de distribution du travail et tous les threads écrivent aux mêmes emplacements dans le diff_vector[]tableau.

D'une part, ce code en général est lié à la mémoire car la quantité de calcul par octet de données est faible - les processeurs modernes peuvent effectuer de nombreuses multiplications et soustractions par cycle tout en extrayant des données de la mémoire et en y écrivant les résultats prend de nombreux cycles. Les problèmes liés à la mémoire ne sont pas plus rapides avec plus de threads puisque le facteur limitant est la bande passante de la mémoire. Ce n'est pas un problème si grave dans votre cas car les entrées de tableau 32K occupent 256 Ko de mémoire et cela tient dans la plupart des caches de processeur, et le cache L3 est extrêmement rapide, mais il est toujours plus grand que le cache L1 le plus rapide Noyau CPU. D'autre part, l'écriture dans les mêmes zones de mémoire à partir de plusieurs threads entraîne un partage vrai et faux, avec l'invalidation de cache inter-thread associée, ce qui entraîne généralement une exécution du code parallèle beaucoup plus lente que la version séquentielle.

Il existe des outils qui peuvent vous aider à analyser les performances de votre code et à repérer les problèmes. Comme je l'ai déjà écrit dans un commentaire, Intel VTune est l'un d'entre eux et est disponible gratuitement dans le cadre de la boîte à outils oneAPI. Intel Inspector en est un autre (encore une fois, gratuit et fait partie de la boîte à outils oneAPI) et il détecte des problèmes tels que des courses de données. Les deux outils fonctionnent très bien ensemble et je ne pourrais pas les recommander assez fortement à tout programmeur parallèle en herbe.

Il existe également une condition de concurrence mineure pour l'écriture numberOfThreads, mais comme toutes les valeurs écrites sont identiques, ce n'est pas vraiment un problème logique. La version correcte du code en question doit être:

#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);
    }
}