c - SSE Intrinsics 算术错误

Question

我一直在试验 SSE 内在函数，我似乎遇到了一个我无法弄清楚的奇怪错误。我正在计算两个 float 组的内积，一次计算 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

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

您在顺序函数中看不到这一点，因为您使用 64 位 double 值来存储结果，并且由于我们在 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 作为最大计算值。