La Multiplication matricielle OpenMP C++ s'exécute plus lentement en parallèle

j'apprends les bases de l'exécution parallèle de for loop en utilisant OpenMP.

malheureusement, mon programme paralel tourne 10 fois plus lentement que la version série. Ce que je fais mal? Je suis pas certains obstacles?

double **basicMultiply(double **A, double **B, int size) {
   int i, j, k;
   double **res = createMatrix(size);
   omp_set_num_threads(4);
   #pragma omp parallel for private(k)
   for (i = 0; i < size; i++) {
      for (j = 0; j < size; j++) {
         for (k = 0; k < size; k++) {
            res[i][j] += A[i][k] * B[k][j];
         }
      }
   }
   return res;
}

Merci beaucoup!

13
demandé sur rerx 2014-03-25 16:07:04

3 réponses

votre problème est dû à une condition de course sur la variable de boucle intérieure j. Il faut que ce soit privé.

Pour C89 je voudrais faire quelque chose comme ceci:

#pragma omp parallel
{
    int i, j, k;
    #pragma omp for
    for(i=0; ...

pour C++ ou C99 utiliser des déclarations mixtes

#pragma omp parallel for
for(int i=0; ...

en Faisant cela, vous n'avez pas à déclarer explicitement quelque chose de partagé ou privé.

quelques autres commentaires à votre code. Votre code threadé simple n'est pas facile à mettre en cache quand vous faites B[k][j]. Il lit un cacheline alors se déplace vers la ligne de cache suivante et ainsi de suite jusqu'à ce que le produit de point est fait à ce moment-là les autres lignes de cache ont été expulsées. Au lieu de cela, vous devriez prendre la transposition premier et accès BT[j][k]. De plus, vous avez attribué des tableaux de tableaux et non un tableau 2D contigu. J'ai corrigé votre code pour utiliser la transposition et un tableau 2D contigu.

Voici les temps que j'obtiens pour la taille=512.

no transpose  no openmp 0.94s
no transpose, openmp    0.23s
tranpose, no openmp     0.27s
transpose, openmp       0.08s

ci-dessous se trouve le code (Voir aussi http://coliru.stacked-crooked.com/a/ee174916fa035f97)

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

void transpose(double *A, double *B, int n) {
    int i,j;
    for(i=0; i<n; i++) {
        for(j=0; j<n; j++) {
            B[j*n+i] = A[i*n+j];
        }
    }
}

void gemm(double *A, double *B, double *C, int n) 
{   
    int i, j, k;
    for (i = 0; i < n; i++) { 
        for (j = 0; j < n; j++) {
            double dot  = 0;
            for (k = 0; k < n; k++) {
                dot += A[i*n+k]*B[k*n+j];
            } 
            C[i*n+j ] = dot;
        }
    }
}

void gemm_omp(double *A, double *B, double *C, int n) 
{   
    #pragma omp parallel
    {
        int i, j, k;
        #pragma omp for
        for (i = 0; i < n; i++) { 
            for (j = 0; j < n; j++) {
                double dot  = 0;
                for (k = 0; k < n; k++) {
                    dot += A[i*n+k]*B[k*n+j];
                } 
                C[i*n+j ] = dot;
            }
        }

    }
}

void gemmT(double *A, double *B, double *C, int n) 
{   
    int i, j, k;
    double *B2;
    B2 = (double*)malloc(sizeof(double)*n*n);
    transpose(B,B2, n);
    for (i = 0; i < n; i++) { 
        for (j = 0; j < n; j++) {
            double dot  = 0;
            for (k = 0; k < n; k++) {
                dot += A[i*n+k]*B2[j*n+k];
            } 
            C[i*n+j ] = dot;
        }
    }
    free(B2);
}

void gemmT_omp(double *A, double *B, double *C, int n) 
{   
    double *B2;
    B2 = (double*)malloc(sizeof(double)*n*n);
    transpose(B,B2, n);
    #pragma omp parallel
    {
        int i, j, k;
        #pragma omp for
        for (i = 0; i < n; i++) { 
            for (j = 0; j < n; j++) {
                double dot  = 0;
                for (k = 0; k < n; k++) {
                    dot += A[i*n+k]*B2[j*n+k];
                } 
                C[i*n+j ] = dot;
            }
        }

    }
    free(B2);
}

int main() {
    int i, n;
    double *A, *B, *C, dtime;

    n=512;
    A = (double*)malloc(sizeof(double)*n*n);
    B = (double*)malloc(sizeof(double)*n*n);
    C = (double*)malloc(sizeof(double)*n*n);
    for(i=0; i<n*n; i++) { A[i] = rand()/RAND_MAX; B[i] = rand()/RAND_MAX;}

    dtime = omp_get_wtime();
    gemm(A,B,C, n);
    dtime = omp_get_wtime() - dtime;
    printf("%f\n", dtime);

    dtime = omp_get_wtime();
    gemm_omp(A,B,C, n);
    dtime = omp_get_wtime() - dtime;
    printf("%f\n", dtime);

    dtime = omp_get_wtime();
    gemmT(A,B,C, n);
    dtime = omp_get_wtime() - dtime;
    printf("%f\n", dtime);

    dtime = omp_get_wtime();
    gemmT_omp(A,B,C, n);
    dtime = omp_get_wtime() - dtime;
    printf("%f\n", dtime);

    return 0;

}
25
répondu Z boson 2014-04-23 07:17:22

En plus. "Z boson", j'ai testé votre code C sur l'ordinateur portable avec intel i5 (2 noyaux physiques ou 4 logiques). Malheureusement, la vitesse de calcul n'est pas très rapide. 2000x2000 double aléatoire matrices j'ai obtenu les résultats suivants (à l'aide de VS 2010 avec OpenMP 2.0):

compilé pour Win64: C = A*B, où A,B sont des matrices avec la taille (2000x2000):

max number of threads = 4
Create random matrices:      = 0.303555 s
no transpose  no openmp = 100.539924 s
no transpose, openmp = 47.876084 s
transpose, no openmp = 27.872169 s
transpose, openmp = 15.821010 s

compilé pour Win32: C = A*B, où A, B sont des matrices avec la taille (2000x2000):

max number of threads = 4
Create random matrices:      = 0.378804 s
no transpose  no openmp = 98.613992 s
no transpose, openmp = 48.233655 s
transpose, no openmp = 29.590350 s
transpose, openmp = 13.678097 s

notez que pour le code" Hynek Blaha " le temps de calcul sur mon système est 739.208 s (226.62 s avec openMP)!

alors que dans Matlab x64:

n = 2000; 
A = rand(n); B = rand(n);

tic
C = A*B;
toc

le temps de calcul est 0.591440 secondes.

Mais en utilisant openBLAS le paquet, j'ai atteint une vitesse de 0.377814 secondes (utilisant minGW avec openMP 4.0). Tatou package fournit un moyen simple (à mon avis) pour la connexion des opérations matricielles avec openBLAS (ou d'autres paquets similaires). Dans ce cas, le code est

#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;

int main(){
    int n = 2000;
    N = 10; // number of repetitions
    wall_clock timer;

    arma_rng::set_seed_random();

    mat A(n, n, fill::randu), B(n, n, fill::randu);

    timer.tic();
    // repeat simulation N times
    for(int n=1;n<N;n++){
      mat C = A*B;
    }
    cout << timer.toc()/double(N) << "s" << endl;

    return 0;
}
4
répondu Alexander Korovin 2017-12-20 16:36:59

Si size est petit, le plafond de la synchronisation de thread va masquer n'importe quel gain de performance du calcul parallèle.

2
répondu rerx 2014-03-25 12:18:54