Ce site est en cours de reconstruction certains liens peuvent ne pas fonctionner ou certaines images peuvent ne pas s'afficher.

5. TP 5
Produit scalaire

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 :

la FPU
la partie basse des registres vectoriels SSE
la vectorisation avec SSE
la vectorisation avec AVX

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 :

dot_product.cpp comprendra la méthode main, gérera les arguments en ligne de commande avec getopt, etc...
dot_product_nasm.asm comprendra le code assembleur

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

Variables et paramètres en ligne de commande
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

le nombre de répétitions est le nombre de fois où le produit scalaire est exécuté
la taille des vecteurs concerne $x$ et $y$ (variables vector_x et vector_y)
la méthode de calcul utilisée est identifiée par son numéro de 1 à $n$

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 :

Warning: file_get_contents(ens/l3/archi/time_it.cpp): Failed to open stream: No such file or directory in /home/jeanmichel.richer/public_html/ez_web.php on line 418

Afficher le code ens/l3/archi/time_it.cpp

Voici le code introductif du programme :

Afficher le code assembly/programs/dot_product.cpp

// ==================================================================
// Program: dot_product.cpp
// Date: July 2020
// Author: Jean-Michel Richer
// Email: jean-michel.richer@univ-angers.fr
// ==================================================================
// Description:
//  This program enables to test different implementations of the
// dot product of two vectors of floats
// ==================================================================
#include <iostream>
#include <string>
#include <stdint.h>
#include <getopt.h>
#include <xmmintrin.h>
#include "cpu_timer.h"
using namespace std;
 
// ------------------------------------------------------------------
//
// TYPES
//
// ------------------------------------------------------------------
typedef float f32;
typedef uint32_t u32;
typedef f32 (*DotProductFn)(f32 *x, f32 *y, u32 size);
 
// ------------------------------------------------------------------
//
// GLOBAL VARIABLES
//
// ------------------------------------------------------------------
 
// use 16 for SSE, 32 for AVX, 64 for AVX512
const int ALIGNMENT = 32;
 
// method to test
u32 method = 1;
// number of times the method is called
u32 zillions = 1000;
// size of each vector
u32 vector_size = 64;
// the two vectors used for the dot product
f32 *vector_x, *vector_y;
 
 
 
// ------------------------------------------------------------------
//
// IMPLEMENTATION OF METHODS
//
// ------------------------------------------------------------------
 
 
/**
 * Method that computes the dot product of two vectors
 * @param x pointer to first vector
 * @param y pointer to second vector
 * @param size size of the vectors
 * @return the dot product of x and y, i.e. x_0 * y_0 + ... + x_(n-1) * y_(n-1)
 */
f32 dp_ref(f32 *x, f32 *y, u32 size) {
  f32 sum = 0.0;
  
  for (u32 i = 0; i < size; ++i) {
    sum += x[i] * y[i];
  }
  
  return sum;
}
 
// ------------------------------------------------------------------
//
// METHODS
//
// ------------------------------------------------------------------
 
// register assembly methods here
extern "C" {
 
f32 dp_fpu(f32 *x, f32 *y, u32 size) ;
 
}
 
typedef struct {
  DotProductFn function;
  string name;
} Method ;
 
#define entry(name) { name, #name }
 
// define methods here
 
Method methods [ ] = {
  { nullptr, "undefined" },
  entry(dp_ref),
  entry(dp_fpu),
  { nullptr, "undefined" }
};  
 
 
/**
 * main subprogram
 */
int main(int argc, char *argv[]) {
 
  // get command line arguments
  // use getopt
  
 
  //
  // create vectors
  //
  
  vector_x = (f32 *) _mm_malloc( vector_size * sizeof(f32), ALIGNMENT);
  vector_y = (f32 *) _mm_malloc( vector_size * sizeof(f32), ALIGNMENT);
  
  //
  // initialize vectors
  //
  const u32 modulo = 17;  
 
  for (u32 i = 0; i < vector_size; ++i) {
    vector_x[i] = ((i % 3) == 0) ? -1.0 : 1.0;
    vector_y[i] = (f32) 1 + (i % modulo);
  }
 
  cout << "method=" << method << endl;
  cout << "method_name=" << methods[ method ].name << endl;
  
  //
  // Measure of the number of cycles of the treatment
  // 
  CPUTimer timer;
  
  timer.start();
  
  float total = 0;
  float result = 0;
  u32 z;
  
  result = methods[ method ].function( vector_x, vector_y, vector_size );
  total += result;
  cout << "result=" << std::fixed << result << endl;
    
  for (z = 2; z <= zillions; ++z) {
    result = methods[ method ].function( vector_x, vector_y, vector_size );
    total += result;
  }
  
  timer.stop();
  
  
  cout << "cycles=" << timer << endl;
  cout << "total=" << total << endl;
  
  return EXIT_SUCCESS;
}
 
 
 

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 :

fldz pour mettre la variable sum à 0.
fld pour charger les données
fmul / fmulp pour réaliser la multiplication
fadd / faddp pour réaliser l'addition
fst /fstp pour stocker les données

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

la variable sum sera remplacée par st0
le registre eax stockera l'adresse du vecteur x
le registre ebx stockera l'adresse du vecteur y
la variable de boucle sera stockée dans ecx
le registre edx stockera size

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 :

dot_product.o
dot_product_nasm.o
cpu_timer.o

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 :

xorps afin de mettre à zéro le vecteur xmm0
movss afin de charger x[i] dans la partie basse de xmm1
mulss pour multiplier xmm1 par y[i]
addss pour additionner xmm1 à xmm0

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 :

xmm0 pour représenter le tableau sum_a
xmm1 pour stocker x[i:i+3] et multiplier xmm1 par y[i:i+3]

Le code équivalent en C doit ressembler à :

Afficher le code assembly/programs/dp_sse_vec.cpp

/**
 * Use of SSE vector units
 */
f32 dp_sse_vec(f32 *x, f32 *y, u32 size) {
  // for example use xmm0 to store v_sum
    f32 v_sum[4];
    // i must be defined here because it will be used
    // by the for loop and then the while loop
    u32 i;
    
    // use xorps xmm0, xmm0
    v_sum[0] = v_sum[1] = v_sum[2] = v_sum[3] = 0.0;
   
    // loop unrolling by a factor of 4
    // use xmm1 to store x[i:i+3]
    for (i = 0; i < (size & ~3); i += 4) {
        v_sum[0] += x[i + 0] * y[i + 0];
        v_sum[1] += x[i + 1] * y[i + 1];
        v_sum[2] += x[i + 2] * y[i + 2];
        v_sum[3] += x[i + 3] * y[i + 3];
    }
    
    // sum of partial sums
    float sum = v_sum[0] + v_sum[1] + v_sum[2] + v_sum[3];
    
    // last iterations
    while (i < size) {
      sum += x[i] * y[i];
      ++i;
    }
   
    return sum;
}
 

Les instructions SSE à utiliser sont les suivantes :

xorps afin de mettre à zéro le vecteur xmm0
movaps pour charger x[i:i+3] dans xmm1
mulps afin de multiplier xmm1 par y[i:i+3]
addps pour additionner xmm1 à xmm0
haddps pour calculer la somme des valeurs de xmm0

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 :

Afficher le code assembly/programs/dp_avx_vec.cpp

/**
 * Use of AVX vector units 
 */
f32 dp_avx_vec(f32 *x, f32 *y, u32 size) {
   // use ymm0 to store v_sum
   f32 v_sum[8];
    u32 i;
    
    v_sum[0] = v_sum[1] = v_sum[2] = v_sum[3] = 0.0;
    v_sum[4] = v_sum[5] = v_sum[6] = v_sum[7] = 0.0;
   
    // loop unrolling by a factor of 8
    // use ymm1 to store x[i:i+7]
    for (i = 0; i < (size & ~7); i += 8) {
        v_sum[0] += x[i + 0] * y[i + 0];
        v_sum[1] += x[i + 1] * y[i + 1];
        v_sum[2] += x[i + 2] * y[i + 2];
        v_sum[3] += x[i + 3] * y[i + 3];
        v_sum[4] += x[i + 4] * y[i + 4];
        v_sum[5] += x[i + 5] * y[i + 5];
        v_sum[6] += x[i + 6] * y[i + 6];
        v_sum[7] += x[i + 7] * y[i + 7];
    }
    
    // sum of partial sums
    float sum = v_sum[0] + v_sum[1] + v_sum[2] + v_sum[3] +;
      v_sum[4] + v_sum[5] + v_sum[6] + v_sum[7];
      
    // last iterations  
    while (i < size) {
      sum += x[i] * y[i];
      ++i;
    }
   
    return sum;
}
 

vxorps afin de mettre à zéro le vecteur ymm0
vmovaps
vmulps
vaddps
vhaddps (voir le cours Section 8.3.3 - Instructions singulières)

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.

__m128 _mm_setzero_ps() fixe le registre à 0
__m128 _mm_load_ps (float const* mem_addr) permet de charger les données
__m128 _mm_add_ps (__m128 a, __m128 b) réalise l'addition de deux vecteurs
__m128 _mm_mul_ps (__m128 a, __m128 b) réalise la multiplication de deux vecteurs
__m128 _mm_hadd_ps(__m128 a, __m128 b) réalise la somme dite horizontale
void _mm_store_ss(float* mem_addr, __m128 a) stocke la partie basse d'un vecteur dans un espace mémoire de type float en indiquant l'adresse de cet espace

5. TP 5Produit scalaire