À propos du produit de matrices et de son optimisation...
Cas particulier d'un serveur Origin 2000

• pour plus d'information concernant la "Performance Tuning Optimization for Origin2000" rendez-vous ici.
• il s'agit ici d'optimiser le produit de deux matrices carrées. Nous étudions successivement plusieurs solutions séquentielles ou parallèles, avec ou sans option de compilation, avec ou sans les BLAS, avec ou sans OpenMP...,
• la commande runon numéroCpu ./a.out listeDesArguments permet de lancer l'exécutable de nom a.out du répertoire courant sur le processeur de numéro numéroCpu. Les processeurs de numéro 16 à 31 et 48 à 63 sont des processeurs rapides cadencés à 250 MHz, tous les autres (de 0 à 15 et de 32 à 47) sont cadencés à 195 MHz.
• comme le montre le tableau ci-après, le choix de la méthode de calcul n'est pas anodin puisque dans le cas d'un produit de matrices carrées à 2000x2000 éléments le même résultat peut-être obtenu en plus de 30 minutes ou en 5 secondes !

• pour obtenir des temps de calcul "optimaux" dans le cas d'une exécution en parallèle sur 8 processeurs, nous avons "réservé" un jeu complet de 8 cpu cadencées à 250 MHz sur le même module (le numéro 8 en l'occurence) à l'aide de la commande miser_cpuset.

1. Produit de matrice - le programme principal

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

extern void produitMatrice(double *A,double *B,double *C,int n);

void initialisation(double *T,int n)
{
  double nn = n*n;
  int i;
  for(i=0 ; i<nn ; i++) 
    T[i] = (double) rand() / (double) RAND_MAX;
}

double trace(double* T,int n)
{
  int i;
  double sommeTermesDiagonaux = 0;
  for(i=0 ; i<n ; i++)
    sommeTermesDiagonaux += T[i*(n+1)];
  return sommeTermesDiagonaux;
}

int main(int argc,char *argv[])
{
  clock_t dateDebut,dateFin;
  int n = atoi(argv[1]);  
  double *A = (double*) calloc(n*n,sizeof(double));
  double *B = (double*) calloc(n*n,sizeof(double));
  double *C = (double*) calloc(n*n,sizeof(double));

  initialisation(A,n);
  initialisation(B,n);
  
  dateDebut = clock();
  produitMatrice(A,B,C,n);
  dateFin = clock();
  
  printf("Execution en %5.2f secondes - Trace(C) = %.4f\n",
	 (double) (dateFin-dateDebut)/(double) CLOCKS_PER_SEC,trace(C,n));

  free(A);
  free(B);
  free(C);    
  return EXIT_SUCCESS;
}

2. Le produit de matrice : version simple avec ou sans optimisation à la compilation

void produitMatrice(double *A,double *B,double *C,int n)
{
  int i,j,k,in=0,kn;
  register double tmp;
  
  for(i=0 ; i<n ; i++) {
    for(j=0 ; j<n ; j++) {
      tmp = 0;
      for(k=0 , kn = 0 ; k<n ; k++ , kn += n)
        tmp += A[in+k]*B[kn+j];
      C[in+j] = tmp;
    }
    in += n;
  }
}

• résultats à l'exécution :

% . /opt/modules/modules/init/sh
% module load modules MIPSpro.7311 scsl
% cc produitMatricesMain.c produitMatricesSimple.c
% runon 10 ./a.out 2000

% runon 52 ./a.out 2000
Execution en 2061.73 secondes - Trace(C) = 999866.6307
% cc produitMatricesMain.c produitMatricesSimple.c -Ofast=ip27
% runon 10 ./a.out 2000
Execution en 63.13 secondes - Trace(C) = 999866.6307
% runon 52 ./a.out 2000
Execution en 48.63 secondes - Trace(C) = 999866.6307

• Conclusions :
• l'ajout de l'option -Ofast=ip27 (IP27 est le type générique des cartes nodales de bacchus) à la compilation permet d'accélérer l'exécution du code d'un facteur 37 dans le cas d'un produit de matrices carrées de taille 2000 !
• il est préférable de travailler sur un processeur cadencé à 250 MHz (accélération de l'ordre de 1.3 par rapport à un processeur cadencé à 195 MHz).


3. Le produit de matrice : version OpenMP

#include <omp.h>

void produitMatrice(double *A,double *B,double *C,int n)
{
  int i,j,k,in=0,kn;
  double tmp;
  /* ATTENTION : ICI ON A MODIFIE LE CALCUL DE LA VALEUR DE LA VARIABLE "in" ! */
#pragma omp parallel for private(j,k,in,kn,tmp) shared(A,B,C)
  for(i=0 ; i<n ; i++) {
    in = i*n; 
    for(j=0 ; j<n ; j++) {
      tmp = 0;
      for(k=0 , kn = 0 ; k<n ; k++ , kn += n)
        tmp += A[in+k]*B[kn+j];
      C[in+j] = tmp;
    }
  }
}

• résultats à l'exécution :

% cc produitMatricesMain.c produitMatricesOMP.c -mp
% OMP_NUM_THREADS=8 ./a.out 2000
Execution en 313.21 secondes - Trace(C) = 999866.6307
% cc produitMatricesMain.c produitMatricesOMP.c -mp -Ofast=ip27
% OMP_NUM_THREADS=8 ./a.out 2000
Execution en  7.59 secondes - Trace(C) = 999866.6307

• Conclusions :
• sans l'optimisation par l'option -Ofast=ip27 ajoutée à la compilation, la version OpenMP sur 8 processeurs n'accélère la version séquentielle que d'un facteur 6 !


4. Le produit de matrice : version CBLAS

#include <cblas.h>

void produitMatrice(double *A,double *B,double *C,int n)
{
  /* Utilisation des BLAS de niveau 3 - man dgemm pour plus de details */
  double alpha = 1., beta = 0.;
  cblas_dgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,n,n,n,alpha,A,n,B,n,beta,C,n);
}

• résultats à l'exécution :

% cc produitMatricesMain.c produitMatricesCBLAS.c -lblas
% runon 40 ./a.out 2000
Execution en 49.54 secondes - Trace(C) = 999866.6307
% runon 62 ./a.out 2000
Execution en 37.84 secondes - Trace(C) = 999866.6307
% cc produitMatricesMain.c produitMatricesCBLAS.c -lblas_mp
% OMP_NUM_THREADS=8 ./a.out 2000
Execution en  5.08 secondes - Trace(C) = 999866.6307

• Conclusions :
• utilisez les BLAS (ou les CBLAS en l'occurence) pour obtenir un code performant !
• pour plus de performance encore, je vous invite à utiliser les BLAS (ou CBLAS) parallèles, en précisant -lblas_mp (et non simplement -lblas) à l'édition de liens.