【发布时间】：2015-09-13 14:23:21
【问题描述】：
我正在编写一个复杂的模拟程序，并且发现最耗时的例程是将四向量 (float4) 与 4x4 矩阵相乘的例程。我需要在几台或多或少旧的计算机上运行这个程序。这就是为什么我尝试在以下代码中检查此类操作的 SIMD 功能：
//#include <xmmintrin.h> // SSE
//#include <pmmintrin.h> // SSE3
//#include <nmmintrin.h> // SSE4.2
#include <immintrin.h> // AVX
#include <iostream>
#include <ctime>
#include <string>
using namespace std;
// 4-vector.
typedef struct
{
float x;
float y;
float z;
float w;
}float4;
// typedef to simplify the pointer of function notation.
typedef void(*Function)(float4&,const float4*,const float4&);
float dot( const float4& in_A, const float4& in_x )
{
return in_A.x*in_x.x + in_A.y*in_x.y + in_A.z*in_x.z + in_A.w*in_x.w; // 7 FLOPS
}
void A_times_x( float4& out_y, const float4* in_A, const float4& in_x )
{
out_y.x = dot(in_A[0], in_x); // 7 FLOPS
out_y.y = dot(in_A[1], in_x); // 7 FLOPS
out_y.z = dot(in_A[2], in_x); // 7 FLOPS
out_y.w = dot(in_A[3], in_x); // 7 FLOPS
}
void A_times_x_SSE( float4& out_y, const float4* in_A, const float4& in_x )
{
// Load matrix A and vector x into SSE registers
__m128 x = _mm_load_ps((const float*)&in_x); // load/store are almost = 0 FLOPS
__m128 A0 = _mm_load_ps((const float*)(in_A + 0));
__m128 A1 = _mm_load_ps((const float*)(in_A + 1));
__m128 A2 = _mm_load_ps((const float*)(in_A + 2));
__m128 A3 = _mm_load_ps((const float*)(in_A + 3));
// Transpose the matrix and re-order the vector.
_MM_TRANSPOSE4_PS( A0,A1,A2,A3 );
__m128 u1 = _mm_shuffle_ps(x,x, _MM_SHUFFLE(0,0,0,0));
__m128 u2 = _mm_shuffle_ps(x,x, _MM_SHUFFLE(1,1,1,1));
__m128 u3 = _mm_shuffle_ps(x,x, _MM_SHUFFLE(2,2,2,2));
__m128 u4 = _mm_shuffle_ps(x,x, _MM_SHUFFLE(3,3,3,3));
// Multiply each matrix row with the vector x
__m128 m0 = _mm_mul_ps(A0, u1); // 4 FLOPS
__m128 m1 = _mm_mul_ps(A1, u2); // 4 FLOPS
__m128 m2 = _mm_mul_ps(A2, u3); // 4 FLOPS
__m128 m3 = _mm_mul_ps(A3, u4); // 4 FLOPS
// Using HADD, we add four floats at a time
__m128 sum_01 = _mm_add_ps(m0, m1); // 4 FLOPS
__m128 sum_23 = _mm_add_ps(m2, m3); // 4 FLOPS
__m128 result = _mm_add_ps(sum_01, sum_23); // 4 FLOPS
// Finally, store the result
_mm_store_ps((float*)&out_y, result);
}
void A_times_x_SSE3( float4& out_y, const float4* in_A, const float4& in_x )
{
// Should be 4 (SSE) x 4 (ALU) = 16 times faster than scalar.
// Load matrix A and vector x into SSE registers
__m128 x = _mm_load_ps((const float*)&in_x); // load/store are almost = 0 FLOPS
__m128 A0 = _mm_load_ps((const float*)(in_A + 0));
__m128 A1 = _mm_load_ps((const float*)(in_A + 1));
__m128 A2 = _mm_load_ps((const float*)(in_A + 2));
__m128 A3 = _mm_load_ps((const float*)(in_A + 3));
// Multiply each matrix row with the vector x
__m128 m0 = _mm_mul_ps(A0, x); // 4 FLOPS
__m128 m1 = _mm_mul_ps(A1, x); // 4 FLOPS
__m128 m2 = _mm_mul_ps(A2, x); // 4 FLOPS
__m128 m3 = _mm_mul_ps(A3, x); // 4 FLOPS
// Using HADD, we add four floats at a time
__m128 sum_01 = _mm_hadd_ps(m0, m1); // 4 FLOPS
__m128 sum_23 = _mm_hadd_ps(m2, m3); // 4 FLOPS
__m128 result = _mm_hadd_ps(sum_01, sum_23); // 4 FLOPS
// Finally, store the result
_mm_store_ps((float*)&out_y, result);
}
void A_times_x_SSE4( float4& out_y, const float4* in_A, const float4& in_x ) // 28 FLOPS
{
// Should be 4 (SSE) x 4 (ALU) = 16 times faster than scalar.
// Load matrix A and vector x into SSE registers
__m128 x = _mm_load_ps((const float*)&in_x); // load/store are almost = 0 FLOPS
__m128 A0 = _mm_load_ps((const float*)(in_A + 0));
__m128 A1 = _mm_load_ps((const float*)(in_A + 1));
__m128 A2 = _mm_load_ps((const float*)(in_A + 2));
__m128 A3 = _mm_load_ps((const float*)(in_A + 3));
// Multiply each matrix row with the vector x
__m128 m0 = _mm_dp_ps(A0, x, 0xFF); // 4 FLOPS
__m128 m1 = _mm_dp_ps(A1, x, 0xFF); // 4 FLOPS
__m128 m2 = _mm_dp_ps(A2, x, 0xFF); // 4 FLOPS
__m128 m3 = _mm_dp_ps(A3, x, 0xFF); // 4 FLOPS
// Using HADD, we add four floats at a time
__m128 mov_01 = _mm_movelh_ps(m0, m1); // 4 FLOPS
__m128 mov_23 = _mm_movelh_ps(m2, m3); // 4 FLOPS
__m128 result = _mm_shuffle_ps(mov_01, mov_23, _MM_SHUFFLE(2, 0, 2, 0)); // 4 FLOPS
// Finally, store the result
_mm_store_ps((float*)&out_y, result);
}
void A_times_x_AVX( float4& out_y, const float4* in_A, const float4& in_x )
{
// Load matrix A and vector x into SSE registers
__m128 x = _mm_load_ps((const float*)&in_x); // load/store are almost = 0 FLOPS
__m256 xx = _mm256_castps128_ps256(x);
xx = _mm256_insertf128_ps(xx,x,1);
__m256 A0 = _mm256_load_ps((const float*)(in_A + 0));
__m256 A2 = _mm256_load_ps((const float*)(in_A + 2));
// Multiply each matrix row with the vector x
__m256 m0 = _mm256_mul_ps(A0, xx); // 4 FLOPS
__m256 m2 = _mm256_mul_ps(A2, xx); // 4 FLOPS
// Using HADD, we add four floats at a time
__m256 sum_00 = _mm256_hadd_ps(m0, m2); // 4 FLOPS
/*__m128 sum_10 = _mm256_extractf128_ps(sum_00,0);
__m128 sum_01 = _mm256_extractf128_ps(sum_00,1);
__m128 result = _mm_hadd_ps(sum_10, sum_01); // 4 FLOPS
// Finally, store the result
_mm_store_ps((float*)&out_y, result);*/
// Finally, store the result (no temp variable: direct HADD, this avoid to copy from ALU128 to ALU256)
_mm_store_ps((float*)&out_y, _mm_hadd_ps(_mm256_extractf128_ps(sum_00,0),
_mm256_extractf128_ps(sum_00,1)));
}
void test_function ( Function f, string simd, unsigned int imax )
{
float4 Y;
float4 X1 = {0.5,1,0.2,0.7};
float4 X2 = {0.7,1,0.2,0.5};
float4 X3 = {0.5,0.2,1,0.7};
float4 X4 = {1,0.7,0.2,0.5};
float4 A[4] = {{0.5,1,0.2,0.7},
{0.6,0.4,0.1,0.8},
{0.3,0.8,0.2,0.5},
{1,0.4,0.6,0.9}};
clock_t tstart = clock();
for( unsigned int i=0 ; i<imax ; i++ )
for( unsigned long int j=0 ; j<250000000 ; j++ )
// Avoid for loop over long long, it is 2 times slower !
{
// Function pointer give a real call, whether the direct
// call is inlined and thus results are overestimated.
f( Y,A,X1 );
f( Y,A,X2 );
f( Y,A,X3 );
f( Y,A,X4 );
}
clock_t tend = clock();
double diff = static_cast<double>(tend - tstart) * 1e-3;
cout << "Time (" << simd << ") = " << diff << " s" << endl;
cout << "Nops (" << simd << ") = " << (double) imax << ".10^9" << endl;
cout << "Power (" << simd << ") = " << (double) imax * 28. / diff << " GFLOPS" << endl; // 28 FLOPS for std.
cout << endl;
}
int main ( int argc, char *argv[] )
{
test_function ( &A_times_x ,"std" , 1 );
test_function ( &A_times_x_SSE ,"SSE" , 2 );
test_function ( &A_times_x_SSE3,"SSE3", 3 );
test_function ( &A_times_x_SSE4,"SSE4", 1 );
test_function ( &A_times_x_AVX ,"AVX" , 3 );
return 0;
}
我对此类问题的改进有一些麻烦。运行代码时，我得到以下结果（Intel Core i5 4670K、3.4GHz、Haswell、Codeblock+MinGW 编译器使用 -O2 -march=corei7-avx）：
Time (std) = 6.287 s
Nops (std) = 1.10^9
Power (std) = 4.45363 GFLOPS
Time (SSE) = 6.661 s
Nops (SSE) = 2.10^9
Power (SSE) = 8.40715 GFLOPS
Time (SSE3) = 8.361 s
Nops (SSE3) = 3.10^9
Power (SSE3) = 10.0466 GFLOPS
Time (SSE4) = 6.131 s
Nops (SSE4) = 1.10^9
Power (SSE4) = 4.56695 GFLOPS
Time (AVX) = 8.767 s
Nops (AVX) = 3.10^9
Power (AVX) = 9.58138 GFLOPS
我的问题如下：
这是否可以进一步提高性能/加速？它应该是 SSE 为 x4（最大值），AVX 为 x8。
为什么 AVX 不比 SSE3 快？
对于那些说：“停止使用你的东西，使用英特尔数学内核库”的人，我回答：“我不会，因为我想要一个小的可执行文件，而且我只需要在这种特定情况下使用 SIMD，而不是其他地方”;-)
【问题讨论】：
-
您使用了什么编译器，您启用了哪些优化开关？
-
问题是你花更多的时间来打乱数据而不是做有用的算术。由于它是大小为 4 的向量 x 矩阵乘法，因此计算量太少而无法单独优化。您可能需要扩大范围并考虑将数据布局更改为更适合 SIMD 的内容。
-
我在 Haswell 上使用 MSVC2013 测试了您的代码，结果不同。我将所有测试设置为 10x10^9，std 为 11.9s，所有其他测试为 2.95s（std 正好慢 x4）。此外，std 已完全优化，因此我在结果中添加了 volatile。另外，似乎所有的调用都是内联的，所以没有调用，这可以解释当你实际调用时速度很慢。
-
我感觉 OP 要么使用旧的编译器，要么没有启用优化，因为我在这里尝试过的每个编译器都会优化所有循环，使测试结果毫无意义。
-
我将 MinGW 与 Codeblock 一起使用（我知道……很遗憾），我的编译器标志是：-Wall -O2 -march=corei7-avx。今晚我将尝试使用 MSVC2013 使用英特尔 C++ 编译器检查这一点。感谢您的有用评论！