g ++ Внутренняя дилемма SSE - значение от встроенных «насыщенных»

Question

Я написал простую программу для реализации SSE-функций для вычисления скалярного произведения двух больших (100000 и более элементов) векторов. Программа сравнивает время выполнения для обоих, внутренний продукт вычисляется обычным способом и использует внутреннюю среду. Все работает отлично, пока я не вставляю (просто для удовольствия) внутренний цикл перед оператором, который вычисляет внутренний продукт. Прежде чем идти дальше, вот код:

//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; 
     } 

Я использую следующий г призывание ++ построить так:

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

Каждая из петель выше, в обеих функциях, работает в общей сложности N^2 раза. Однако, учитывая, что arr1 и arr2 (два вектора с плавающей запятой) загружаются со значением 2,5, длина массива равна 100 000, результат в обоих случаях должен быть 6.25e + 10. Результаты я получаю:

Результаты:
result_Conventional = 6.25e + 10
result_Intrinsics = 5.36871e + 08

Это не все. Кажется, что значение, возвращаемое функцией, которая использует intrinsics, «насыщается» по значению выше. Я попытался поставить другие значения для элементов массива и разных размеров. Но кажется, что любое значение выше 1.0 для содержимого массива и любого размера выше 1000 соответствует тому же значению, которое мы видим выше.

Первоначально я думал, что это может быть потому, что все операции в SSE находятся в плавающей запятой, но с плавающей запятой следует хранить число, имеющее порядок e + 08.

Я пытаюсь понять, где я могу ошибиться, но, похоже, не могу понять. Я использую версию g ++: g ++ (GCC) 4.4.1 20090725 (Red Hat 4.4.1-2).

Любая помощь в этом приветствуется.

Thanks,
Sriram.

Answer 1

Проблема, с которой вы сталкиваетесь, заключается в том, что в то время как float может хранить 6.25e + 10, он имеет только несколько значащих цифр точности.

Это означает, что, когда вы строите большое количество, добавляя множество небольших чисел друг к другу понемногу, вы достигаете точки, где меньшее число меньше, чем самая низкая цифра точности в большем количестве, поэтому добавляя ее вверх не имеет никакого эффекта.

Относительно того, почему вы не получаете это поведение в неинтерминированной версии, вероятно, что переменная result хранится в регистре, которая использует более высокую точность, чем фактическое хранение поплавка, поэтому оно не усекается с точностью до float на каждой итерации цикла. Вам нужно будет посмотреть на сгенерированный код ассемблера.