SIMD Programming(1. Basic)

SIMD

就是所谓一条指令多个数据（Single Instruction Multiple Data）。

上个图，比如以下操作就可以一条指令完成。

Intel Intrisic

我们使用Intel Intrisic来简化SIMD指令的编写。

不然在不能使用OpenMP的MSVC的编译器上只能手写汇编了。

float a[4] = { 1,2,3,4 };
float b[4] = { 5,6,7,8 };
float res[4];

__asm 
{
    movups xmm0, [a];     // 将a所指内存的128位数据放入xmm0寄存器
    movups xmm1, [b];     // 将b所指内存的128位数据放入xmm0寄存器
    mulps xmm0, xmm1;     // 计算 xmm0 * xmm1，结果放入 xmm0
    movups [res], xmm0;   // 将xmm0的数据放入res所指内存
}

对于Intel Intrisic而言，每个头文件代表不同的指令集，你可以使用CPU-Z查看你自己的CPU支持什么指令集。然后在编译器选项上选择类似于march=native类似的选项即可。

#include <mmintrin.h>     // MMX
#include <xmmintrin.h>    // SSE
#include <emmintrin.h>    // SSE2
#include <pmmintrin.h>    // SSE3
#include <tmmintrin.h>    // SSSE3
#include <smmintrin.h>    // SSE4.1
#include <nmmintrin.h>    // SSE4.2
#include <wmmintrin.h>    // AES
#include <immintrin.h>    // AVX
#include <intrin.h>

我们使用Intel Intrisic以后，之前的汇编代码就可以写成

#include <xmmintrin.h>

float a[4] = { 1,2,3,4 };
float b[4] = { 5,6,7,8 };
float res[4];

__m128 A = _mm_loadu_ps(a);
__m128 B = _mm_loadu_ps(b);
__m128 RES = _mm_mul_ps(A, B);

_mm_storeu_ps(res, RES);

查看你写的SIMD代码

在这里你可以把你手写的代码粘贴上去，并且选择合适的编译器指令和编译器。事实上，这个网站编写本科生的大作业完全没问题。

另外，你可以试试用gcc，clang或intel-cpp编译器的时候打开-fopenmp，然后把你想要优化的部分加入杂注#pragma omp simd，看看你的代码和OpenMP相比哪个优化的更好。下面给出一个例子

void sumVec4(float* __restrict dest, float* __restrict one, float *__restrict another)
{
#pragma omp simd
    for (auto i = 0; i < 4; ++i)
        dest[i] = one[i] + another[i];
}

sumVec4(float*, float*, float*):
        vmovups (%rsi), %xmm0
        vaddps  (%rdx), %xmm0, %xmm0
        vmovups %xmm0, (%rdi)
        ret

读写数据

#include <xmmintrin.h>
#include <stdio.h>

int main() {
    float a[4] = {1.0f, 2.0f, 3.0f, 4.0f};
    __m128 values = _mm_loadu_ps(a); // load from a float array

    float result[4];
    _mm_storeu_ps(result, values); // store to a float array

    for (int i = 0; i < 4; i++) {
        printf("%f ", result[i]);
    }
    printf("\n");

    return 0;
}

你首先应该注意到，啥是ps，p是packed，就是向量化的，与之相对的是标量的scalar，至于第二个s，意思就是single，意思是单浮点数，如果是双浮点数就是d。

此外你注意到，为什么load前面有一个u，这代表着你的输入/输出到的是未对齐的内存。如果你保证你输入/输出的对齐方式和__m128相同，那么可以用没u版本的load与store.

以下是C++栈上分配对齐内存的方式

alignas(alignof(__m128)) float fvec4_data[4];

对于C++而言，堆上分配对齐内存，不同编译器有不同的实现方式，对于gcc系列的编译器而言，可以使用C++17定义的std::aligned_alloc，对于其他编译器而言，比如MSVC，就有自己的一套实现方式。

• MSVC 不支持 std::aligned_alloc 函数，因为 C11 标准中指定的 aligned_alloc 函数的实现方式与 Microsoft 实现的 free 函数不兼容。C11 标准要求 free 函数必须能够处理高度对齐的分配，但 Microsoft 的实现并不满足这一要求1。

• 作为替代，MSVC 提供了 _aligned_malloc 函数（需要使用 _aligned_free 函数来释放内存）。你可以使用这个函数来分配对齐的内存。

那么为了实现跨平台的aligned_alloc和aligned_free我们使用宏来实现这两个功能。

#include <cstdlib>

#if defined(_MSC_VER)
    #include <malloc.h>
    void* aligned_alloc(size_t alignment, size_t size) {
        return _aligned_malloc(size, alignment);
    }
    void aligned_free(void* ptr) {
        _aligned_free(ptr);
    }
#elif defined(__MINGW32__) || defined(__MINGW64__)
    #include <malloc.h>
    void* aligned_alloc(size_t alignment, size_t size) {
        return __mingw_aligned_malloc(size, alignment);
    }
    void aligned_free(void* ptr) {
        __mingw_aligned_free(ptr);
    }
#else
    void aligned_free(void* ptr) {
        free(ptr);
    }
#endif

但是如果你使用SIMD，他有一个自带的分配对齐内存的分配器_mm_alloc。使用方式类似于C的malloc和free.

#include <stdio.h>
#include <xmmintrin.h>

int main() {
    // 分配一个大小为1024字节且对齐到32字节边界的内存块
    void* ptr = _mm_malloc(1024, 32);

    // 检查内存分配是否成功
    if (ptr == NULL) {
        printf("Memory allocation failed!\n");
        return -1;
    }

    // 释放内存
    _mm_free(ptr);

    return 0;
}

计算简单四则运算

四则运算是最简单SIMD指令了，只需要加载+加减乘除+输出就好。

#include <xmmintrin.h>
#include <iostream>

int main() {
    float a[4] = {1.0f, 2.0f, 3.0f, 4.0f};
    float b[4] = {5.0f, 6.0f, 7.0f, 8.0f};
    __m128 a_values = _mm_loadu_ps(a);
    __m128 b_values = _mm_loadu_ps(b);

    // 加法
    __m128 add_result = _mm_add_ps(a_values, b_values);
    float add[4];
    _mm_storeu_ps(add, add_result);
    std::cout << "Addition result: ";
    for (int i = 0; i < 4; i++) {
        std::cout << add[i] << " ";
    }
    std::cout << std::endl;

    // 减法
    __m128 sub_result = _mm_sub_ps(a_values, b_values);
    float sub[4];
    _mm_storeu_ps(sub, sub_result);
    std::cout << "Subtraction result: ";
    for (int i = 0; i < 4; i++) {
        std::cout << sub[i] << " ";
    }
    std::cout << std::endl;

    // 乘法
    __m128 mul_result = _mm_mul_ps(a_values, b_values);
    float mul[4];
    _mm_storeu_ps(mul, mul_result);
    std::cout << "Multiplication result: ";
    for (int i = 0; i < 4; i++) {
        std::cout << mul[i] << " ";
    }
    std::cout << std::endl;

    // 除法
    __m128 div_result = _mm_div_ps(a_values, b_values);
    float div[4];
    _mm_storeu_ps(div, div_result);
    std::cout << "Division result: ";
    for (int i = 0; i < 4; i++) {
        std::cout << div[i] << " ";
    }
    std::cout << std::endl;

    return 0;
}

点乘

点乘需要的是两个向量读入并且相乘，然后横向地加起来。

这里我们直观地使用hadd这个函数，它的原理是这样的

__m128 _mm_hadd_ps (__m128 a, __m128 b)
dst[31:0] := a[63:32] + a[31:0]
dst[63:32] := a[127:96] + a[95:64]
dst[95:64] := b[63:32] + b[31:0]
dst[127:96] := b[127:96] + b[95:64]

那么我们使用两次hadd就可以让四个横向的值加在一起了。

#include <pmmintrin.h> // SSE3 support hadd
#include <xmmintrin.h>
#include <iostream>

int main() {
    float a[4] = {1.0f, 2.0f, 3.0f, 4.0f};
    float b[4] = {5.0f, 6.0f, 7.0f, 8.0f};
    __m128 a_values = _mm_loadu_ps(a);
    __m128 b_values = _mm_loadu_ps(b);

    // 点积
    __m128 mul_result = _mm_mul_ps(a_values, b_values);
    __m128 sum_result = _mm_hadd_ps(mul_result, mul_result);
    sum_result = _mm_hadd_ps(sum_result, sum_result);
    float dot_product;
    _mm_store_ss(&dot_product, sum_result);

    std::cout << "Dot product: " << dot_product << std::endl;

    return 0;
}

但是我在这里设置点乘的教学还有另外一个目的，就是使用shuffle函数，这是个十分有用的函数。

shuffle的作用大概是这样的，通过imm8这个无符号数的后8位控制新向量的某个数来自a和b的哪个数，因为是单浮点数，那么只需要2位（即最大为4）就可以指示。

__m128 _mm_shuffle_ps (__m128 a, __m128 b, unsigned int imm8)
// 不用担心，这里只是个宏而已
DEFINE SELECT4(src, control) {
    CASE(control[1:0]) OF
    0:  tmp[31:0] := src[31:0]
    1:  tmp[31:0] := src[63:32]
    2:  tmp[31:0] := src[95:64]
    3:  tmp[31:0] := src[127:96]
    ESAC
    RETURN tmp[31:0]
}
dst[31:0] := SELECT4(a[127:0], imm8[1:0])
dst[63:32] := SELECT4(a[127:0], imm8[3:2])
dst[95:64] := SELECT4(b[127:0], imm8[5:4])
dst[127:96] := SELECT4(b[127:0], imm8[7:6])

SSE头文件内有_MM_SHUFFLE来方便的产生这样的数字（如果宏内的数字都是常量表达式，这不会产生任何性能损失）。但是这个宏和直观违背，它的顺序是相反的！我们可以使用一个新的宏来代替它。

#include <stdio.h>
#include <xmmintrin.h>

#define X_SHUFFLE(fa, sa, fb, sb) _MM_SHUFFLE(sb, fb, sa, fa)

int main() {
    float values[4] = { 0.0f, 1.0f, 2.0f, 3.0f };

    __m128 vec = _mm_loadu_ps(values);

    // 使用_mm_shuffle对寄存器中的元素进行重新排列
    vec = _mm_shuffle_ps(vec, vec, X_SHUFFLE(1, 1, 2, 3));

    _mm_storeu_ps(values, vec);

    printf("%f %f %f %f\n", values[0], values[1], values[2], values[3]); // 1, 1, 2, 3

    return 0;
}

那么既然可以shuffle，我们计算点乘的思路就是这样的，我们先计算出两个向量的乘法，然后一直做“前一半向量加后一半向量”这一操作，直到得到的向量长度为1。折半操作的次数是\log_2 n，所以我们的次数是\log_2 4 = 2.

#include <stdio.h>
#include <xmmintrin.h>

#define X_SHUFFLE(fa, sa, fb, sb) _MM_SHUFFLE(sb, fb, sa, fa)

int main() {
    float values[4] = { 0.0f, 1.0f, 2.0f, 3.0f };
    float another_values[5] = { 1.f, 2.f, 3.f, 4.f };

    auto vec = _mm_loadu_ps(values);
    auto vec_another = _mm_loadu_ps(another_values);

    // 点乘
    auto ans = _mm_mul_ps(vec, vec_another);

    // or ans = _mm_add_ps(ans, _mm_movehl_ps(ans, ans));
    ans = _mm_add_ps(ans, _mm_shuffle_ps(ans, ans, X_SHUFFLE(2, 3, 0, 0)));
    ans = _mm_add_ps(ans, _mm_shuffle_ps(ans, ans, X_SHUFFLE(1, 0, 0, 0)));

    // convert to float32
    printf("%f", _mm_cvtss_f32(ans));
    return 0;
}

叉乘

考虑到三维向量的叉乘对SIMD并不友好（更喜欢2的指数！），我们使用四个浮点数表示一个三维向量（类似于齐次坐标的）。

float vec_one[4] {1.f, 2.f, 3.f, 4.f}, vec_another[4]{5.f, 6.f, 7.f, 8.f};

然后我们就需要计算叉乘了

根据叉乘的定义

注意到\vec i并没有对应x，所以我们需要shuffle一下再shuffle回来。

#include <stdio.h>
#include <xmmintrin.h>

int main() {
    float vec1[4] = {1.0f, 2.0f, 3.0f, 0.0f};
    float vec2[4] = {4.0f, 5.0f, 6.0f, 0.0f};

    __m128 v1 = _mm_loadu_ps(vec1);
    __m128 v2 = _mm_loadu_ps(vec2);

    // 计算两个向量的叉积
    __m128 v1_yzx = _mm_shuffle_ps(v1, v1, _MM_SHUFFLE(3, 0, 2, 1));
    __m128 v2_yzx = _mm_shuffle_ps(v2, v2, _MM_SHUFFLE(3, 0, 2, 1));
    __m128 c1 = _mm_mul_ps(v1_yzx, v2);
    __m128 c2 = _mm_mul_ps(v1, v2_yzx);
    __m128 cross_product = _mm_sub_ps(c1, c2);
    cross_product = _mm_shuffle_ps(cross_product, cross_product, _MM_SHUFFLE(3, 0, 2, 1));

    float values[4];
    _mm_storeu_ps(values, cross_product);

    printf("%f %f %f\n", values[0], values[1], values[2]);

    return 0;
}

预读取

对于一个读取连续内存，我们希望接下来读取的内容已经被存储在了较高级别的cache中以优化IO性能。这时我们需要预读取指令。一个预读取指令由希望读取的内存地址和希望读取到的cache级别构成（这只是一个HINT，具体优化取决于编译器）。以下是一个例子。

#include <stdio.h>
#include <xmmintrin.h>

int main() {
    #define SIZE 1000000
    float values[SIZE];

    for (int i = 0; i < SIZE; i++) {
        values[i] = (float)i;
    }

    // 计算数组元素的和
    float sum = 0.0f;
    for (int i = 0; i < SIZE; i++) {
        // 预取下一个缓存行中的数据
        _mm_prefetch((const char*)(&values[i + 64]), _MM_HINT_T0); // 64B-CACHE-LINE

        sum += values[i];
    }

    printf("Sum: %f\n", sum);

    return 0;
}

_mm_prefetch指令的延迟（即执行所需的时钟周期数）取决于许多因素，包括处理器型号、内存子系统的性能以及数据所在的内存位置等。通常情况下，_mm_prefetch指令的延迟会比许多其他SSE指令要长，因为它需要从内存中获取数据。

然而，_mm_prefetch指令通常不会阻塞后续指令的执行。这意味着，在调用_mm_prefetch指令后，处理器可以继续执行后续指令，而不必等待预取操作完成。因此即使_mm_prefetch指令本身的延迟较长，它也不一定会对程序的总体性能产生负面影响。

存储流

考虑到prefetch对cache的要求，所以我们尽量不要将输出的短期内不再访问的数据放到缓存（CPU会自动把要输出的数据写入到cache以防止短期内再次访问）。所以一个直接输出到内存的指令是必要的。

下面是一个例子（非典型，只说明如何使用）。

#include <stdio.h>
#include <xmmintrin.h>

int main() {
    float values[4] = {1.0f, 2.0f, 3.0f, 4.0f};

    __m128 vec = _mm_loadu_ps(values);

    __m128 ones = _mm_set1_ps(1.0f);
    vec = _mm_add_ps(vec, ones);

    // 使用_mm_stream_ps将结果直接写入内存
    _mm_stream_ps(values, vec);

    printf("%f %f %f %f\n", values[0], values[1], values[2], values[3]);
    return 0;
}