代码之家  ›  专栏  ›  技术社区  ›  Till

C中的快速4x4矩阵乘法

  •  16
  • Till  · 技术社区  · 15 年前

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

    目前,我使用的是一种相当标准的方法——只是展开一个小循环。

    #define O(y,x) (y + (x<<2))
    
    static inline 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算法中获益吗?

    5 回复  |  直到 15 年前
        1
  •  36
  •   Nils Pipenbrinck    15 年前

    不,Strassen或Coppersmith-Winograd算法在这里不会有什么不同。他们开始只为更大的矩阵付出代价。

    如果你的矩阵乘法真的是一个瓶颈,你可以用氖单指令多数据指令重写算法。这只对ARMV7有帮助,因为ARMV6没有这个扩展。

    我希望系数3在为您的案例编译的C代码上加速。

    编辑:您可以在arm-neon中找到一个很好的实现: http://code.google.com/p/math-neon/

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

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

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

        2
  •  20
  •   Patrick Schlüter    15 年前

    只是吹毛求疵。我想知道为什么人们仍然故意混淆他们的代码?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]; 
    };
    
        3
  •  3
  •   fortran    15 年前

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

    事实上,我敢打赌,编译器从一个编写良好的循环中自动发出simd指令的可能性比从一系列“无关”语句中发出指令的可能性要大…

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

        4
  •  2
  •   Paul    15 年前

    这些矩阵是任意的还是对称的?如果是这样的话,这些对称性常常可以用来提高性能(例如在旋转矩阵中)。

    此外,我同意上面的Fortran,并将运行一些计时测试来验证手动展开的代码比优化编译器创建的代码要快。至少,您可以简化代码。

    保罗

        5
  •  2
  •   Ira Baxter    15 年前

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

    您的矩阵太小,无法克服人们听到的用显式索引和分区代码管理传统形式的strassen乘法;您可能会失去对优化的任何影响。

    但是如果你想快点的话,如果可以的话,我会用SIMD指令。如果这几天胳膊上没有碎片,我会很惊讶的。如果是这样,您可以在一条指令中以行/列的形式管理所有产品;如果simd是8宽的,您可以管理 行/列在一条指令中相乘。将操作数设置为执行该指令可能需要一些交替;SIMD指令将很容易拾取行(相邻值),但不会拾取列(非相邻)。计算行/列中的乘积和可能需要一些努力。