g++上证所内在困境——内在“饱和”带来的价值

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++ 调用来构建此代码：

 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。我得到的结果是：

结果：
结果_常规 = 6.25e+10
result_Intrinsics = 5.36871e+08

这还不是全部。似乎从使用内在函数的函数返回的值在上述值处“饱和”。我尝试为数组的元素设置其他值以及不同的大小。但似乎数组内容大于 1.0 的任何值以及大于 1000 的任何大小都符合我们上面看到的相同值。

最初，我认为这可能是因为SSE中的所有操作都是浮点型的，但是浮点型应该能够存储e+08量级的数字。

我试图找出哪里可能出错，但似乎无法弄清楚。我使用的 g++ 版本：g++ (GCC) 4.4.1 20090725 (Red Hat 4.4.1-2)。

非常欢迎对此提供任何帮助。

谢谢，
斯里拉姆。

Answer 1

您遇到的问题是，虽然 float 可以存储 6.25e+10，但它只有几个有效数字的精度。

这意味着，当您通过一次一点地添加大量小数字来构建一个大数字时，您会达到一个点，即较小的数字小于较大数字中的最低精度数字，因此将其相加没有效果。

至于为什么在非内在版本中没有得到这种行为，很可能 result 变量被保存在一个寄存器中，该寄存器使用比浮点的实际存储更高的精度，因此它是在循环的每次迭代中不会被截断为 float 的精度。您必须查看生成的汇编代码才能确定。