内积实验：CPU、FPU 与 SSE*

• 下载 inner_demo.zip - 9.6 KB
• 下载 inner_src.zip - 4.38 KB

引言

内积（或点积、标量积）运算是数字信号处理领域中的主要运算。它被广泛应用于傅里叶变换（FFT、DCT）、小波分析、滤波运算等。借助 SSE 技术的进步，您可以并行化此运算，以便立即对多个数字执行乘法和加法运算。但是，在计算中选择哪种精度，整数、浮点数还是双精度浮点数？ 在本文中，我演示了对 shorts、ints、floats、doubles 进行内积运算，分别使用 CPU 和 SSE/SSE2/SSE3 优化的版本。

背景（可选）

您需要了解内积运算和 SSE 技术。我喜欢维基百科，因为它提供了所有问题的答案，请参阅 内积。 简而言之，您需要获取 2 个 std::vector 数组（floats、ints、shorts、doubles），长度相等，逐元素相乘，然后将结果向量中的条目相加，生成一个数字。对于 SSE 编程，codeproject 上有一篇不错的文章 SSE 编程入门。

使用代码

只需运行控制台应用程序，并将数组的长度作为第一个参数传递给它。它会创建 2 个相同长度且包含随机条目的向量，并计算它们的内积，打印结果和处理时间，分别针对 chars、shorts/shorts SSE2、ints、floats/floats SSE/floats SSE3、doubles/doubles SSE2。

>inner.exe 5000000
 chars          processing time: 13 ms
 -119
 shorts         processing time: 17 ms
 3241
 shorts sse2    processing time: 6 ms
 3241
 ints           processing time: 16 ms
 786540
 floats         processing time: 30 ms
 1339854
 sse2 intrin    processing time: 13 ms
 1339854
 sse3 assembly  processing time: 12 ms
 1339854
 doubles        processing time: 22 ms
 1107822
 doubles sse2   processing time: 25 ms
 1107822

这是对大小为 5000000 的 2 个向量的内积提供的结果。每种类型后的第二行表示运算的四舍五入结果。我是在我的 AMD turion 64 2.2Ghz 处理器上运行的。我们可以看到，使用 SSE2 的 short 精度是最快的。SSE3 中的新型 haddps 指令允许比 SSE2 更快的处理。使用 FPU 执行的双精度运算优于 FPU 浮点数。但是，SSE2 优化双精度浮点数会降低计算速度，与浮点数相比，使用 SSE 可以提高超过 2 倍的速度。

如果您对图像数据执行 DSP，请使用 SSE2 和 shorts，使用整数滤波器，或使用浮点数。不要优化双精度浮点数，这会导致性能下降。如果您不打算使用 SSE 优化，请使用双精度浮点数，它比浮点数更快。

我没有在其他处理器上尝试过，如果您的结果不同，请告诉我。下面是执行内积运算的代码

在浮点数上使用 FPU

std::vector<float< v1;
std::vector<float> v2;
v1.resize(size);
v2.resize(size);
float* pv1 = &v1[0];
float* pv2 = &v2[0];

for (unsigned int i = 0; i < size; i++) {
        pv1[i] = float(rand() % 64 - 32);
        pv2[i] = float(rand() % 64 - 32);
}

float sum = 0;
for (unsigned int i = 0; i < size; i++)
    sum += pv1[i] * pv2[i];        
wprintf(L" %d\n", (int)sum);

SSE2 优化的 shorts

short sse2_inner_s(const short* p1, const short* p2, unsigned int size)
{
        __m128i* mp1 = (__m128i *)p1;
        __m128i* mp2 = (__m128i *)p2;
        __m128i mres = _mm_set_epi16(0, 0, 0, 0, 0, 0, 0, 0);
        
        for(unsigned int i = 0; i < size/8; i++) {                 
                __m128i mtmp = _mm_mullo_epi16(_mm_loadu_si128(mp1),
                    _mm_loadu_si128(mp2)); 
                mres = _mm_add_epi16(mres, mtmp);
                mp1++;
                mp2++;
        }

        short res[8];
        __m128i* pmres = (__m128i *)res;
        _mm_storeu_si128(pmres, mres);

        return res[0]+res[1]+res[2]+res[3]+res[4]+res[5]+res[6]+res[7];
}

SSE 优化的浮点数

float sse_inner(const float* a, const float* b, unsigned int size)
{
        float z = 0.0f, fres = 0.0f;
        __declspec(align(16)) float ftmp[4] = { 0.0f, 0.0f, 0.0f, 0.0f };
        __m128 mres;

        if ((size / 4) != 0) {
                mres = _mm_load_ss(&z);
                for (unsigned int i = 0; i < size / 4; i++)
                        mres = _mm_add_ps(mres, _mm_mul_ps(_mm_loadu_ps(&a[4*i]),
                        _mm_loadu_ps(&b[4*i])));

                //mres = a,b,c,d
                __m128 mv1 = _mm_movelh_ps(mres, mres);     //a,b,a,b
                __m128 mv2 = _mm_movehl_ps(mres, mres);     //c,d,c,d
                mres = _mm_add_ps(mv1, mv2);                //res[0],res[1]

                _mm_store_ps(ftmp, mres);                

                fres = ftmp[0] + ftmp[1];
        }

        if ((size % 4) != 0) {
                for (unsigned int i = size - size % 4; i < size; i++)
                        fres += a[i] * b[i];
        }

        return fres;
}

SSE3 优化的浮点数

float sse3_inner(const float* a, const float* b, unsigned int size)
{
        float z = 0.0f, fres = 0.0f;
        
        if ((size / 4) != 0) {
                const float* pa = a;
                const float* pb = b;
                __asm {
                        movss   xmm0, xmmword ptr[z]
                }
                for (unsigned int i = 0; i < size / 4; i++) {
                        __asm {
                                mov     eax, dword ptr[pa]
                                mov     ebx, dword ptr[pb]
                                movups  xmm1, [eax]
                                movups  xmm2, [ebx]
                                mulps   xmm1, xmm2
                                addps   xmm0, xmm1
                        }
                        pa += 4;
                        pb += 4;
                }  
                __asm {
                        haddps  xmm0, xmm0
                        haddps  xmm0, xmm0
                        movss   dword ptr[fres], xmm0                        
                }                
        }

        return fres;
}

SSE 优化的双精度浮点数

double sse_inner_d(const double* a, const double* b, unsigned int size)
{
        double z = 0.0, fres = 0.0;
        __declspec(align(16)) double ftmp[2] = { 0.0, 0.0 };
        __m128d mres;
        
        if ((size / 2) != 0) {
                mres = _mm_load_sd(&z);
                for (unsigned int i = 0; i < size / 2; i++)
                        mres = _mm_add_pd(mres, _mm_mul_pd(_mm_loadu_pd(&a[2*i]),
                        _mm_loadu_pd(&b[2*i])));                

                _mm_store_pd(ftmp, mres);                

                fres = ftmp[0] + ftmp[1];
        }

        if ((size % 2) != 0) {
                for (unsigned int i = size - size % 2; i < size; i++)
                        fres += a[i] * b[i];
        }

        return fres;
}