优化数组转置函数答案

【问题标题】：Optimizing array transposing function优化数组转置函数
【发布时间】：2012-06-04 08:22:21
【问题描述】：

我正在处理一项家庭作业，但我的解决方案已经被困了好几个小时。我们遇到的问题是优化下面的代码，让它运行得更快，不管它变得多么混乱。我们应该使用诸如利用缓存块和循环展开之类的东西。

问题：

//transpose a dim x dim matrix into dist by swapping all i,j with j,i
void transpose(int *dst, int *src, int dim) {
    int i, j;

    for(i = 0; i < dim; i++) {
        for(j = 0; j < dim; j++) {
                dst[j*dim + i] = src[i*dim + j];
        }
    }
}

到目前为止我所拥有的：

//attempt 1
void transpose(int *dst, int *src, int dim) {
    int i, j, id, jd;

    id = 0;
    for(i = 0; i < dim; i++, id+=dim) {
        jd = 0;
        for(j = 0; j < dim; j++, jd+=dim) {
                dst[jd + i] = src[id + j];
        }
    }
}

//attempt 2
void transpose(int *dst, int *src, int dim) {
    int i, j, id;
    int *pd, *ps;
    id = 0;
    for(i = 0; i < dim; i++, id+=dim) {
        pd = dst + i;
        ps = src + id;
        for(j = 0; j < dim; j++) {
                *pd = *ps++;
                pd += dim;
        }
    }
}

一些想法，如果我错了，请纠正我：

我考虑过循环展开，但我认为这不会有帮助，因为我们不知道 NxN 矩阵是否具有素数维度。如果我对此进行检查，它将包含过多的计算，这只会减慢函数的速度。

缓存块不会很有用，因为无论如何，我们将线性访问一个数组 (1,2,3,4)，而另一个我们将在 N 的跳转中访问。虽然我们可以得到函数滥用缓存并更快地访问 src 块，将它们放入 dst 矩阵仍然需要很长时间。

我也尝试过使用指针而不是数组访问器，但我认为这实际上不会以任何方式加速程序。

任何帮助将不胜感激。

谢谢

【问题讨论】：

标签： c caching optimization loops matrix

【解决方案1】：

展开对于大型矩阵很有用。
如果矩阵大小不是展开次数的倍数，您将需要一些代码来处理多余的元素。但这将超出最关键的循环，因此对于大型矩阵来说这是值得的。

关于访问的方向 - 线性读取和写入 N 的跳跃可能会更好，而不是相反。这是因为读取操作会阻塞 CPU，而写入操作不会（达到限制）。

其他建议：
1. 可以使用并行化吗？ OpenMP 可以提供帮助（但如果您希望提供单 CPU 性能，那就不好了）。
2. 拆解函数并阅读，重点关注最里面的循环。您可能会发现在 C 代码中不会注意到的东西。
3. 使用递减计数器（停止在 0）可能比递增计数器更有效。
4. 编译器必须假定src 和dst 可能有别名（指向相同或重叠的内存），这限制了它的优化选项。如果你能以某种方式告诉编译器它们不能重叠，那可能会有很大的帮助。但是，我不确定该怎么做（也许使用 restrict 限定符）。

【讨论】：

【解决方案2】：

只是一个想法如何实现展开：

void transpose(int *dst, int *src, int dim) {
    int i, j;
    const int dim1 = (dim / 4) * 4;

    for(i = 0; i < dim; i++) {
        for(j = 0; j < dim1; j+=4) {
                dst[j*dim + i]     = src[i*dim + j];
                dst[(j+1)*dim + i] = src[i*dim + (j+1)];
                dst[(j+2)*dim + i] = src[i*dim + (j+2)];
                dst[(j+3)*dim + i] = src[i*dim + (j+3)];
        }
        for( ; j < dim; j++) {
                dst[j*dim + i] = src[i*dim + j];
        }
        __builtin_prefetch (&src[(i+1)*dim], 0, 1);
    }
}

当然，您应该从内部循环中删除计数（例如 i*dim），就像您在尝试中所做的那样。

缓存预取可用于源矩阵。

【讨论】：

【解决方案3】：

缓存阻塞很有用。举个例子，假设我们有一个 64 字节的缓存线大小（这就是 x86 现在使用的）。因此，对于一个足够大的矩阵，它大于缓存大小，那么如果我们转置一个 16x16 块（因为 sizeof(int) == 4，因此 16 个整数适合缓存行，假设矩阵在缓存行边界上对齐) 我们需要从内存中加载 32 个（源矩阵中的 16 个，目标矩阵中的 16 个）缓存行并存储另外 16 个行（即使存储不是顺序的）。相反，在没有缓存阻塞的情况下，转置等效的 16*16 元素需要我们从源矩阵加载 16 条缓存线，但要加载 16*16=256 条缓存线，然后为目标矩阵存储。

【讨论】：

这是要走的路。 “缓存遗忘矩阵转置”是谷歌短语。注意：通过使用 2*2 块的 16*16 缓存行，可以填充 4096 字节，这是（大多数）x86 机器上的内存页。
是的！！！根据我的经验，优化内存访问可以带来几倍的改进。
那么我该如何在这段代码中实现缓存阻塞呢？我只是在想如果给我一个 100x100 的数组，我认为这不会受益。可能我完全听不懂你在说什么。
@Glen：你有外部循环，你在块上循环，然后内部循环循环在块中的每个元素上。所以你有4个嵌套循环。对于 100x100 阵列，好处可能不会那么大；尝试并测量！

【解决方案4】：

您可能知道这一点，但register int（您告诉编译器将其放入寄存器会很聪明）。并且使 int 的 unsigned 可能会使事情进展得更快一些。

【讨论】：

注册关键字在那里并没有真正的帮助。问题是面向缓存/内存的微优化寄存器使用无济于事。
int 和 unsigned int 指令级别相同。

【解决方案5】：

凌乱不是问题，所以：我会在每个矩阵中添加一个transposed 标志。该标志指示矩阵的存储数据数组是以正常顺序还是转置顺序解释。

除了每个矩阵参数之外，所有矩阵运算都应该接收这些新标志。在每个操作中实现所有可能的标志组合的代码。也许宏可以在这里省去多余的写法。

在这个新实现中，矩阵转置只是切换标志：转置操作所需的空间和时间是恒定的。

【讨论】：