如何使用 SSE 执行 8 x 8 矩阵运算？答案

【问题标题】：How do I perform 8 x 8 matrix operation using SSE?如何使用 SSE 执行 8 x 8 矩阵运算？
【发布时间】：2012-01-07 07:44:13
【问题描述】：

我最初的尝试是这样的（假设我们想乘）

  __m128 mat[n]; /* rows */
  __m128 vec[n] = {1,1,1,1};
  float outvector[n];
   for (int row=0;row<n;row++) {
       for(int k =3; k < 8; k = k+ 4)
       {
           __m128 mrow = mat[k];
           __m128 v = vec[row];
           __m128 sum = _mm_mul_ps(mrow,v);
           sum= _mm_hadd_ps(sum,sum); /* adds adjacent-two floats */
       }
           _mm_store_ss(&outvector[row],_mm_hadd_ps(sum,sum));
 }

但这显然行不通。我该如何处理？

我应该一次加载 4 个....

另一个问题是：如果我的数组非常大（比如 n = 1000），我怎样才能使它 16 字节对齐？这可能吗？

【问题讨论】：

你期待什么结果？我没有看到任何矩阵，只有向量乘法。还有，3、8、4是哪来的？
@user963889，尺寸没有任何意义。你想做什么？将 8x1 向量或向量数组乘以 8x8 矩阵？
@BrettHale 假设我们有 8x8 倍的向量 8x1。结果我想得到8x1。我被困住了。你们能引导我走向正确的方向吗？谢谢。
@user963889，好的 - 提供了 [8x8] x [8x1] 的答案 ...

标签： c++ sse intrinsics

【解决方案1】：

好的...我将使用行优先矩阵约定。 [m] 的每一行需要 (2) 个 __m128 元素来产生 8 个浮点数。 8x1 向量v 是一个列向量。由于您使用的是 haddps 指令，我假设 SSE3 可用。寻找r = [m] * v：

void mul (__m128 r[2], const __m128 m[8][2], const __m128 v[2])
{
    __m128 t0, t1, t2, t3, r0, r1, r2, r3;

    t0 = _mm_mul_ps(m[0][0], v[0]);
    t1 = _mm_mul_ps(m[1][0], v[0]);
    t2 = _mm_mul_ps(m[2][0], v[0]);
    t3 = _mm_mul_ps(m[3][0], v[0]);

    t0 = _mm_hadd_ps(t0, t1);
    t2 = _mm_hadd_ps(t2, t3);
    r0 = _mm_hadd_ps(t0, t2);

    t0 = _mm_mul_ps(m[0][1], v[1]);
    t1 = _mm_mul_ps(m[1][1], v[1]);
    t2 = _mm_mul_ps(m[2][1], v[1]);
    t3 = _mm_mul_ps(m[3][1], v[1]);

    t0 = _mm_hadd_ps(t0, t1);
    t2 = _mm_hadd_ps(t2, t3);
    r1 = _mm_hadd_ps(t0, t2);

    t0 = _mm_mul_ps(m[4][0], v[0]);
    t1 = _mm_mul_ps(m[5][0], v[0]);
    t2 = _mm_mul_ps(m[6][0], v[0]);
    t3 = _mm_mul_ps(m[7][0], v[0]);

    t0 = _mm_hadd_ps(t0, t1);
    t2 = _mm_hadd_ps(t2, t3);
    r2 = _mm_hadd_ps(t0, t2);

    t0 = _mm_mul_ps(m[4][1], v[1]);
    t1 = _mm_mul_ps(m[5][1], v[1]);
    t2 = _mm_mul_ps(m[6][1], v[1]);
    t3 = _mm_mul_ps(m[7][1], v[1]);

    t0 = _mm_hadd_ps(t0, t1);
    t2 = _mm_hadd_ps(t2, t3);
    r3 = _mm_hadd_ps(t0, t2);

    r[0] = _mm_add_ps(r0, r1);
    r[1] = _mm_add_ps(r2, r3);
}

关于对齐，__m128 类型的变量应该在堆栈上自动对齐。对于动态内存，这不是一个安全的假设。一些 malloc / new 实现可能只返回保证为 8 字节对齐的内存。

内在函数头提供 _mm_malloc 和 _mm_free。在这种情况下，align 参数应为 (16)。

【讨论】：

谢谢。经过两天的工作，我的代码现在看起来和你的很相似......但你的代码很清楚。我学会了。谢谢。

【解决方案2】：

英特尔为大小从 1×1 到 6×6 的矩阵开发了一个Small Matrix Library。应用笔记AP-930 Streaming SIMD Extensions - Matrix Multiplication 详细描述了将两个 6×6 矩阵相乘的算法。这应该可以通过一些努力适应其他尺寸矩阵。

【讨论】：