G++ SSE intrinsics dilemma-значение от intrinsic " насыщает"


Я написал простую программу для реализации SSE intrinsics для вычисления внутреннего произведения двух больших (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++ для построения это:

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

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

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

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

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

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

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

Спасибо,
Шрирам.

1   5  

1 ответ:

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

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

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