为什么转换512x512的矩阵要比转换513x513的矩阵慢得多？

Question

在对不同尺寸的方形矩阵进行一些实验后,出现了一种模式.转换一个大小的矩阵2^n2^n+1总是比转换一个大小的矩阵慢.对于较小的值n,差异并不重要.

然而,在512的值上会出现很大的差异.(至少对我而言)

免责声明:我知道由于元素的双重交换,函数实际上并没有转置矩阵,但它没有任何区别.

遵循代码:

#define SAMPLES 1000
#define MATSIZE 512

#include <time.h>
#include <iostream>
int mat[MATSIZE][MATSIZE];

void transpose()
{
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
   {
       int aux = mat[i][j];
       mat[i][j] = mat[j][i];
       mat[j][i] = aux;
   }
}

int main()
{
   //initialize matrix
   for ( int i = 0 ; i < MATSIZE ; i++ )
   for ( int j = 0 ; j < MATSIZE ; j++ )
       mat[i][j] = i+j;

   int t = clock();
   for ( int i = 0 ; i < SAMPLES ; i++ )
       transpose();
   int elapsed = clock() - t;

   std::cout << "Average for a matrix of " << MATSIZE << ": " << elapsed / SAMPLES;
}

改变MATSIZE让我们改变大小(呃!).我在ideone上发布了两个版本:

• 大小512 - 平均2.46毫秒 - http://ideone.com/1PV7m
• 大小513 - 平均0.75毫秒 - http://ideone.com/NShpo

在我的环境中(MSVS 2010,完全优化),差异是相似的:

• 大小512 - 平均2.19毫秒
• 大小513 - 平均0.57毫秒

为什么会这样？

Answer 1

解释来自Agner Fog在C++中优化软件,它简化了数据访问和存储在缓存中的方式.

有关术语和详细信息,请参阅有关缓存的wiki条目,我将在此处缩小范围.

缓存按组和行组织.一次只使用一组,其中可以使用其中包含的任何行.一行可以镜像的行数乘以行数给出了缓存大小.

对于特定的内存地址,我们可以使用以下公式计算哪个集合应该反映它:

set = ( address / lineSize ) % numberOfsets

这种公式理想地在整个集合中提供均匀分布,因为每个存储器地址都可能被读取(我理想地说).

很明显,可能会发生重叠.在高速缓存未命中的情况下,在高速缓存中读取存储器并替换旧值.请记住,每个集合都有许多行,其中最近最少使用的行会被新读取的内存覆盖.

我会尝试在某种程度上遵循Agner的例子:

假设每组有4行,每行包含64个字节.我们首先尝试读取0x2710集合中的地址28.然后,我们也尝试读取地址0x2F00,0x3700,0x3F00和0x4700.所有这些都属于同一套.在阅读之前0x4700,集合中的所有行都将被占用.读取该内存会驱逐集合中的现有行,即最初持有的行0x2710.问题在于我们读取的地址(对于此示例)是0x800分开的.这是关键的进步(再次,对于这个例子).

关键步幅也可以计算:

criticalStride = numberOfSets * lineSize

间隔criticalStride或多个间隔的变量争用相同的高速缓存行.

这是理论部分.接下来,解释(也是Agner,我正密切关注以避免犯错):

假设一个64x64的矩阵(记住,效果根据缓存而变化)具有8kb缓存,每组4行*64字节的行大小.每行可以容纳矩阵中的8个元素(64位int).

临界步幅为2048字节,对应于矩阵的4行(在存储器中是连续的).

假设我们正在处理第28行.我们正在尝试获取该行的元素并将它们与第28列中的元素交换.该行的前8个元素构成一个缓存行,但它们将分为8个不同的第28列中的缓存行.请记住,临界步幅相隔4行(列中的4个连续元素).

当列中到达元素16时(每组4个高速缓存行和4行间隔=故障),ex-0元素将从高速缓存中逐出.当我们到达列的末尾时,所有先前的高速缓存行都将丢失,并且在访问下一个元素时需要重新加载(整行被覆盖).

具有不是关键步幅的倍数的大小会使这个完美的灾难场景变得混乱,因为我们不再处理在垂直方向上重要的元素,因此缓存重新加载的数量会大大减少.

另一个免责声明 - 我只是解释了解释并希望我钉它,但我可能会弄错.无论如何,我正在等待Mysticial的回复(或确认).:)

Answer 2

Luchian解释了为什么会出现这种情况,但我认为向这个问题展示一种可能的解决方案并同时展示一些缓存遗忘算法是个不错的主意.

你的算法基本上做了:

for (int i = 0; i < N; i++) 
   for (int j = 0; j < N; j++) 
        A[j][i] = A[i][j];

这对现代CPU来说太可怕了.一种解决方案是了解有关缓存系统的详细信息并调整算法以避免这些问题.只要你知道那些细节就行得很好..不是特别便携.

我们能做得更好吗？是的,我们可以:解决这个问题的一般方法是缓存遗忘算法,正如名称所说,避免依赖于特定的缓存大小[1]

解决方案如下所示:

void recursiveTranspose(int i0, int i1, int j0, int j1) {
    int di = i1 - i0, dj = j1 - j0;
    const int LEAFSIZE = 32; // well ok caching still affects this one here
    if (di >= dj && di > LEAFSIZE) {
        int im = (i0 + i1) / 2;
        recursiveTranspose(i0, im, j0, j1);
        recursiveTranspose(im, i1, j0, j1);
    } else if (dj > LEAFSIZE) {
        int jm = (j0 + j1) / 2;
        recursiveTranspose(i0, i1, j0, jm);
        recursiveTranspose(i0, i1, jm, j1);
    } else {
    for (int i = i0; i < i1; i++ )
        for (int j = j0; j < j1; j++ )
            mat[j][i] = mat[i][j];
    }
}

略微复杂,但是一个简短的测试在我的古老的e8400上展示了一些非常有趣的东西,带有VS2010 x64版本,测试代码为 MATSIZE 8192

int main() {
    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);
    QueryPerformanceCounter(&start);
    recursiveTranspose(0, MATSIZE, 0, MATSIZE);
    QueryPerformanceCounter(&end);
    printf("recursive: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));

    QueryPerformanceCounter(&start);
    transpose();
    QueryPerformanceCounter(&end);
    printf("iterative: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));
    return 0;
}

results: 
recursive: 480.58ms
iterative: 3678.46ms

编辑:关于大小的影响:虽然在某种程度上仍然很明显,但它更不明显,这是因为我们将迭代解决方案用作叶子节点而不是递归到1(通常的递归算法优化).如果我们设置LEAFSIZE = 1,缓存对我没有影响[ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms- 在误差范围内,波动在100ms区域; 如果我们想要完全准确的价值观,这个"基准"并不是我太舒服的事情])

[1]这些东西的来源:如果你不能从与Leiserson合作的人那里得到一个讲座,那么我认为他们的论文是一个很好的起点.这些算法仍然很少被描述 - CLR有一个关于它们的脚注.它仍然是让人们大吃一惊的好方法.

编辑(注意:我不是发布此答案的人;我只想添加此内容):
这是上述代码的完整C++版本:

template<class InIt, class OutIt>
void transpose(InIt const input, OutIt const output,
    size_t const rows, size_t const columns,
    size_t const r1 = 0, size_t const c1 = 0,
    size_t r2 = ~(size_t) 0, size_t c2 = ~(size_t) 0,
    size_t const leaf = 0x20)
{
    if (!~c2) { c2 = columns - c1; }
    if (!~r2) { r2 = rows - r1; }
    size_t const di = r2 - r1, dj = c2 - c1;
    if (di >= dj && di > leaf)
    {
        transpose(input, output, rows, columns, r1, c1, (r1 + r2) / 2, c2);
        transpose(input, output, rows, columns, (r1 + r2) / 2, c1, r2, c2);
    }
    else if (dj > leaf)
    {
        transpose(input, output, rows, columns, r1, c1, r2, (c1 + c2) / 2);
        transpose(input, output, rows, columns, r1, (c1 + c2) / 2, r2, c2);
    }
    else
    {
        for (ptrdiff_t i1 = (ptrdiff_t) r1, i2 = (ptrdiff_t) (i1 * columns);
            i1 < (ptrdiff_t) r2; ++i1, i2 += (ptrdiff_t) columns)
        {
            for (ptrdiff_t j1 = (ptrdiff_t) c1, j2 = (ptrdiff_t) (j1 * rows);
                j1 < (ptrdiff_t) c2; ++j1, j2 += (ptrdiff_t) rows)
            {
                output[j2 + i1] = input[i2 + j1];
            }
        }
    }
}

Answer 3

作为Luchian Grigore答案中解释的说明,下面是64x64和65x65矩阵的两种情况下矩阵缓存的存在情况(有关数字的详细信息,请参见上面的链接).

以下动画中的颜色表示以下内容:

• - 不在缓存中,
• - 在缓存中,
• - 缓存命中,
• - 只需从RAM读取,
• - 缓存未命中.

64x64案例:

请注意,几乎每次访问新行都会导致缓存未命中.现在它如何寻找正常情况,一个65x65矩阵:

在这里,您可以看到初始预热后的大多数访问都是缓存命中.这是CPU缓存一般用于工作的方式.