g ++ SSE intrínsecos dilema - valor de intrínsecas "satura"

Question

I escribió un programa simple de implementar características intrínsecas de ESS para calcular el producto interno de dos grandes (100000 o más elementos) vectores. El programa compara el tiempo de ejecución para ambos, el producto interno calculado de la manera convencional y el uso de intrínsecos. Todo funciona bien, hasta que inserte (solo por diversión) un bucle interno antes de la declaración que calcula el producto interno. Antes de seguir adelante, aquí está el código:

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

utilizo el siguiente g invocación ++ para generar esto:

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

Cada uno de los bucles anteriormente, en tanto las funciones, se ejecuta un total de N^2 veces. Sin embargo, dado que arr1 y arr2 (los dos vectores de coma flotante) se cargan con un valor de 2.5, la longitud de la matriz es 100.000, el resultado en ambos casos debería ser 6.25e + 10. Los resultados que obtengo son:

Resultados:
result_Conventional = 6.25e + 10 =
result_Intrinsics 5.36871e + 08

Esto no es todo. Parece que el valor devuelto por la función que usa intrínsecos "satura" en el valor anterior. Traté de poner otros valores para los elementos de la matriz y diferentes tamaños también. Pero parece que cualquier valor por encima de 1.0 para los contenidos de la matriz y cualquier tamaño superior a 1000 cumple con el mismo valor que se ve arriba.

Inicialmente, pensé que podría ser debido a que todas las operaciones dentro SSE están en coma flotante, pero coma flotante deben ser capaces de almacenar un número que es del orden de e + 08.

estoy tratando de ver donde podía ir mal, pero parece que no puede entenderlo. Estoy usando la versión de g ++: g ++ (GCC) 4.4.1 20090725 (Red Hat 4.4.1-2).

Cualquier ayuda al respecto es bienvenida.

Gracias,
Sriram.

Answer 1

El problema que está teniendo es que mientras que un float puede almacenar 6.25e + 10, solo tiene unos pocos dígitos significativos de precisión.

Esto significa que cuando construye un número grande sumando muchos números pequeños juntos de a poco, llega a un punto donde el número más pequeño es menor que el dígito de precisión más bajo en el número más grande, por lo que lo suma no tiene efecto.

En cuanto a por qué no obtiene este comportamiento en la versión no intrínseca, es probable que result variable se mantenga en un registro que utiliza una mayor precisión que el almacenamiento real de un flotante para que no se trunque a la precisión de float en cada iteración del ciclo. Tendría que mirar el código de ensamblador generado para estar seguro.