SSE整数算法误差

Question

我一直在尝试SSE的本质，我似乎遇到了一个奇怪的错误，我无法弄清楚。我正在计算两个浮点数组的内积，一次计算4个元素。

为了进行测试，我将两个数组的每个元素设置为1，因此产品应该是==大小。

它正确运行，但是每当我运行大小大于68000000的代码时，使用sse本质的代码就会开始计算错误的内部积。它似乎被卡在一定的金额上，而且永远不会超过这个数字。下面是一个运行示例：

joe:~$./test_sse 70000000
sequential inner product: 70000000.000000
sse        inner product: 67108864.000000
sequential          time: 0.417932
sse                 time: 0.274255

汇编：

gcc -fopenmp test_sse.c -o test_sse -std=c99

这个错误在我测试过的几台计算机中似乎是一致的。下面是代码，也许有人能帮我弄清楚到底发生了什么：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>
#include <math.h>
#include <assert.h>

#include <xmmintrin.h>

double inner_product_sequential(float * a, float * b, unsigned int size) {

  double sum = 0;

  for(unsigned int i = 0; i < size; i++) {
    sum += a[i] * b[i];
  }

  return sum;

}

double inner_product_sse(float * a, float * b, unsigned int size) {

  assert(size % 4 == 0);

  __m128 X, Y, Z;

  Z = _mm_set1_ps(0.0f);

  float arr[4] __attribute__((aligned(sizeof(float) * 4)));

  for(unsigned int i = 0; i < size; i += 4) {

    X = _mm_load_ps(a+i);
    Y = _mm_load_ps(b+i);

    X = _mm_mul_ps(X, Y);
    Z = _mm_add_ps(X, Z);

  }

  _mm_store_ps(arr, Z);

  return arr[0] + arr[1] + arr[2] + arr[3];

}

int main(int argc, char ** argv) {

  if(argc < 2) {
    fprintf(stderr, "usage: ./test_sse <size>\n");
    exit(EXIT_FAILURE);
  }

  unsigned int size = atoi(argv[1]);

  srand(time(0));

  float *a = (float *) _mm_malloc(size * sizeof(float), sizeof(float) * 4);
  float *b = (float *) _mm_malloc(size * sizeof(float), sizeof(float) * 4);

  for(int i = 0; i < size; i++) {
    a[i] = b[i] = 1;
  }



  double start, time_seq, time_sse;


  start = omp_get_wtime();

  double inner_seq = inner_product_sequential(a, b, size);

  time_seq = omp_get_wtime() - start;


  start = omp_get_wtime();

  double inner_sse = inner_product_sse(a, b, size);

  time_sse = omp_get_wtime() - start;


  printf("sequential inner product: %f\n", inner_seq);
  printf("sse        inner product: %f\n", inner_sse);
  printf("sequential          time: %f\n", time_seq);
  printf("sse                 time: %f\n", time_sse);




  _mm_free(a);
  _mm_free(b);
}

Answer 1

您正遇到单个精确浮点数的精度限制。数字16777216 (2^ 24 )是向量Z达到“极限”内积时各分量的值，用32位浮点表示为十六进制0x4b800000或二进制0 10010111 00000000000000000000000，即23位尾数全部为零(隐式前导1位)，8位指数部分为151表示指数151-127=24。如果将1加到该值中，这将需要增加指数，但不能再在尾数中表示所加的指数，因此在单精度浮点算术2^24 +1= 2^24。

您没有在顺序函数中看到这一点，因为您使用64位双精度值来存储结果，而且当我们在x86平台上工作时，内部最可能使用的是80位超额精度寄存器。

您可以在顺序代码中强制使用单个精度，方法是将其重写为

float sum;

float inner_product_sequential(float * a, float * b, unsigned int size) {
  sum = 0;

  for(unsigned int i = 0; i < size; i++) {
    sum += a[i] * b[i];
  }

  return sum;
}

您将看到16777216.000000是计算值的最大值。