SIMD 运算
SIMD(Single Instruction Multiple Data)即单指令流多数据流，是一种采用一个控制器来控制多个处理器，同时对一组数据（又称“数据向量”）中的每一个分别执行相同的操作从而实现空间上的并行性的技术。简单来说就是一个指令能够同时处理多个数据。
SSE/AVX指令主要定义于以下一些头文件中：
: SSE, 支持同时对4个32位单精度浮点数的操作。 : SSE 2, 支持同时对2个64位双精度浮点数的操作。 : SSE 3, 支持对SIMD寄存器的水平操作(horizontal operation)，如hadd, hsub等...。 : SSSE 3, 增加了额外的instructions。 : SSE 4.1, 支持点乘以及更多的整形操作。 : SSE 4.2, 增加了额外的instructions。 : AVX, 支持同时操作8个单精度浮点数或4个双精度浮点数。
一步一步使用simd
一个简单的例子
#include <iostream>
#include <immintrin.h>
int main(int argc, const char** argv)
{
double input1[4] = {1, 1, 1, 1};
double input2[4] = {1, 2, 3, 4};
double result[4];
std::cout << "address of input1: " << input1 << std::endl;
std::cout << "address of input2: " << input2 << std::endl;
__m256d a = _mm256_load_pd(input1);
__m256d b = _mm256_load_pd(input2);
__m256d c = _mm256_add_pd(a, b);
_mm256_store_pd(result, c);
std::cout << result[0] << " " << result[1] << " " << result[2] << " " << result[3] << std::endl;
return 0;
}
使用g++编译运行.
g++ -std=c++11 -march=native main.cpp -o simd
在PC上默认参数可以看到是\(march=x86\_64\).
g++ -Q --help=target
发现会出现crash。这里使用了AVX指令集运算，要求32字节对齐，我们可以把这两个输入变量的地址除以32，看是否能够整除。结果发现此时的地址都不能整除。当然，其实直接看倒数第二位是否是偶数即可，是偶数就可以被32整除，是奇数则不能被32整除。
对于局部变量内存地址是编译器决定的，需要告诉编译器声明变量需要内存对齐，gcc的语法为__attribute__((aligned(32)))，MSVC的语法为 __declspec(align(32))。
__attribute__ ((aligned (32))) double input1[4] = {1, 1, 1, 1};
__attribute__ ((aligned (32))) double input2[4] = {1, 2, 3, 4};
__attribute__ ((aligned (32))) double result[4];
此时编译运行，发现程序可以正常运行，并且内存地址都是32字节对齐。
封装
Eigen中会将向量封装起来，此时可以仿照封装一个向量加法。
#include <immintrin.h>
#include <iostream>
class Vector4d {
using aligned_double4 = __attribute__ ((aligned (32))) double[4];
public:
Vector4d(double d1, double d2, double d3, double d4)
{
data[0] = d1;
data[1] = d2;
data[2] = d3;
data[3] = d4;
std::cout << "address: " << data << std::endl;
}
aligned_double4 data;
};
Vector4d operator+ (const Vector4d& v1, const Vector4d& v2)
{
__m256d data1 = _mm256_load_pd(v1.data);
__m256d data2 = _mm256_load_pd(v2.data);
__m256d data3 = _mm256_add_pd(data1, data2);
Vector4d result = {0, 0, 0, 0};
_mm256_store_pd(result.data, data3);
return result;
}
std::ostream& operator<< (std::ostream& o, const Vector4d& v)
{
o << "(" << v.data[0] << ", " << v.data[1] << ", " << v.data[2] << ", " << v.data[3] << ")";
return o;
}
int main(int argc, const char** argv)
{
Vector4d input1 = {1, 1, 1, 1};
Vector4d input2 = {1, 2, 3, 4};
Vector4d result = input1 + input2;
std::cout << result << std::endl;
return 0;
}
此时程序正常运行，但是当我们动态实例化对象就会发现内存对齐出错了，此时程序crash。
int main(int argc, const char** argv)
{
/*Vector4d input1 = {1, 1, 1, 1};
Vector4d input2 = {1, 2, 3, 4};
Vector4d result = input1 + input2;*/
Vector4d* input1 = new Vector4d{1, 1, 1, 1};
Vector4d* input2 = new Vector4d{1, 2, 3, 4};
Vector4d result = *input1 + *input2;
std::cout << result << std::endl;
std::cout << result << std::endl;
return 0;
}
问题的根源在于直接声明的对象是存储在栈上的，其内存地址由编译器在编译时确定，因此预编译指令会生效。但用new动态创建的对象则存储在堆中，其地址在运行时确定。C++的运行时库并不会关心预编译指令声明的对齐方式，需要重载new方法。
void* operator new (std::size_t count) {
void* original = ::operator new(count + 32);
void* aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~size_t(32 - 1)) + 32);
*(reinterpret_cast<void**>(aligned) - 1) = original;
return aligned;
}
void operator delete (void* ptr) {
::operator delete(*(reinterpret_cast<void**>(ptr) - 1));
}
此时发现编译运行，程序又正常运行。
作为成员变量
当我们需要Vector4d作为成员变量的时候，bug又出现了，内存又出现了不对齐的情况。
class Point {
public:
Point(const Vector4d& position) : position(position) {}
Vector4d position;
};
int main(int argc, const char** argv)
{
Vector4d* input1 = new Vector4d{1, 1, 1, 1};
Vector4d* input2 = new Vector4d{1, 2, 3, 4};
Point* point1 = new Point{*input1};
Point* point2 = new Point{*input2};
Vector4d result = point1->position + point2->position;
std::cout << result << std::endl;
delete input1;
delete input2;
delete point1;
delete point2;
return 0;
}
此时需要在Point类中同样重载new和delete方法，将次操作封装为宏，方便其他地方调用。
#define ALIGNED_OPERATOR_NEW \
void* operator new (std::size_t count) { \
void* original = ::operator new(count + 32); \
void* aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~size_t(32 - 1)) + 32); \
*(reinterpret_cast<void**>(aligned) - 1) = original; \
return aligned;\
} \
void operator delete (void* ptr) { \
::operator delete(*(reinterpret_cast<void**>(ptr) - 1)); \
}
完整程序如下：
#include <immintrin.h>
#include <iostream>
#define ALIGNED_OPERATOR_NEW \
void* operator new (std::size_t count) { \
void* original = ::operator new(count + 32); \
void* aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~size_t(32 - 1)) + 32); \
*(reinterpret_cast<void**>(aligned) - 1) = original; \
return aligned;\
} \
void operator delete (void* ptr) { \
::operator delete(*(reinterpret_cast<void**>(ptr) - 1)); \
}
class Vector4d {
using aligned_double4 = __attribute__ ((aligned (32))) double[4];
public:
ALIGNED_OPERATOR_NEW
Vector4d(double d1, double d2, double d3, double d4)
{
data[0] = d1;
data[1] = d2;
data[2] = d3;
data[3] = d4;
std::cout << "address: " << data << std::endl;
}
aligned_double4 data;
};
Vector4d operator+ (const Vector4d& v1, const Vector4d& v2)
{
__m256d data1 = _mm256_load_pd(v1.data);
__m256d data2 = _mm256_load_pd(v2.data);
__m256d data3 = _mm256_add_pd(data1, data2);
Vector4d result = {0, 0, 0, 0};
_mm256_store_pd(result.data, data3);
return result;
}
std::ostream& operator<< (std::ostream& o, const Vector4d& v)
{
o << "(" << v.data[0] << ", " << v.data[1] << ", " << v.data[2] << ", " << v.data[3] << ")";
return o;
}
class Point {
public:
ALIGNED_OPERATOR_NEW
Point(const Vector4d& position) : position(position) {}
Vector4d position;
};
int main(int argc, const char** argv)
{
/*Vector4d input1 = {1, 1, 1, 1};
Vector4d input2 = {1, 2, 3, 4};
Vector4d result = input1 + input2;*/
/*Vector4d* input1 = new Vector4d{1, 1, 1, 1};
Vector4d* input2 = new Vector4d{1, 2, 3, 4};
Vector4d result = *input1 + *input2;*/
Vector4d input1 = {1, 1, 1, 1};
Vector4d input2 = {1, 2, 3, 4};
Point* point1 = new Point{input1};
Point* point2 = new Point{input2};
Vector4d result = point1->position + point2->position;
std::cout << result << std::endl;
delete point1;
delete point2;
return 0;
}
Reference
[1] https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html#expand=18,1155,1286,450,2913,3836,6065,6065&text=_mm_unpacklo_ps
[2] https://en.cppreference.com/w/cpp/memory/new/operator_new
[3] https://zhuanlan.zhihu.com/p/93824687