g ++ dilemme SSE - valeur intrinsèque de « sature »

Question

J'ai écrit un programme simple à mettre en œuvre intrinsics ESS pour le calcul du produit intérieur de deux grands (100000 éléments ou plus) vecteurs. Le programme compare le temps d'exécution pour les deux, le produit scalaire calculé de manière classique et en utilisant des valeurs intrinsèques. Tout fonctionne bien, jusqu'à ce que j'insérer (juste pour le plaisir de celui-ci) une boucle interne avant l'instruction qui calcule le produit intérieur. Avant d'aller plus loin, voici le code:

    //this is a sample Intrinsics program to compute inner product of two vectors and compare Intrinsics with traditional method of doing things.

        #include <iostream>
        #include <iomanip>
        #include <xmmintrin.h>
        #include <stdio.h>
        #include <time.h>
        #include <stdlib.h>
        using namespace std;

        typedef float v4sf __attribute__ ((vector_size(16)));

        double innerProduct(float* arr1, int len1, float* arr2, int len2) {  //assume len1 = len2.

          float result = 0.0;
          for(int i = 0; i < len1; i++) {
            for(int j = 0; j < len1; j++) {
              result += (arr1[i] * arr2[i]);
            }
          }

         //float y = 1.23e+09;
         //cout << "y = " << y << endl;
         return result;
        }

        double sse_v4sf_innerProduct(float* arr1, int len1, float* arr2, int len2) { //assume that len1 = len2.

          if(len1 != len2) {
            cout << "Lengths not equal." << endl;
            exit(1);
          }

          /*steps:
         * 1. load a long-type (4 float) into a v4sf type data from both arrays.
         * 2. multiply the two.
         * 3. multiply the same and store result.
         * 4. add this to previous results.
         */

          v4sf arr1Data, arr2Data, prevSums, multVal, xyz;
          //__builtin_ia32_xorps(prevSums, prevSums);   //making it equal zero.
         //can explicitly load 0 into prevSums using loadps or storeps (Check).

          float temp[4] = {0.0, 0.0, 0.0, 0.0};
          prevSums = __builtin_ia32_loadups(temp);
          float result = 0.0;

          for(int i = 0; i < (len1 - 3); i += 4) {
            for(int j = 0; j < len1; j++) {
            arr1Data = __builtin_ia32_loadups(&arr1[i]);
            arr2Data = __builtin_ia32_loadups(&arr2[i]);  //store the contents of two arrays.
            multVal = __builtin_ia32_mulps(arr1Data, arr2Data);   //multiply.
            xyz = __builtin_ia32_addps(multVal, prevSums);
            prevSums = xyz;
          }
         }
          //prevSums will hold the sums of 4 32-bit floating point values taken at a time. Individual entries in prevSums also need to be added.
          __builtin_ia32_storeups(temp, prevSums);  //store prevSums into temp.

           cout << "Values of temp:" << endl;
           for(int i = 0; i < 4; i++)
             cout << temp[i] << endl;

          result += temp[0] + temp[1] + temp[2] + temp[3];

        return result;
        }

        int main() {
          clock_t begin, end;
          int length = 100000;
          float *arr1, *arr2;
          double result_Conventional, result_Intrinsic;

 //         printStats("Allocating memory.");
          arr1 = new float[length];
          arr2 = new float[length];
 //         printStats("End allocation.");

          srand(time(NULL));  //init random seed.
 //         printStats("Initializing array1 and array2");
          begin = clock();
          for(int i = 0; i < length; i++) {
         //   for(int j = 0; j < length; j++) {
          //    arr1[i] = rand() % 10 + 1;
                arr1[i] = 2.5;
           //    arr2[i] = rand() % 10 - 1;
                arr2[i] = 2.5;
         //   }
          }
          end = clock();
          cout << "Time to initialize array1 and array2 = " << ((double) (end - begin)) / CLOCKS_PER_SEC << endl;
  //        printStats("Finished initialization.");

    //      printStats("Begin inner product conventionally.");
          begin = clock();
          result_Conventional = innerProduct(arr1, length, arr2, length);
          end = clock();
          cout << "Time to compute inner product conventionally = " << ((double) (end - begin)) / CLOCKS_PER_SEC << endl;
    //      printStats("End inner product conventionally.");

      //    printStats("Begin inner product using Intrinsics.");
          begin = clock();
          result_Intrinsic = sse_v4sf_innerProduct(arr1, length, arr2, length);
          end = clock();
          cout << "Time to compute inner product with intrinsics = " << ((double) (end - begin)) / CLOCKS_PER_SEC << endl;
          //printStats("End inner product using Intrinsics.");

          cout << "Results: " << endl;
          cout << " result_Conventional = " << result_Conventional << endl;
          cout << " result_Intrinsics = " << result_Intrinsic << endl;
        return 0;
        }

J'utilise les éléments suivants g ++ invocation de construire ceci:

 g++ -W -Wall -O2 -pedantic -march=i386 -msse intrinsics_SSE_innerProduct.C -o innerProduct  

Chacune des boucles ci-dessus, dans les deux fonctions, gère un total de N ^ 2 fois. Cependant, étant donné que arr1 et arr2 (les deux vecteurs à virgule flottante) sont chargés avec une valeur de 2,5, la longueur de la matrice est 100 000, le résultat dans les deux cas doit être 6.25e + 10. Les résultats que je reçois sont:

Résultats:
 result_Conventional = 6.25e + 10
 result_Intrinsics = 5.36871e + 08

Ce n'est pas tout. Il semble que la valeur retournée par la fonction qui utilise « sature » intrinsèques à la valeur ci-dessus. J'ai essayé de mettre d'autres valeurs pour les éléments du tableau et de tailles différentes aussi. Mais il semble que toute valeur supérieure à 1,0 pour le contenu du tableau et une taille supérieure à 1000 rencontre la même valeur que nous voyons ci-dessus.

Dans un premier temps, je pensais que ce serait peut-être parce que toutes les opérations au sein de l'ESS sont en virgule flottante, mais virgule flottante devrait être en mesure de stocker un nombre qui est de l'ordre de + 08 e.

Je suis en train de voir où je pouvais aller mal, mais ne peut pas sembler comprendre. J'utilise g ++ version: g ++ (GCC) 4.4.1 20.090.725 (Red Hat 4.4.1-2).

Toute aide à ce sujet est la bienvenue.

Merci,
Sriram.

Answer 1

Le problème que vous rencontrez est que si un float peut stocker 6.25e + 10, il n'a que quelques chiffres significatifs de précision.

Cela signifie que lorsque vous construisez un grand nombre en ajoutant beaucoup de petits nombres ensemble un peu à la fois, vous atteignez un point où le plus petit nombre est inférieur au chiffre de précision le plus bas dans le plus grand nombre afin de l'ajouter en a aucun effet.

Quant à savoir pourquoi vous ne recevez pas ce comportement dans la version non-intrinsèque, il est probable que variable result se tient dans un registre qui utilise une plus grande précision que le stockage réelle d'un flotteur de sorte qu'il ne soit pas tronqué la précision d'un float à chaque itération de la boucle. Vous devrez regarder le code assembleur généré pour être sûr.