对SSE3矩阵乘法优化的认识

Question

参考http://blogs.msdn.com/b/xiangfan/archive/2009/04/28/optimize-your-code-matrix-multiplication.aspx。

template<>
void SeqMatrixMult4(int size, float** m1, float** m2, float** result)
{
    Transpose(size, m2);
    for (int i = 0; i < size; i++) {
        for (int j = 0; j < size; j++) {
            __m128 c = _mm_setzero_ps();

            for (int k = 0; k < size; k += 4) {
                c = _mm_add_ps(c, _mm_mul_ps(_mm_load_ps(&m1[i][k]), _mm_load_ps(&m2[j][k])));
            }
            c = _mm_hadd_ps(c, c);
            c = _mm_hadd_ps(c, c);
            _mm_store_ss(&result[i][j], c);
        }
    }
    Transpose(size, m2);
}

为什么在内部最易循环之后还会有两个_mm_hadd_ps(c, c)？为了验证我的理解:这段代码从m1加载了4个浮点数，从m2加载了另外4个浮点数，然后将它们相乘，生成4个浮点数(__m128)。然后我将它们相加为c (此时，它仍然是4个浮点数？)。那么在for循环i之后，hadd这个结果会出现两次吗？那有什么用？

我的代码稍加重写就会产生错误的结果。

long long start, end;
__m128 v1, v2, vMul, vRes;
vRes = _mm_setzero_ps();

start = wall_clock_time();
transpose_matrix(m2);
for (int i = 0; i < SIZE; i++) {
    for (int j = 0; j < SIZE; j++) {
        float tmp = 0;
        for (int k = 0; k < SIZE; k+=4) {
            v1 = _mm_load_ps(&m1[i][k]);
            v2 = _mm_load_ps(&m2[j][k]);
            vMul = _mm_mul_ps(v1, v2);

            vRes = _mm_add_ps(vRes, vMul);
        }
        vRes = _mm_hadd_ps(vRes, vRes);
        _mm_store_ss(&result[i][j], vRes);
    }
}
end = wall_clock_time();
fprintf(stderr, "Optimized Matrix multiplication took %1.2f seconds\n", ((float)(end - start))/1000000000);

// reverse the transposition
transpose_matrix(m2);

Answer 1

haddps并不是将所有四个元素相加在一个向量中。需要两个haddps指令才能得到完全的水平和。

如果对向量{c0,c1,c2,c3}的元素进行编号，第一个haddps将生成{c0+c1, c2+c3, c0+c1, c2+c3}。第二种是产生{c0+c1+c2+c3, <same thing in the other lanes>}。