矩阵乘法:矩阵大小差异小,时序差异大

Question

我有一个矩阵乘法代码,如下所示:

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

这里,矩阵的大小由表示dimension.现在,如果矩阵的大小是2000,运行这段代码需要147秒,而如果矩阵的大小是2048,则需要447秒.所以虽然差别没有.乘法是(2048*2048*2048)/(2000*2000*2000)= 1.073,时间上的差异是447/147 = 3.有人可以解释为什么会发生这种情况吗？我预计它会线性扩展,但这不会发生.我不是要尝试制作最快的矩阵乘法代码,只是试图理解它为什么会发生.

规格:AMD Opteron双核节点(2.2GHz),2G RAM,gcc v 4.5.0

程序编译为 gcc -O3 simple.c

我也在英特尔的icc编译器上运行了这个,并看到了类似的结果.

编辑:

正如评论/答案中所建议的那样,我运行了维度= 2060的代码,需要145秒.

继承完整的计划:

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

/* change dimension size as needed */
const int dimension = 2048;
struct timeval tv; 

double timestamp()
{
        double t;
        gettimeofday(&tv, NULL);
        t = tv.tv_sec + (tv.tv_usec/1000000.0);
        return t;
}

int main(int argc, char *argv[])
{
        int i, j, k;
        double *A, *B, *C, start, end;

        A = (double*)malloc(dimension*dimension*sizeof(double));
        B = (double*)malloc(dimension*dimension*sizeof(double));
        C = (double*)malloc(dimension*dimension*sizeof(double));

        srand(292);

        for(i = 0; i < dimension; i++)
                for(j = 0; j < dimension; j++)
                {   
                        A[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                        B[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                        C[dimension*i+j] = 0.0;
                }   

        start = timestamp();
        for(i = 0; i < dimension; i++)
                for(j = 0; j < dimension; j++)
                        for(k = 0; k < dimension; k++)
                                C[dimension*i+j] += A[dimension*i+k] *
                                        B[dimension*k+j];

        end = timestamp();
        printf("\nsecs:%f\n", end-start);

        free(A);
        free(B);
        free(C);

        return 0;
}

Answer 1

这是我疯狂的猜测:缓存

可能是你可以将2行2000 double秒放入缓存中.这比32kb的L1缓存要小得多.(同时留出房间其他必要的东西)

但是当你把它提升到2048时,它会使用整个缓存(而你会溢出一些因为你需要其他东西的空间)

假设缓存策略是LRU,将缓存溢出一小部分将导致整个行被重复刷新并重新加载到L1缓存中.

另一种可能性是由于二次幂而导致的缓存关联性.虽然我认为处理器是双向L1关联的,所以在这种情况下我认为不重要.(但无论如何我会把想法抛到那里)

可能的解释2:由于L2缓存上的超对齐而导致冲突缓存未命中.

您的B数组正在列上进行迭代.因此访问是跨越式的.您的总数据大小2k x 2k约为每个矩阵32 MB.这比你的二级缓存大得多.

当数据未完全对齐时,您将在B上具有良好的空间局部性.虽然您正在跳跃行并且每个高速缓存行只使用一个元素,但是高速缓存行保留在L2高速缓存中以便在中间循环的下一次迭代中重用.

然而,当数据完全对齐时(2048),这些跳跃将全部落在相同的"缓存方式"上并且将远远超过L2缓存关联性.因此,访问的高速缓存行B不会停留在高速缓存中以进行下一次迭代.相反,它们需要从ram中完全拉入.

Answer 2

你肯定得到我称之为缓存共振的东西.这类似于别名,但不完全相同.让我解释.

高速缓存是硬件数据结构,它提取地址的一部分并将其用作表中的索引,与软件中的数组不同.(实际上,我们将它们称为硬件中的数组.)缓存数组包含数据缓存行和标记 - 有时在数组中每个索引有一个这样的条目(直接映射),有时是几个这样的(N路集合关联).提取地址的第二部分并与存储在数组中的标记进行比较.索引和标记一起唯一地标识高速缓存行存储器地址.最后,其余的地址位标识高速缓存行中的哪些字节以及访问的大小.

通常索引和标记是简单的位域.所以内存地址看起来像

  ...Tag... | ...Index... | Offset_within_Cache_Line

(有时索引和标记是哈希值,例如,其他位的几个XOR到作为索引的中间位中.更少见,有时是索引,更少有标记,就像采用高速缓存行地址模数一样质量数.这些更复杂的指数计算试图解决共振问题,我在这里解释.所有都受到某种形式的共振,但最简单的位域提取方案在共同的访问模式上遭受共振,正如你所发现的那样.)

所以,典型的价值......有许多不同的"Opteron双核"模型,我在这里看不到任何指定你拥有的模型.随机选择一本,我在AMD的网站上看到的最新手册,AMD家庭15h型号00h-0Fh,2012年3月12日的Bios和内核开发人员指南(BKDG).

(家庭15h = Bulldozer家族,最新的高端处理器--BKDG提到双核心,虽然我不知道你所描述的产品编号.但是,无论如何,同样的共振理念适用于所有处理器,只是缓存大小和关联性等参数可能会有所不同.)

从第33页开始:

AMD系列15h处理器包含一个16 KB,4路预测的L1数据缓存和两个128位端口.这是一个直写式缓存,每个周期最多支持两个128字节的加载.它分为16个库,每个宽16个字节.[...]在一个周期内,只能从L1高速缓存的给定存储区执行一次加载.

总结一下:

• 64字节高速缓存行=>高速缓存行中的6个偏移位

• 16KB/4路=>共振为4KB.

  即地址位0-5是高速缓存行偏移.

• 16KB/64B高速缓存行=> 2 ^ 14/2 ^ 6 = 2 ^ 8 =高速缓存中的256个高速缓存行.
  (修正:我最初把它误算为128.我修复了所有依赖项.)

• 4路关联=> 256/4 =缓存阵列中的64个索引.我(英特尔)称这些"集合".

  即,您可以将缓存视为包含32个条目或集合的数组,每个条目包含4个缓存行和其标记.(这比这更复杂,但没关系).

(顺便说一下,术语"集合"和"方式"有不同的定义.)

• 在最简单的方案中有6个索引位,位6-11.

  这意味着在索引位(位6-11)中具有完全相同值的任何高速缓存行将映射到同一组高速缓存.

现在看看你的程序.

C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

循环k是最里面的循环.基类型为double,8个字节.如果dimension = 2048,即2K,则B[dimension*k+j]循环访问的连续元素将是2048*8 = 16K字节.它们都将映射到同一组L1缓存 - 它们在缓存中都具有相同的索引.这意味着,缓存中没有256个缓存行可供使用,而只有4个 - 缓存的"4向关联性".

也就是说,你可能会在这个循环中每4次迭代得到一次缓存未命中.不好.

(实际上,事情有点复杂,但上面的是一个良好的了解.上面提到的B的条目的地址是一个虚拟地址,所以有可能会略有不同的物理地址.此外,推土机有办法预测缓存,可能使用虚拟地址位,因此它不必等待虚拟到物理地址的转换.但是,在任何情况下:你的代码都有16K的"共振".L1数据缓存的共振为16K.不好.)]

如果稍微更改维度,例如更改为2048 + 1,则阵列B的地址将分布在所有缓存集中.而且你会得到更少的缓存未命中.

填充阵列是一种相当常见的优化,例如更改2048到2049,以避免这种共振.但"缓存阻塞是一个更重要的优化.http://suif.stanford.edu/papers/lam-asplos91.pdf

除了缓存线共振之外,还有其他事情在这里发生.例如,L1高速缓存具有16个库,每个库宽16个字节.当维度= 2048时,内循环中的连续B访问将始终到达同一个库.因此他们无法并行 - 如果A访问恰好进入同一家银行,您将失败.

我不认为,看着它,这和缓存共振一样大.

而且,是的,可能会出现混叠现象.例如,STLF(存储到加载转发缓冲区)可能仅使用小位域进行比较,并且获得错误匹配.

(其实,如果你仔细想想,共振缓存就像是走样,与使用位字段.共振是由多个缓存线映射同一组,不被传播arond引起的.Alisaing通过匹配基于不完整的地址造成的位.)

总的来说,我的调整建议:

1. 尝试缓存阻止而无需进一步分析.我这样说是因为缓存阻塞很容易,很可能这就是你需要做的.

2. 之后,使用VTune或OProf.或者Cachegrind.要么 ...

3. 更好的是,使用调整良好的库例程来进行矩阵乘法.

Answer 3

有几种可能的解释.一个可能的解释是Mysticial建议:耗尽有限的资源(缓存或TLB).另一种可能的可能性是伪混叠停顿,当连续的存储器访问被一些2的幂(通常是4KB)的倍数分开时,可能会发生这种错误.

您可以通过为一系列值绘制时间/维度^ 3来开始缩小工作范围.如果你已经吹掉缓存或疲惫的TLB范围,你会看到一个或多或少的平坦部分,然后是2000年到2048年之间的急剧上升,接着是另一个平坦部分.如果你看到与锯齿相关的档位,你会看到一个或多或少的平面图,在2048处有一个窄的尖峰.

当然,这具有诊断能力,但并不是决定性的.如果您想最终知道减速的来源是什么,您将需要了解性能计数器,它可以明确地回答这类问题.

Answer 4

一些答案提到了L2 Cache问题.

您可以使用缓存模拟来实际验证这一点.Valgrind的cachegrind工具可以做到这一点.

valgrind --tool=cachegrind --cache-sim=yes your_executable

设置命令行参数,使其与CPU的L2参数匹配.

使用不同的矩阵大小进行测试,您可能会看到L2未命中率的突然增加.

Answer 5

我知道这太久了,但我会咬一口.正如所说的那样(缓存问题)是导致两次幂减速的原因.但是还有另外一个问题:它太慢了.如果你看看你的计算循环.

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

最内层循环每次迭代将k改变1,这意味着您只能访问距您使用的A的最后一个元素1个双倍但整个"维度"距离B的最后一个元素两倍.这不会带来任何优势缓存B的元素

如果您将其更改为:

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];

你得到完全相同的结果(模数加法相关性误差),但它更加缓存友好(本地).我试了一下,它给了很大的改进.这可以概括为

不要按定义乘以矩阵,而是乘以行

加速示例(我更改了代码以将维度作为参数)

$ diff a.c b.c
42c42
<               C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
---
>               C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];
$ make a
cc     a.c   -o a
$ make b
cc     b.c   -o b
$ ./a 1024

secs:88.732918
$ ./b 1024

secs:12.116630

作为奖励(以及与此问题相关的原因)是这个循环不会受到前一个问题的影响.

如果你已经知道了这一切,那么我道歉!