Boost : multiplications de matrice 1000*1000 lentes

**bossonm** · 08/07/2008, 11h55

Bonjour.
Je dois développer un prog qui fait du calcul sur des matrices le plus performants possibles.
Je suis donc à la recherche des meilleurs librairies pour ca.
Boost me semblait être une bonne solution, en tout cas simple d'utilisation en comparaison des librairies fortran (blas, lapack).
En testant simplement la multiplication de matrice, je me suis rendu compte que cela est très lent. Ma question est donc est-ce que j'ai oublié quelque chose ou ces performances sont normales.

Code c++ : 7 à 8 secondes

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 
#define NDEBUG
#include <stdio.h>
#include <math.h>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/timer.hpp> 
 
 
using namespace boost::numeric::ublas;
using namespace std;
 
typedef  matrix <double> Matrix;
boost::timer t; 
 
int main () {
    Matrix m(1000, 1000);
    Matrix r(1000,1000);
    for (unsigned i = 0; i < m.size1 (); ++ i){
        for (unsigned j = 0; j < m.size2 (); ++ j){
            m (i, j) = rand();
		}
    }
 
   t.restart(); 
   r=prod(m,m);
   cout<< t.elapsed(); 
 
   return 0;
}

Compilation : "c++ -I /Developer/SDKs/MacOSX10.5.sdk/usr/include main.cpp -o test -O3"

Code Matlab : 0.3 seconde environ

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
function r=test_mult(m,n)
   tic;
   r=m*n;
   toc

Merci de votre aide pour obtenir des performances similaire à Matlab

**pyknite** · 08/07/2008, 12h42

Hello...

Bon je vais peut etre pas aider, mais bon on sait jamais

Déjà la différence entre ton code c++ et matlab, c'est normal... Matlab est fait expès pour les math et calcul en tout genre... Ils ont des algo hyper pointu (et surtout très secret

)... En plus, la faut te dire, que tu lui demande pas grand chose a ton programme... Imagine si tu devais faire la multiplication des deux matrice

Autrement, je sais que tu peux remplacer ca:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
 
for (unsigned i = 0; i < m.size1 (); ++ i){
        for (unsigned j = 0; j < m.size2 (); ++ j){
            m (i, j) = rand();
		}
    }

par ca (mais tu vas pas vraiment gagner bcp):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
 
unsigned int Taille1 = m.size1();
unsigned int Taille2 = m.size2();
 
for (register unsigned int i = 0; i <= Taille1; ++ i)
     for (register unsigned int j = 0; j <= Taille2; ++ j)
            m (i, j) = rand();

déjà tu enregistre les taille dans des var, ca evite a chaque qu'il doivent recalculer la taille... ensuite, le register unsigned int permet d'utiliser directement un registre plutot qu'une variable..

Mais bon, comme je l'ai dit, c'est un peu du chipottage... ca va pas changer grand chose

Ahh oui, et quand tu n'as qu'une seul instruction dans un "sous-bloc" (genre if, while, for,...) pas besoin d'accolade

**bossonm** · 08/07/2008, 14h03

Merci pour les détailles d'optimisation. Ca peut toujours servir car je suis plutôt débutant en c++.
Mais mon problème est réellement au niveau des opérations sur matrice.
Je croyais que u-blas était plus ou moin équivalent à BLAS au niveau des performances. Forcément un peu moin bien, mais pas de l'ordre de 20 fois moin bien.

En fait je crois que Matlab utilise directement BLAS aussi dans le cas de 2 matrices pleines.

J'ai vérifié, en faisant directement appel à BLAS :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
t.restart(); 
cblas_dgemm ( CblasRowMajor, CblasNoTrans, CblasNoTrans, 1000, 1000, 1000, alpha, A, 1000, B, 1000, beta, C, 1000 );
cout<< t.elapsed();

Ca mes aussi environ 0.3 secondes (similaire à Matlab).

Mon problème, c'est que, par la suite je dois faire tourner mon prog sur un cluster avec 1.000.000 d'entrées différentes. Donc pour moi c'était bien pratique d'avoir un truc aussi simple à installer et à compiler que la librairie Boost.
Car je ne sais pas ce qu'ils ont installé sur leur cluster, l'installation et la compilation ont l'air fastidieuse avec BLAS (déjà installer dans mon macbook avec xcode).

**pyknite** · 08/07/2008, 17h55

autrement regarde un peu sur le net... Il y a des algo qui permettent d'optimiser certaines opérations sur les matrice... Et au lieu d'optimiser le code, tu optimise ton algo

(et c'est le plus efficace en général )

**bossonm** · 08/07/2008, 21h35

Oui bien sur je peux réimplémenter des algorithmes comme celui de Volker Strassen mais je pensais que c'était le rôle de boost de fournir un algo optimal. D'autant plus que ca permet de paralleliser.
Je voulais essentiellement savoir si c'était normal d'avoir des perf 20 fois inférieurs à BLAS.

**loufoque** · 09/07/2008, 04h45

Boost.ublas peut être configuré pour utiliser des bibliothèques de blas.

**Charlemagne** · 09/07/2008, 17h01

L'une des implémentation BLAS/LAPACK les plus rapides est celle d'ATLAS
http://math-atlas.sourceforge.net/

Ma bibliothèque est également très rapide mais je n'ai implémenté qu'une partie des fonctions BLAS, dont la multiplication matricielle.
J'ai mis des benchs en ligne, même des benchs en multi-thread.
http://www.ient.rwth-aachen.de/team/...al/genial.html.

PS: C'est normal qu'une implémentation "naïve" de la multiplication matricielle soit beaucoup plus lente que Matlab qui il me semble a abandonné sa propre implémentation pour y intégrer celle d'ATLAS (de même avec la fft de FFTW).
Que l'implémentation de Boost soit 20 fois plus lente montre déjà qu'il y a eu un effort d'optimisation, car sinon le rapport aurait vraisemblablement plutôt été de l'ordre de 100, si si)

PS2: Je ne connais pas vraiment les algos comme Volker Strassen (c'est bien un algo de l'ordre n^2.x, non?), mais je crois bien qu'en pratique ces algos ne valent pas grand chose pour la vitesse. Ces algos ne permettent pas d'utiliser des instructions SIMD, et surtout ne permettent pas d'optimiser les acces mémoire.

**Matthieu Brucher** · 10/07/2008, 08h13

Envoyé par pyknite

Matlab est fait expès pour les math et calcul en tout genre... Ils ont des algo hyper pointu (et surtout très secret

)... En plus, la faut te dire, que tu lui demande pas grand chose a ton programme... Imagine si tu devais faire la multiplication des deux matrice

Ou pas.
Matlab délègue ce genre de chose à une bibliothèque du commerce que eux ne développent pas (ACML, MKL ou anciennement ATLAS, FFTW pour la FFT, ...). Et pour le reste, ce sont des algorithmes standard trouvés dans la littérature.
Hyper pointus ? potentiellement, selon les toolbox. Très secrets ? Non.

Envoyé par bossonm

Oui bien sur je peux réimplémenter des algorithmes comme celui de Volker Strassen mais je pensais que c'était le rôle de boost de fournir un algo optimal. D'autant plus que ca permet de paralleliser.
Je voulais essentiellement savoir si c'était normal d'avoir des perf 20 fois inférieurs à BLAS.

Boost ne fournit aucun algorithme et surtout pas les algorithmes un peu compliqués.

**pyknite** · 10/07/2008, 11h09

Envoyé par Matthieu Brucher

Ou pas.
Matlab délègue ce genre de chose à une bibliothèque du commerce que eux ne développent pas (ACML, MKL ou anciennement ATLAS, FFTW pour la FFT, ...). Et pour le reste, ce sont des algorithmes standard trouvés dans la littérature.
Hyper pointus ? potentiellement, selon les toolbox. Très secrets ? Non.

Je ne connais pas très bien matlab, je me suis peut etre avancé un peu trop vite

merci pour les precision

**bossonm** · 11/07/2008, 15h29

Merci pour les précisions.
Ayant une formation surtout orienté math, je ne me rendais pas compte à quel point l'optimisation du code à un aussi bas niveau de programmation pouvait transfigurer les performances.

Peut-on forcer certaines variables à être stocker en mémoire cache, ou c'est la façon de programmer qui va entraîner une meilleur utilisation de la mémoire cache ?

De même pour les utilisations des instructions spéciales, c'est le programmeur qui précise à la main quels instructions utiliser ou la procédure est un minimum automatique. Il me semble que gcc essaye de produire des instructions optimisés mais peut-être faut-il l'aider un peu ?

Je suis totalement ignorant dans ce genre de problème. Merci d'éclairer un peu ma lanterne.

**bredelet** · 11/07/2008, 16h05

Est-ce que le C++ est plus rapide si tu utilises d'autres donnees dans ta matrice?

rand() (qui par ailleurs ne devrait pas etre utilise) renvoie un entier, alors tu multiplies des valeurs genre 34245545 * 1034345630 pour le moment.

D'autre part tu as:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

   r=prod(m,m);

Donc ce n'est pas necessaire d'allouer de la memoire d'avance pour r, puisque tu lui donne une nouvelle valeur ici. Tu peux juste le declarer:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

   Matrix r;

Il se peut que C++ perde du temps a allouer la memoire pour le resultat de la multiplication qui fait 1000*1000*8 octets. Est-ce que Boost n'a pas de produit qui donne un resultat "sur-place", comme ca tu peux reserver la memoire a l'avance?

**Charlemagne** · 11/07/2008, 17h01

Envoyé par bossonm

Peut-on forcer certaines variables à être stocker en mémoire cache, ou c'est la façon de programmer qui va entraîner une meilleur utilisation de la mémoire cache ?

Pas vraiment possible de forcer le stockage de certaines variables en mémoire cache. Dans l'ensemble, c'est le processeur qui retient lui-même en cache suivant diverse implémentations la mémoire récemment utilisée. Ceci dit, quand on connait les principes, il est possible d'adapter les algos pour travailler par blocs. (c'est le principe des algo d'ATLAS ou de ma bibliothèque)
Certes il existe 3-4 instructions pour aider le processeur à prendre les bonnes décisions (prefetch), mais en pratique ça ne vaut pas un bon algo.

Envoyé par bossonm

De même pour les utilisations des instructions spéciales, c'est le programmeur qui précise à la main quels instructions utiliser ou la procédure est un minimum automatique. Il me semble que gcc essaye de produire des instructions optimisés mais peut-être faut-il l'aider un peu ?

C'est la question qui revient tout le temps pour ceux qui sont étonnés de ne pas voir de gain de perfs quand ils indiquent aux compilos d'utiliser les instructions SIMD. Je considère ces options comme commerciales. Il y a de nombreuses raisons qui les rendent inéfficaces.
Les seules possibilités efficaces que je connaisse pour utiliser les instructions vectorielles:
- assembleur: ce que fait plus ou moins directement ATLAS, c'est pour ça qu'ATLAS ne se laisse compiler qu'avec GCC
- intrinsic functions: fonctions de bas niveau

**ac_wingless** · 16/07/2008, 14h39

L'implémentation de prod (et des autres produits) dans Boost est délibérément naïve (cf. code, r(i,j)+=m(i,k)*m(k,j)), avec en plus une légère perte due à l'utilisation de foncteurs. En revanche, il n'y a pas de copie inutile de temporaires si on utilise la 2ème version de prod: prod(m,m,r,true).

**bossonm** · 16/07/2008, 17h03

Merci pour les détails.

**Lucien63** · 10/08/2008, 15h58

J'ai testé avec de vulguaire tableau C de double.

Le temps tombe à 31ms sur un P4 2GHz

**ac_wingless** · 11/08/2008, 11h23

Envoyé par Lucien63

J'ai testé avec de vulguaire tableau C de double.

Le temps tombe à 31ms sur un P4 2GHz

J'avoue être assez curieux d'une solution allant 10 fois plus vite que les solutions industrielles classiques. Peut-on voir le code?

**Lucien63** · 11/08/2008, 11h59

J'avoue être assez curieux d'une solution allant 10 fois plus vite que les solutions industrielles classiques. Peut-on voir le code?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include<iostream>
#include<windows.h> // Pour GetTickCount
 
using namespace std;
 
int main(){
 
  const unsigned taille = 1000;
  int temps;
 
  //Allocation des deux tableaux dynamiquement.
  double **m = new double* [taille];
  for ( int i = 0; i < taille; ++i)    m[i] = new double[taille];
 
  double **r = new double* [taille];
  for ( int i = 0; i < taille; ++i)    r[i] = new double[taille];
 
  //Remplissage de m avec rand comme dans l'exemple.
  for (unsigned i = 0; i < taille; ++ i){
        for (unsigned j = 0; j < taille; ++ j){
            m[i][j] = rand();
		}
    }
 
  temps = GetTickCount();   //Pour mesurer le temps
 
  //Produit des cellules
  for (unsigned i = 0; i < taille; ++ i){
        for (unsigned j = 0; j < taille; ++ j){
            r[i][j] = m[i][j]*m[i][j];
		}
    }
 
  temps = GetTickCount() - temps;
  cout << temps << endl;
 
 
//Libération mémoire
  for (int i = 0; i < taille; ++i)   delete [] r[i] ;
  delete [] r;
  for (int i = 0; i < taille; ++i)   delete [] m[i] ;
  delete [] m;
 
 
  return 0;
}

Avec des tableaux de 2000 x 2000 ---> 93 ms
(les mesures en dessous de 50ms avec GetTickCount ne sont pas terribles)

**Matthieu Brucher** · 11/08/2008, 12h20

Mais c'est pas une multiplication matricielle, ça...

**Lucien63** · 11/08/2008, 15h03

C'est peut-être plutôt ça :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include<iostream>
#include<windows.h> // Pour GetTickCount
 
using namespace std;
 
int main(){
 
  const unsigned taille = 1000;
  int temps;
 
  //Allocation des deux tableaux dynamiquement.
  double **m = new double* [taille];
  for ( int i = 0; i < taille; ++i)    m[i] = new double[taille];
 
  double **r = new double* [taille];
  for ( int i = 0; i < taille; ++i)    r[i] = new double[taille];
 
  //Remplissage de m avec rand comme dans l'exemple.
  for (unsigned i = 0; i < taille; ++ i){
        for (unsigned j = 0; j < taille; ++ j){
            m[i][j] = rand();
		}
    }
 
  temps = GetTickCount();   //Pour mesurer le temps
 
  //Produit matriciel
  for (unsigned i = 0; i < taille; ++i){
     for (unsigned j = 0; j < taille; ++j){
        r[i][j]=0;
        for (unsigned n = 0; n < taille; ++n){
            r[i][j] += m[i][n]*m[n][j];
        }
     }
  }
 
  temps = GetTickCount() - temps;
  cout << temps << endl;
 
 
  //Libération mémoire
  for (int i = 0; i < taille; ++i)   delete [] r[i] ;
  delete [] r;
  for (int i = 0; i < taille; ++i)   delete [] m[i] ;
  delete [] m;
 
        return 0;
}