我用CUDA, c++, c#, Java做了一些基准测试,并使用MATLAB进行验证和矩阵生成。当我用MATLAB执行矩阵乘法时,2048x2048甚至更大的矩阵几乎立即被相乘。

             1024x1024   2048x2048   4096x4096
             ---------   ---------   ---------
CUDA C (ms)      43.11      391.05     3407.99
C++ (ms)       6137.10    64369.29   551390.93
C# (ms)       10509.00   300684.00  2527250.00
Java (ms)      9149.90    92562.28   838357.94
MATLAB (ms)      75.01      423.10     3133.90

只有CUDA是有竞争力的,但我认为至少c++会有点接近,而不是慢60倍。我也不知道如何看待c#的结果。算法与c++和Java一样,但从1024年到2048年有了巨大的飞跃。

MATLAB是如何如此快速地执行矩阵乘法的?

c++代码:

float temp = 0;
timer.start();
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * matice2[m][k];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

当前回答

取决于你的Matlab版本,我相信它可能已经在使用你的GPU了。

另一件事;Matlab可以跟踪矩阵的许多性质;无论是对角线的,还是赫尔密斯的,等等,并在此基础上专门设计算法。也许它的专门化是基于你传递给它的0矩阵,或者类似的东西?也许它正在缓存重复的函数调用,这会打乱您的计时?也许它优化了重复未使用的矩阵积?

为了防止这样的事情发生,使用一个随机数矩阵,并确保通过将结果打印到屏幕或磁盘或其他地方来强制执行。

其他回答

当做矩阵乘法时,你使用朴素乘法,它需要O(n^3)的时间。

有一个矩阵乘法算法,它需要O(n^2.4)。这意味着当n=2000时,你的算法需要的计算量是最佳算法的100倍。 你真的应该去维基百科上查看矩阵乘法的页面,以获得关于有效实现矩阵乘法的进一步信息。

它在c++中很慢,因为你没有使用多线程。本质上,如果A = B C,其中它们都是矩阵,则A的第一行可以独立于第二行计算,等等。如果A、B和C都是n × n矩阵,您可以将乘法运算速度提高一个因子n^2,如

A_ {i,j} = sum_{k} b_{i,k} c_{k,j}

如果您使用Eigen [http://eigen.tuxfamily.org/dox/GettingStarted.html],多线程是内置的,线程的数量是可调的。

取决于你的Matlab版本,我相信它可能已经在使用你的GPU了。

另一件事;Matlab可以跟踪矩阵的许多性质;无论是对角线的,还是赫尔密斯的,等等,并在此基础上专门设计算法。也许它的专门化是基于你传递给它的0矩阵,或者类似的东西?也许它正在缓存重复的函数调用,这会打乱您的计时?也许它优化了重复未使用的矩阵积?

为了防止这样的事情发生,使用一个随机数矩阵,并确保通过将结果打印到屏幕或磁盘或其他地方来强制执行。

以下是我在一台特斯拉C2070上使用MATLAB R2011a +并行计算工具箱的结果:

>> A = rand(1024); gA = gpuArray(A);
% warm up by executing the operations a couple of times, and then:
>> tic, C = A * A; toc
Elapsed time is 0.075396 seconds.
>> tic, gC = gA * gA; toc
Elapsed time is 0.008621 seconds.

MATLAB使用高度优化的矩阵乘法库,这就是为什么简单的MATLAB矩阵乘法如此之快。gpuArray版本使用MAGMA。

更新了在特斯拉K20c的机器上使用R2014a,以及新的timeit和gputimeit函数:

>> A = rand(1024); gA = gpuArray(A);
>> timeit(@()A*A)
ans =
    0.0324
>> gputimeit(@()gA*gA)
ans =
    0.0022

在拥有16个物理核和特斯拉V100的WIN64机器上使用R2018b进行更新:

>> timeit(@()A*A)
ans =
    0.0229
>> gputimeit(@()gA*gA)
ans =
   4.8019e-04

(注意:在某些时候(我忘记确切的时间)gpuArray从MAGMA切换到cuBLAS -岩浆仍然用于一些gpuArray操作)

在有32个物理核和A100 GPU的WIN64机器上使用R2022a更新:

>> timeit(@()A*A)
ans =
    0.0076
>> gputimeit(@()gA*gA)
ans =
   2.5344e-04

这种鲜明的对比不仅是由于Matlab的惊人优化(正如许多其他答案已经讨论过的那样),而且是由于您将矩阵作为一个对象来表述的方式。

看起来你把矩阵变成了列表的列表?列表的列表包含指向列表的指针,然后包含您的矩阵元素。所包含列表的位置是任意分配的。当循环遍历第一个索引(行号?)时,内存访问时间非常重要。相比之下,为什么不尝试实现矩阵作为一个单一的列表/向量使用下面的方法?

#include <vector>

struct matrix {
    matrix(int x, int y) : n_row(x), n_col(y), M(x * y) {}
    int n_row;
    int n_col;
    std::vector<double> M;
    double &operator()(int i, int j);
};

And

double &matrix::operator()(int i, int j) {
    return M[n_col * i + j];
}

应该使用相同的乘法算法,以使触发器的数量相同。(大小为n的方阵为n^3)

我要求您对它进行计时,以便结果与前面(在同一台机器上)的结果相比较。通过比较,您将准确地显示内存访问时间有多么重要!