SSE Intrinsics арифметическая ошибка

Question

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

Для тестирования я установил каждый элемент обоих массивов в 1, поэтому продукт должен быть размером ==.

Он работает правильно, но всякий раз, когда я запускаю код размером ~ ~ 68000000, код, использующий sse intrinsics, начинает вычислять неправильный внутренний продукт. Кажется, что он застрял на определенной сумме и никогда не превышает этого числа. Вот пример запуска:

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 как максимальное вычисленное значение.