16

我正在尝试找到一个优化的 C 或汇编程序实现，该函数将两个 4x4 矩阵相乘。该平台是基于 ARM6 或 ARM7 的 iPhone 或 iPod。

目前，我正在使用一种相当标准的方法——只是展开了一点循环。

#define O(y,x) (y + (x<<2))

静态内联 void Matrix4x4MultiplyBy4x4 (float *src1, float *src2, float *dest)
{
    *(dest+O(0,0)) = (*(src1+O(0,0)) * *(src2+O(0,0))) + (*(src1+O(0,1)) * *(src2+O(1,0))) + (*(src1+O(0,2)) * *(src2+O(2,0))) + (*(src1+O(0,3) )) * *(src2+O(3,0)));
    *(dest+O(0,1)) = (*(src1+O(0,0)) * *(src2+O(0,1))) + (*(src1+O(0,1)) * *(src2+O(1,1))) + (*(src1+O(0,2)) * *(src2+O(2,1))) + (*(src1+O(0,3) )) * *(src2+O(3,1)));
    *(dest+O(0,2)) = (*(src1+O(0,0)) * *(src2+O(0,2))) + (*(src1+O(0,1)) * *(src2+O(1,2))) + (*(src1+O(0,2)) * *(src2+O(2,2))) + (*(src1+O(0,3) )) * *(src2+O(3,2)));
    *(dest+O(0,3)) = (*(src1+O(0,0)) * *(src2+O(0,3))) + (*(src1+O(0,1)) * *(src2+O(1,3))) + (*(src1+O(0,2)) * *(src2+O(2,3))) + (*(src1+O(0,3) )) * *(src2+O(3,3)));
    *(dest+O(1,0)) = (*(src1+O(1,0)) * *(src2+O(0,0))) + (*(src1+O(1,1)) * *(src2+O(1,0))) + (*(src1+O(1,2)) * *(src2+O(2,0))) + (*(src1+O(1,3) )) * *(src2+O(3,0)));
    *(dest+O(1,1)) = (*(src1+O(1,0)) * *(src2+O(0,1))) + (*(src1+O(1,1)) * *(src2+O(1,1))) + (*(src1+O(1,2)) * *(src2+O(2,1))) + (*(src1+O(1,3) )) * *(src2+O(3,1)));
    *(dest+O(1,2)) = (*(src1+O(1,0)) * *(src2+O(0,2))) + (*(src1+O(1,1)) * *(src2+O(1,2))) + (*(src1+O(1,2)) * *(src2+O(2,2))) + (*(src1+O(1,3) )) * *(src2+O(3,2)));
    *(dest+O(1,3)) = (*(src1+O(1,0)) * *(src2+O(0,3))) + (*(src1+O(1,1)) * *(src2+O(1,3))) + (*(src1+O(1,2)) * *(src2+O(2,3))) + (*(src1+O(1,3) )) * *(src2+O(3,3)));
    *(dest+O(2,0)) = (*(src1+O(2,0)) * *(src2+O(0,0))) + (*(src1+O(2,1)) * *(src2+O(1,0))) + (*(src1+O(2,2)) * *(src2+O(2,0))) + (*(src1+O(2,3) )) * *(src2+O(3,0)));
    *(dest+O(2,1)) = (*(src1+O(2,0)) * *(src2+O(0,1))) + (*(src1+O(2,1)) * *(src2+O(1,1))) + (*(src1+O(2,2)) * *(src2+O(2,1))) + (*(src1+O(2,3) )) * *(src2+O(3,1)));
    *(dest+O(2,2)) = (*(src1+O(2,0)) * *(src2+O(0,2))) + (*(src1+O(2,1)) * *(src2+O(1,2))) + (*(src1+O(2,2)) * *(src2+O(2,2))) + (*(src1+O(2,3) )) * *(src2+O(3,2)));
    *(dest+O(2,3)) = (*(src1+O(2,0)) * *(src2+O(0,3))) + (*(src1+O(2,1)) * *(src2+O(1,3))) + (*(src1+O(2,2)) * *(src2+O(2,3))) + (*(src1+O(2,3) )) * *(src2+O(3,3)));
    *(dest+O(3,0)) = (*(src1+O(3,0)) * *(src2+O(0,0))) + (*(src1+O(3,1)) * *(src2+O(1,0))) + (*(src1+O(3,2)) * *(src2+O(2,0))) + (*(src1+O(3,3) )) * *(src2+O(3,0)));
    *(dest+O(3,1)) = (*(src1+O(3,0)) * *(src2+O(0,1))) + (*(src1+O(3,1)) * *(src2+O(1,1))) + (*(src1+O(3,2)) * *(src2+O(2,1))) + (*(src1+O(3,3) )) * *(src2+O(3,1)));
    *(dest+O(3,2)) = (*(src1+O(3,0)) * *(src2+O(0,2))) + (*(src1+O(3,1)) * *(src2+O(1,2))) + (*(src1+O(3,2)) * *(src2+O(2,2))) + (*(src1+O(3,3) )) * *(src2+O(3,2)));
    *(dest+O(3,3)) = (*(src1+O(3,0)) * *(src2+O(0,3))) + (*(src1+O(3,1)) * *(src2+O(1,3))) + (*(src1+O(3,2)) * *(src2+O(2,3))) + (*(src1+O(3,3) )) * *(src2+O(3,3)));
};

我会从使用 Strassen 或 Coppersmith-Winograd 算法中受益吗？

4

5 回答 5

36

不，Strassen 或 Coppersmith-Winograd 算法在这里不会有太大的不同。他们开始只为更大的矩阵带来回报。

如果您的矩阵乘法确实是一个瓶颈，您可以使用 NEON SIMD 指令重写算法。这只会对 ARMv7 有所帮助，因为 ARMv6 没有这个扩展。

对于您的情况，我预计编译的 C 代码的速度会提高 3 倍。

编辑：您可以在这里找到 ARM-NEON 中的一个不错的实现：http ://code.google.com/p/math-neon/

对于您的 C 代码，您可以做两件事来加快代码速度：

  1. 不要内联函数。您的矩阵乘法在展开时会生成相当多的代码，而 ARM 只有一个非常小的指令缓存。过多的内联会使您的代码变慢，因为 CPU 将忙于将代码加载到缓存中而不是执行它。

  2. 使用restrict 关键字告诉编译器源指针和目标指针在内存中不重叠。目前，每当写入结果时，编译器都被迫从内存中重新加载每个源值，因为它必须假设源和目标可能重叠甚至指向相同的内存。

于 2009-11-04T14:21:20.920 回答
20

只是吹毛求疵。我想知道为什么人们仍然自愿混淆他们的代码？C 已经很难读了，不需要再添加了。

static inline void Matrix4x4MultiplyBy4x4 (float src1[4][4], float src2[4][4], float dest[4][4])
{
dest[0][0] = src1[0][0] * src2[0][0] + src1[0][1] * src2[1][0] + src1[0][2] * src2[2][0] + src1[0][3] * src2[3][0]; 
dest[0][1] = src1[0][0] * src2[0][1] + src1[0][1] * src2[1][1] + src1[0][2] * src2[2][1] + src1[0][3] * src2[3][1]; 
dest[0][2] = src1[0][0] * src2[0][2] + src1[0][1] * src2[1][2] + src1[0][2] * src2[2][2] + src1[0][3] * src2[3][2]; 
dest[0][3] = src1[0][0] * src2[0][3] + src1[0][1] * src2[1][3] + src1[0][2] * src2[2][3] + src1[0][3] * src2[3][3]; 
dest[1][0] = src1[1][0] * src2[0][0] + src1[1][1] * src2[1][0] + src1[1][2] * src2[2][0] + src1[1][3] * src2[3][0]; 
dest[1][1] = src1[1][0] * src2[0][1] + src1[1][1] * src2[1][1] + src1[1][2] * src2[2][1] + src1[1][3] * src2[3][1]; 
dest[1][2] = src1[1][0] * src2[0][2] + src1[1][1] * src2[1][2] + src1[1][2] * src2[2][2] + src1[1][3] * src2[3][2]; 
dest[1][3] = src1[1][0] * src2[0][3] + src1[1][1] * src2[1][3] + src1[1][2] * src2[2][3] + src1[1][3] * src2[3][3]; 
dest[2][0] = src1[2][0] * src2[0][0] + src1[2][1] * src2[1][0] + src1[2][2] * src2[2][0] + src1[2][3] * src2[3][0]; 
dest[2][1] = src1[2][0] * src2[0][1] + src1[2][1] * src2[1][1] + src1[2][2] * src2[2][1] + src1[2][3] * src2[3][1]; 
dest[2][2] = src1[2][0] * src2[0][2] + src1[2][1] * src2[1][2] + src1[2][2] * src2[2][2] + src1[2][3] * src2[3][2]; 
dest[2][3] = src1[2][0] * src2[0][3] + src1[2][1] * src2[1][3] + src1[2][2] * src2[2][3] + src1[2][3] * src2[3][3]; 
dest[3][0] = src1[3][0] * src2[0][0] + src1[3][1] * src2[1][0] + src1[3][2] * src2[2][0] + src1[3][3] * src2[3][0]; 
dest[3][1] = src1[3][0] * src2[0][1] + src1[3][1] * src2[1][1] + src1[3][2] * src2[2][1] + src1[3][3] * src2[3][1]; 
dest[3][2] = src1[3][0] * src2[0][2] + src1[3][1] * src2[1][2] + src1[3][2] * src2[2][2] + src1[3][3] * src2[3][2]; 
dest[3][3] = src1[3][0] * src2[0][3] + src1[3][1] * src2[1][3] + src1[3][2] * src2[2][3] + src1[3][3] * src2[3][3]; 
};
于 2009-11-04T15:35:34.683 回答
3

您确定展开的代码比基于显式循环的方法更快吗？请注意，编译器通常比人类执行优化要好！

事实上，我敢打赌，编译器从编写良好的循环中自动发出 SIMD 指令的机会比从一系列“不相关”的语句中发出的机会更多……

您还可以在参数声明中指定矩阵大小。然后你可以使用普通的括号语法来访问元素，它也可以是编译器进行优化的一个很好的提示。

于 2009-11-04T15:08:59.380 回答
2

您完全展开的传统产品可能非常快。

您的矩阵太小，无法克服使用显式索引和分区代码以传统形式管理 Strassen 乘法的意外情况；您可能会对该开销的优化失去任何影响。

但是，如果您想要快速，我会使用 SIMD 指令（如果它们可用）。如果这些天的 ARM 芯片没有它们，我会感到惊讶。如果是这样，您可以在一条指令中管理行/列中的所有产品；如果 SIMD 为 8 宽，您可以在一条指令中管理2行/列乘法。设置操作数来执行该指令可能需要一些跳舞；SIMD 指令将轻松获取您的行（相邻值），但不会获取列（非连续）。计算一行/列中产品的总和可能需要一些努力。

于 2009-11-04T14:25:05.880 回答
2

这些是任意矩阵还是它们有任何对称性？如果是这样，通常可以利用这些对称性来提高性能（例如在旋转矩阵中）。

另外，我同意上面的 fortran，并且会运行一些计时测试来验证您的手动展开代码是否比优化编译器可以创建的更快。至少，您可以简化代码。

保罗

于 2009-11-04T16:09:07.980 回答