5. TP 5
Produit scalaire



Vous pouvez retrouver une vidéo de ce TP sur youtube.

5.1. Corps du programme

Créer un exécutable (dot_product.exe) en architecture 32 bits qui permet de calculer le produit scalaire de deux vecteurs de floats en utilisant les méthodes suivantes :

On rappelle que le produit scalaire de deux vecteurs est calculé par :

$$ ps(x,y) = ∑↙{i=0}↖{n-1} x_i × y_i $$

Le programme sera composé de codes sources en C++ et en assembleur :

Le programme devra prendre en entrée 3 paramètres :

 Description   Variable   Option (getopt) 
 nombre de répétitions du calcul   zillions   -z entier 
 taille des vecteurs de réels   vector_size   -s entier 
 méthode de calcul utilisée   method   -m entier 
Variables et paramètres en ligne de commande

On utilisera deux vecteurs de floats dont la taille sera multiple de 8 (car nous utiliserons la technologie vectorielle AVX) et qui seront alloués dynamiquement en utilisant _mm_malloc ou posix_memalign.

On calculera le nombre de cycles CPU utilisés par chaque version en utilisant la classe CPUTimer qui se compose de 2 fichiers :

Note : on peut également mesurer le temps en utilisant la librairie chrono du C++11 :

  1. #include <iostream>
  2. #include <chrono>
  3. #include <unistd.h>
  4.  
  5. using namespace std;
  6. using namespace std::chrono;
  7. typedef high_resolution_clock::time_point event;
  8.  
  9.  
  10. int main() {
  11.  
  12.     event start = high_resolution_clock::now();
  13.     sleep(1);
  14.     event stop = high_resolution_clock::now();
  15.    
  16.     duration<double> time_span = duration_cast<duration<double>>(stop - start);
  17.     auto elapsed_ns = std::chrono::duration_cast<std::chrono::nanoseconds>(stop - start);
  18.     auto elapsed_ms = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
  19.  
  20.     cout << "time " << time_span.count() << " seconds." << endl;
  21.     cout << "time " << elapsed_ns.count() << " ns." << endl;
  22.     cout << "time " << elapsed_ms.count() << " ms." << endl;
  23.    
  24.     return 0;
  25. }
  26.  

Voici le code introductif du programme :

  1. // ==================================================================
  2. // Program: dot_product.cpp
  3. // Date: July 2020
  4. // Author: Jean-Michel Richer
  5. // Email: jean-michel.richer@univ-angers.fr
  6. // ==================================================================
  7. // Description:
  8. //  This program enables to test different implementations of the
  9. // dot product of two vectors of floats
  10. // ==================================================================
  11. #include <iostream>
  12. #include <string>
  13. #include <stdint.h>
  14. #include <getopt.h>
  15. #include <xmmintrin.h>
  16. #include "cpu_timer.h"
  17. using namespace std;
  18.  
  19. // ------------------------------------------------------------------
  20. //
  21. // TYPES
  22. //
  23. // ------------------------------------------------------------------
  24. typedef float f32;
  25. typedef uint32_t u32;
  26. typedef f32 (*DotProductFn)(f32 *x, f32 *y, u32 size);
  27.  
  28. // ------------------------------------------------------------------
  29. //
  30. // GLOBAL VARIABLES
  31. //
  32. // ------------------------------------------------------------------
  33.  
  34. // use 16 for SSE, 32 for AVX, 64 for AVX512
  35. const int ALIGNMENT = 32;
  36.  
  37. // method to test
  38. u32 method = 1;
  39. // number of times the method is called
  40. u32 zillions = 1000;
  41. // size of each vector
  42. u32 vector_size = 64;
  43. // the two vectors used for the dot product
  44. f32 *vector_x, *vector_y;
  45.  
  46.  
  47.  
  48. // ------------------------------------------------------------------
  49. //
  50. // IMPLEMENTATION OF METHODS
  51. //
  52. // ------------------------------------------------------------------
  53.  
  54.  
  55. /**
  56.  * Method that computes the dot product of two vectors
  57.  * @param x pointer to first vector
  58.  * @param y pointer to second vector
  59.  * @param size size of the vectors
  60.  * @return the dot product of x and y, i.e. x_0 * y_0 + ... + x_(n-1) * y_(n-1)
  61.  */
  62. f32 dp_ref(f32 *x, f32 *y, u32 size) {
  63.     f32 sum = 0.0;
  64.    
  65.     for (u32 i = 0; i < size; ++i) {
  66.         sum += x[i] * y[i];
  67.     }
  68.    
  69.     return sum;
  70. }
  71.  
  72. // ------------------------------------------------------------------
  73. //
  74. // METHODS
  75. //
  76. // ------------------------------------------------------------------
  77.  
  78. // register assembly methods here
  79. extern "C" {
  80.  
  81. f32 dp_fpu(f32 *x, f32 *y, u32 size) ;
  82.  
  83. }
  84.  
  85. typedef struct {
  86.     DotProductFn function;
  87.     string name;
  88. } Method ;
  89.  
  90. #define entry(name) { name, #name }
  91.  
  92. // define methods here
  93.  
  94. Method methods [ ] = {
  95.     { nullptr, "undefined" },
  96.     entry(dp_ref),
  97.     entry(dp_fpu),
  98.     { nullptr, "undefined" }
  99. }; 
  100.  
  101.  
  102. /**
  103.  * main subprogram
  104.  */
  105. int main(int argc, char *argv[]) {
  106.  
  107.     // get command line arguments
  108.     // use getopt
  109.    
  110.  
  111.     //
  112.     // create vectors
  113.     //
  114.    
  115.     vector_x = (f32 *) _mm_malloc( vector_size * sizeof(f32), ALIGNMENT);
  116.     vector_y = (f32 *) _mm_malloc( vector_size * sizeof(f32), ALIGNMENT);
  117.    
  118.     //
  119.     // initialize vectors
  120.     //
  121.     const u32 modulo = 17; 
  122.  
  123.     for (u32 i = 0; i < vector_size; ++i) {
  124.         vector_x[i] = ((i % 3) == 0) ? -1.0 : 1.0;
  125.         vector_y[i] = (f32) 1 + (i % modulo);
  126.     }
  127.  
  128.     cout << "method=" << method << endl;
  129.     cout << "method_name=" << methods[ method ].name << endl;
  130.    
  131.     //
  132.     // Measure of the number of cycles of the treatment
  133.     //
  134.     CPUTimer timer;
  135.    
  136.     timer.start();
  137.    
  138.     float total = 0;
  139.     float result = 0;
  140.     u32 z;
  141.    
  142.     result = methods[ method ].function( vector_x, vector_y, vector_size );
  143.     total += result;
  144.     cout << "result=" << std::fixed << result << endl;
  145.        
  146.     for (z = 2; z <= zillions; ++z) {
  147.         result = methods[ method ].function( vector_x, vector_y, vector_size );
  148.         total += result;
  149.     }
  150.    
  151.     timer.stop();
  152.    
  153.    
  154.     cout << "cycles=" << timer << endl;
  155.     cout << "total=" << total << endl;
  156.    
  157.     return EXIT_SUCCESS;
  158. }
  159.  
  160.  
  161.  

5.2. Méthode utilisant la FPU

Ecrire la méthode dp_fpu qui réalise le calcul du produit scalaire en utilisant la FPU et l'intégrer au programme. Pour cela on utilisera les instructions :

Afin de coder la méthode dp_ref du listing précédent, on utilisera les conventions suivantes :

5.3. Automatisation de la compilation

Mettre au point un makefile afin d'automatiser la compilation. On doit notamment prendre en compte lors de l'édition de liens les fichiers objets suivants :

Comme on va utiliser les instructions SSE dans leur version intrinsics plus tard, il est nécessaire d'utiliser l'option de compilation -msse4.2 de gcc/g++.

5.4. Méthode utilisant la partie basse des unités vectorielles SSE

Ecrire la méthode dp_sse_low qui réalise le calcul du produit scalaire en utilisant les 32 premiers bits du vecteur SSE xmm0 qui représentera la variable sum. Utiliser dans ce cas les instructions :

A la fin du calcul il faut placer la partie basse de xmm0 dans st0 puisque la convention d'appel du C en 32 bits impose que le résultat soit placé dans st0. Pour cela, on réserve de l'espace mémoire dans la pile en abaissant le sommet de pile de 4 octets. On oubliera pas, après le transfert de la donnée, de remonter le sommet de pile de 4 octets.

5.5. Méthode utilisant les unités vectorielles SSE

Ecrire la méthode dp_sse_vec qui réalise le calcul du produit scalaire en utilisant $4 × 32$ bits d'un vecteur SSE. Vous pouvez utiliser :

Le code équivalent en C doit ressembler à :

  1. /**
  2.  * Use of SSE vector units
  3.  */
  4. f32 dp_sse_vec(f32 *x, f32 *y, u32 size) {
  5.     // for example use xmm0 to store v_sum
  6.     f32 v_sum[4];
  7.     // i must be defined here because it will be used
  8.     // by the for loop and then the while loop
  9.     u32 i;
  10.    
  11.     // use xorps xmm0, xmm0
  12.     v_sum[0] = v_sum[1] = v_sum[2] = v_sum[3] = 0.0;
  13.    
  14.     // loop unrolling by a factor of 4
  15.     // use xmm1 to store x[i:i+3]
  16.     for (i = 0; i < (size & ~3); i += 4) {
  17.         v_sum[0] += x[i + 0] * y[i + 0];
  18.         v_sum[1] += x[i + 1] * y[i + 1];
  19.         v_sum[2] += x[i + 2] * y[i + 2];
  20.         v_sum[3] += x[i + 3] * y[i + 3];
  21.     }
  22.    
  23.     // sum of partial sums
  24.     float sum = v_sum[0] + v_sum[1] + v_sum[2] + v_sum[3];
  25.    
  26.     // last iterations
  27.     while (i < size) {
  28.         sum += x[i] * y[i];
  29.         ++i;
  30.     }
  31.    
  32.     return sum;
  33. }
  34.  

Les instructions SSE à utiliser sont les suivantes :

5.6. Méthode utilisant les unités vectorielles AVX

Ecrire la méthode dp_avx_vec qui réalise le calcul du produit scalaire en utilisant $8 × 32$ bits d'un vecteur AVX. Utiliser dans ce cas les instructions :

  1. /**
  2.  * Use of AVX vector units
  3.  */
  4. f32 dp_avx_vec(f32 *x, f32 *y, u32 size) {
  5.    // use ymm0 to store v_sum
  6.    f32 v_sum[8];
  7.     u32 i;
  8.    
  9.     v_sum[0] = v_sum[1] = v_sum[2] = v_sum[3] = 0.0;
  10.     v_sum[4] = v_sum[5] = v_sum[6] = v_sum[7] = 0.0;
  11.    
  12.     // loop unrolling by a factor of 8
  13.     // use ymm1 to store x[i:i+7]
  14.     for (i = 0; i < (size & ~7); i += 8) {
  15.         v_sum[0] += x[i + 0] * y[i + 0];
  16.         v_sum[1] += x[i + 1] * y[i + 1];
  17.         v_sum[2] += x[i + 2] * y[i + 2];
  18.         v_sum[3] += x[i + 3] * y[i + 3];
  19.         v_sum[4] += x[i + 4] * y[i + 4];
  20.         v_sum[5] += x[i + 5] * y[i + 5];
  21.         v_sum[6] += x[i + 6] * y[i + 6];
  22.         v_sum[7] += x[i + 7] * y[i + 7];
  23.     }
  24.    
  25.     // sum of partial sums
  26.     float sum = v_sum[0] + v_sum[1] + v_sum[2] + v_sum[3] +;
  27.         v_sum[4] + v_sum[5] + v_sum[6] + v_sum[7];
  28.        
  29.     // last iterations 
  30.     while (i < size) {
  31.         sum += x[i] * y[i];
  32.         ++i;
  33.     }
  34.    
  35.     return sum;
  36. }
  37.  

5.7. Méthode utilisant les intrinsics

Ecrire une version intrinsics dp_sse_intrin qui sera traduite par le compilateur C++. On utilisera le site Intel Intrinsics Guide afin d'obtenir une correspondance entre instructions assembleur et méthodes intrinsics.