AVX2指令集優化浮點陣列求和演演算法

2022-05-18 16:00:57

一、AVX2指令集介紹

AVX2是SIMD(單指令多資料流)指令集，支援在一個指令週期內同時對256位記憶體進行操作。包含乘法，加法，位運算等功能。下附Intel官網使用檔案。

我們本次要用到的指令有 __m256i _mm256_add_pd(__m256i a, __m256i b), __m256i _mm256_add_ps等，（p代表精度precision，s代表single，d代表double）

它們可以一次取256位的記憶體，並按32/64位元一個浮點進行加法運算。下附官網描述。

Synopsis

__m256d _mm256_add_pd (__m256d a, __m256d b)

#include <immintrin.h>

Instruction: vaddpd ymm, ymm, ymm

CPUID Flags: AVX

Description

Add packed double-precision (64-bit) floating-point elements in a and b, and store the results in dst.

Operation

FOR j := 0 to 3
	i := j*64
	dst[i+63:i] := a[i+63:i] + b[i+63:i]
ENDFOR
dst[MAX:256] := 0

Performance

Architecture	Latency	Throughput (CPI)
Icelake	4	0.5
Skylake	4	0.5
Broadwell	3	1
Haswell	3	1
Ivy Bridge	3	1

二、程式碼實現

0. 資料生成

為了比較結果，我們生成從1到N的等差數列。這裡利用模版相容不同資料型別。由於AVX2指令集一次要操作多個資料，為了防止訪存越界，我們將大小擴充套件到256的整數倍位位元，也就是32位元組的整數倍。

uint64_t lowbit(uint64_t x)
{
    return x & (-x);
}
uint64_t extTo2Power(uint64_t n, int i)//arraysize datasize
{
    while(lowbit(n) < i)
        n += lowbit(n);
    return n;
}

template <typename T>
T* getArray(uint64_t size)
{
    uint64_t ExSize = extTo2Power(size, 32/sizeof(T));
    T* arr = new T[ExSize];
    for (uint64_t i = 0; i < size; i++)
        arr[i] = i+1;
    for (uint64_t i = size; i < ExSize; i++)
        arr[i] = 0;
    return arr;
}

1. 普通陣列求和

為了比較效能差異，我們先實現一份普通的陣列求和。這裡也使用模版。

template <typename T>
T simpleSum(T* arr, uint64_t size)
{
    T sum = 0;
    for (uint64_t i = 0; i < size; i++)
        sum += arr[i];
    return sum;
}

2. AVX2指令集求和：單精度浮點(float)

這裡我們預開一個avx2的整形變數，每次從陣列中取8個32位元浮點，加到這個變數上，最後在對這8個32位元浮點求和。

float avx2Sum(float* arr, uint64_t size)
{
    float sum[8] = {0};
    __m256 sum256 = _mm256_setzero_ps();
    __m256 load256 = _mm256_setzero_ps();
    for (uint64_t i = 0; i < size; i += 8)
    {
        load256 = _mm256_loadu_ps(&arr[i]);
        sum256 = _mm256_add_ps(sum256, load256);
    }
    sum256 = _mm256_hadd_ps(sum256, sum256);
    sum256 = _mm256_hadd_ps(sum256, sum256);
    _mm256_storeu_ps(sum, sum256);
    sum[0] += sum[4];
    return sum[0];
}

這裡的hadd是橫向加法，具體實現類似下圖，可以幫我們實現陣列內求和：

3. AVX2指令集求和：雙精度浮點(double)

double avx2Sum(double* arr, uint64_t size)
{
    double sum[4] = {0};
    __m256d sum256 = _mm256_setzero_pd();
    __m256d load256 = _mm256_setzero_pd();
    for (uint64_t i = 0; i < size; i += 4)
    {
        load256 = _mm256_loadu_pd(&arr[i]);
        sum256 = _mm256_add_pd(sum256, load256);
    }
    sum256 = _mm256_hadd_pd(sum256, sum256);
    _mm256_storeu_pd(sum, sum256);
    sum[0] += sum[2];
    return sum[0];
}

三、效能測試

測試環境

Device	Description
CPU	Intel Core i9-9880H 8-core 2.3GHz
Memory	DDR4-2400MHz Dual-Channel 32GB
complier	Apple Clang-1300.0.29.30

計時方式

利用chrono庫獲取系統時鐘計算執行時間，精確到毫秒級

uint64_t getTime()
{
    uint64_t timems = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
    return timems;
}

測試內容

對1到1e9求和，答案應該為500000000500000000, 分別測試float和double。

	uint64_t N = 1e9;
    // compare the performance of normal add and avx2 add
    uint64_t start, end;
    // test float
    cout << "compare float sum: " << endl;
    float* arr3 = getArray<float>(N);
    start = getTime();
    float sum3 = simpleSum(arr3, N);
    end = getTime();
    cout << "float simpleSum time: " << end - start << endl;
    cout << "float simpleSum sum: " << sum3 << endl;
    start = getTime();
    sum3 = avx2Sum(arr3, N);
    end = getTime();
    cout << "float avx2Sum time: " << end - start << endl;
    cout << "float avx2Sum sum: " << sum3 << endl;
    delete[] arr3;
    cout << endl << endl;
    // test double
    cout << "compare double sum: " << endl;
    double* arr4 = getArray<double>(N);
    start = getTime();
    double sum4 = simpleSum(arr4, N);
    end = getTime();
    cout << "double simpleSum time: " << end - start << endl;
    cout << "double simpleSum sum: " << sum4 << endl;
    start = getTime();
    sum4 = avx2Sum(arr4, N);
    end = getTime();
    cout << "double avx2Sum time: " << end - start << endl;
    cout << "double avx2Sum sum: " << sum4 << endl;
    delete[] arr4;
    cout << endl << endl;

進行效能測試

第一次測試

測試命令

g++ -mavx2 avx_big_integer.cpp 
./a.out

測試結果

方法	耗時(ms)
AVX2加法單精度	615
普通加法單精度	2229
AVX2加法雙精度	1237
普通加法雙精度	2426

這裡能看到單精度下已經出現了比較明顯的誤差，並且由於普通求和和avx2求和的加法順序不一樣，導致誤差值也不一樣。

第二次測試