首页 文章

C中的优化矩阵乘法

提问于
浏览 次
22

我正在尝试比较矩阵乘法的不同方法 . 第一个是常规方法：

do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            suma = 0;
            for (l = 0; l < i; l++)
                suma += MatrixA[j][l]*MatrixB[l][k];
                MatrixR[j][k] = suma;
            }
        }
    }
    c++;
} while (c<iteraciones);

第二个包括首先转置矩阵B然后按行进行乘法运算：

int f, co;
for (f = 0; f < i; f++) {
    for ( co = 0; co < i; co++) {
        MatrixB[f][co] = MatrixB[co][f];
    }
}

c = 0;
do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            suma = 0;
            for (l = 0; l < i; l++)
                suma += MatrixA[j][l]*MatrixB[k][l];
                MatrixR[j][k] = suma;
            }
        }
     }
     c++;
} while (c<iteraciones);

第二种方法应该更快，因为我们正在访问连续的内存插槽，但我没有在性能上获得显着改进 . 难道我做错了什么？

我可以发布完整的代码，但我认为不需要 .

c matrix

13 回答

  • 1

    Ulrich Drepper提供了很多关于内存效率的好主意，但他特别使用矩阵乘法作为了解内存和使用这些知识可以加快这一过程的一个例子 . 请参阅他的论文中的附录A.1，并阅读6.2.1节 . 文章中的表6.2显示，对于1000x1000矩阵，他可以从一个天真的实现时间中获得10％的运行时间 .

    当然，他的最终代码非常多毛，并且使用了许多系统特定的东西和编译时调整，但是，如果你真的需要速度，阅读那篇论文并阅读他的实现绝对是值得的 .

    回复于
  • 1

    做到这一点非常重要 . 一种对大型矩阵特别重要的优化是平铺乘法以将内容保留在缓存中 . 我曾经测量过12倍的性能差异，但是我特意选择了一个消耗了我的缓存倍数的矩阵大小（大约'97，因此缓存很小） .

    有很多关于这个主题的文献 . 一个起点是：

    http://en.wikipedia.org/wiki/Loop_tiling

    有关更深入的研究，以下参考资料，尤其是Banerjee书籍，可能会有所帮助：

    [Ban93] Banerjee，Utpal，Loop Transformations for Restructuring Compilers：the Foundations，Kluwer Academic Publishers，Norwell，MA，1993 .

    [Ban94] Banerjee，Utpal，Loop Parallelization，Kluwer Academic Publishers，Norwell，MA，1994 .

    [BGS93] Bacon，David F.，Susan L. Graham和Oliver Sharp，加利福尼亚大学伯克利分校计算机科学部高性能计算编译器转换技术报告No UCB / CSD-93-781 .

    [LRW91] Lam，Monica S.，Edward E. Rothberg和Michael E Wolf . 阻塞算法的缓存性能和优化，在第四届编程语言的建筑支持国际 Session 上，于1991年4月在加利福尼亚州圣克拉拉举行，63-74 .

    [LW91] Lam，Monica S.和Michael E Wolf . 循环变换理论和最大化并行性的算法，在IEEE并行和分布式系统中，1991,2（4）：452-471 .

    [PW86] Padua，David A.和Michael J. Wolfe，“超级计算机高级编译器优化”，ACM通讯，29（12）：1184-1201,1986 .

    [Wolfe89] Wolfe，Michael J.优化超级计算机超级编译器，麻省理工学院出版社，剑桥，MA，1989 .

    [Wolfe96] Wolfe，Michael J.，并行计算高性能编译器，Addison-Wesley，CA，1996 .

    回复于
  • 0

    ATTENTION: You have a BUG in your second implementation

    for (f = 0; f < i; f++) {
    for (co = 0; co < i; co++) {
        MatrixB[f][co] = MatrixB[co][f];
    }
}

    当你做f = 0时，c = 1

    MatrixB[0][1] = MatrixB[1][0];

    你覆盖 MatrixB[0][1] 并失去那个 Value ！当循环达到f = 1时，c = 0

    MatrixB[1][0] = MatrixB[0][1];

    复制的值与已存在的值相同 .

    回复于
  • 13

    如果矩阵不够大或者您不重复操作很多次，您将看不到明显的差异 .

    如果矩阵是，例如1,000x1,000，你将开始看到改进，但我会说，如果它低于100x100，你不应该担心它 .

    此外，任何“改进”可能是毫秒级，除非yoy要么使用非常大的矩阵要么重复操作数千次 .

    最后，如果您更换的计算机速度更快，差异会更小！

    回复于
  • 4

    您获得多大的改进将取决于：

    • 缓存的大小

    • 缓存行的大小

    • 缓存的关联度

    对于小矩阵大小和现代处理器，很有可能在第一次触摸后， MatrixA 和 MatrixB 的数据几乎完全保留在缓存中 .

    回复于
  • 24

    只需要你尝试一下（但这只会对大型矩阵产生影响）：从内部循环中的乘法逻辑中分离出你的加法逻辑，如下所示：

    for (k = 0; k < i; k++)
{
    int sums[i];//I know this size declaration is illegal in C. consider 
            //this pseudo-code.
    for (l = 0; l < i; l++)
        sums[l] = MatrixA[j][l]*MatrixB[k][l];

    int suma = 0;
    for(int s = 0; s < i; s++)
       suma += sums[s];
}

    这是因为当你写suma时，你最终会停止你的管道 . 当然，在寄存器重命名等方面，大部分工作都得到了解决，但我的数量有限理解硬件，如果我想从代码中挤出每一盎司的性能，我会这样做，因为现在你不必停止管道等待写入suma . 由于乘法比加法更昂贵，你想让机器尽可能地与它并行化，因此保存加法的档位意味着你在加法循环中等待的时间少于在乘法循环中等待的时间 .

    这只是我的逻辑 . 在该领域拥有更多知识的其他人可能不同意 .

    回复于
  • 1

    您是否可以发布一些数据，比较两种矩阵大小的方法？可能是您的期望不切实际，而您的第二个版本更快，但您还没有完成测量 .

    在测量执行时间时，不要忘记包括转置矩阵B的时间 .

    您可能想要尝试的其他事情是将代码的性能与BLAS库中的等效操作的性能进行比较 . 这可能无法直接回答您的问题，但它可以让您更好地了解您对代码的期望 .

    回复于
  • 1

    两个N * N矩阵相乘的计算复杂度为O（N ^ 3） . 如果使用可能已被MATLAB采用的O（N ^ 2.73）算法，性能将得到显着提高 . 如果安装了MATLAB，请尝试将两个1024 * 1024矩阵相乘 . 在我的计算机上，MATLAB以0.7s完成它，但像你这样的天真算法的C \ C实现需要20秒 . 如果您真的关心性能，请参阅低复杂算法 . 我听说存在O（N ^ 2.4）算法，但它需要一个非常大的矩阵，以便可以忽略其他操作 .

    回复于
  • 7

    你不应该写矩阵乘法 . 你应该依赖外部库 . 特别是你应该使用 BLAS 库中的 GEMM 例程 . GEMM通常提供以下优化

    Blocking

    高效矩阵乘法依赖于阻塞矩阵并执行几个较小的阻塞乘法 . 理想情况下，选择每个块的大小以很好地适应缓存，从而大大提高性能 .

    Tuning

    理想的块大小取决于底层的内存层次结构（缓存有多大？） . 因此，应针对每台特定计算机调整和编译库 . 除其他外，这是通过 BLAS 的 ATLAS 实现的 .

    Assembly Level Optimization

    矩阵乘法很常见，开发人员会手动优化它 . 特别是这在 GotoBLAS 中完成 .

    Heterogeneous(GPU) Computing

    Matrix Multiply具有非常大的FLOP /计算能力，使其成为在GPU上运行的理想选择 . cuBLAS 和 MAGMA 是很好的候选人 .

    简而言之，密集线性代数是一个研究得很好的主题 . 人们致力于改进这些算法 . 你应该使用他们的工作;它会让他们开心 .

    回复于
  • 1

    没那么特别但更好：

    c = 0;
do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            sum = 0; sum_ = 0;
            for (l = 0; l < i; l++) {
                MatrixB[j][k] = MatrixB[k][j];
                sum += MatrixA[j][l]*MatrixB[k][l];
                l++;
                MatrixB[j][k] = MatrixB[k][j];
                sum_ += MatrixA[j][l]*MatrixB[k][l];

                sum += sum_;
            }
            MatrixR[j][k] = sum;
        }
     }
     c++;
} while (c<iteraciones);
    回复于
  • 0

    如果您正在处理小数字，那么您提到的改进可以忽略不计 . 此外，性能将根据您运行的硬件而有所不同 . 但如果你正在处理数百万的数字，那么它将会产生影响 . 来到该程序，你可以粘贴你写的程序 .

    回复于
  • 1

    非常古老的问题，但继承了我的opengl项目的当前实现：

    typedef float matN[N][N];

inline void matN_mul(matN dest, matN src1, matN src2)
{
    unsigned int i;
    for(i = 0; i < N^2; i++)
    {
        unsigned int row = (int) i / 4, col = i % 4;
        dest[row][col] = src1[row][0] * src2[0][col] +
                         src1[row][1] * src2[1][col] +
                         ....
                         src[row][N-1] * src3[N-1][col];
    }
}

    其中N被矩阵的大小替换 . 因此，如果您乘以4x4矩阵，那么您使用：

    typedef float mat4[4][4];    

inline void mat4_mul(mat4 dest, mat4 src1, mat4 src2)
{
    unsigned int i;
    for(i = 0; i < 16; i++)
    {
        unsigned int row = (int) i / 4, col = i % 4;
        dest[row][col] = src1[row][0] * src2[0][col] +
                         src1[row][1] * src2[1][col] +
                         src1[row][2] * src2[2][col] +
                         src1[row][3] * src2[3][col];
    }
}

    这个函数主要是最小化循环，但模数可能很重要......在我的计算机上，这个函数比三重循环乘法函数快大约50％ .

    缺点：

    • 需要很多代码（例如mat3 x mat3，mat5 x mat5 ......的不同功能）

    • 不规则乘法所需的调整（例如mat3 x mat4）.....

    回复于
  • 0

    一般来说，转置B应该比天真的实现快得多，但代价是浪费另一个NxN值的内存 . 我只花了一周的时间来挖掘矩阵乘法优化，到目前为止绝对不屈不挠的赢家是这样的：

    for (int i = 0; i < N; i++)
    for (int k = 0; k < N; k++)
        for (int j = 0; j < N; j++)
            if (likely(k)) /* #define likely(x) __builtin_expect(!!(x), 1) */
                C[i][j] += A[i][k] * B[k][j];
            else
                C[i][j] = A[i][k] * B[k][j];

    这甚至比早期评论中提到的Drepper方法更好，因为无论底层CPU的缓存属性如何，它都能以最佳方式工作 . 诀窍在于重新排序循环，以便以行主顺序访问所有三个矩阵 .

    回复于

相关问题