高效缓存矩阵转置程序？

问题描述 投票：0回答：8

因此转置矩阵的明显方法是使用：

  for( int i = 0; i < n; i++ )

    for( int j = 0; j < n; j++ )

      destination[j+i*n] = source[i+j*n];

但我想要一些能够利用局部性和缓存阻塞的东西。我正在查找它，但找不到可以执行此操作的代码，但我被告知这应该是对原始代码的非常简单的修改。有什么想法吗？

编辑：我有一个 2000x2000 矩阵，我想知道如何使用两个 

for
循环更改代码，基本上将矩阵分割成单独转置的块，例如 2x2 块或 40x40 块，然后查看哪个块尺寸是最有效的。

Edit2：矩阵按列主要顺序存储，也就是说对于矩阵

a1 a2    
a3 a4

存储为 

a1 a3 a2 a4
。

algorithm caching matrix
8个回答
49
投票

您可能需要四个循环 - 两个循环遍历块，然后另外两个循环执行单个块的转置复制。为了简单起见，假设块大小除以矩阵的大小，我认为是这样的，尽管我想在信封背面画一些图片以确保：

for (int i = 0; i < n; i += blocksize) {
    for (int j = 0; j < n; j += blocksize) {
        // transpose the block beginning at [i,j]
        for (int k = i; k < i + blocksize; ++k) {
            for (int l = j; l < j + blocksize; ++l) {
                dst[k + l*n] = src[l + k*n];
            }
        }
    }
}

一个重要的进一步见解是，实际上有一个缓存忽略算法（请参阅 http://en.wikipedia.org/wiki/Cache-oblivious_algorithm，它使用这个确切的问题作为示例）。 “缓存遗忘”的非正式定义是，您无需尝试调整任何参数（在本例中为块大小）即可达到良好/最佳的缓存性能。这种情况下的解决方案是通过将矩阵递归地分成两半来进行转置，然后将两半转置到目标中的正确位置。

无论缓存实际大小是多少，此递归都会利用它。我预计与您的策略相比，会有一些额外的管理开销，即使用性能实验实际上直接跳到递归中缓存真正启动的点，并且不再继续。另一方面，您的性能实验可能会给您一个在您的机器上有效但在客户的机器上无效的答案。

15
投票

我昨天遇到了完全相同的问题。 我最终得到了这个解决方案：

void transpose(double *dst, const double *src, size_t n, size_t p) noexcept {
    THROWS();
    size_t block = 32;
    for (size_t i = 0; i < n; i += block) {
        for(size_t j = 0; j < p; ++j) {
            for(size_t b = 0; b < block && i + b < n; ++b) {
                dst[j*n + i + b] = src[(i + b)*p + j];
            }
        }
    }
}

这比我机器上明显的解决方案快 4 倍。

此解决方案处理尺寸不是块大小的倍数的矩形矩阵。

如果 dst 和 src 是相同的方阵，则应该使用就地函数：

void transpose(double*m,size_t n)noexcept{
    size_t block=0,size=8;
    for(block=0;block+size-1<n;block+=size){
        for(size_t i=block;i<block+size;++i){
            for(size_t j=i+1;j<block+size;++j){
                std::swap(m[i*n+j],m[j*n+i]);}}
        for(size_t i=block+size;i<n;++i){
            for(size_t j=block;j<block+size;++j){
                std::swap(m[i*n+j],m[j*n+i]);}}}
    for(size_t i=block;i<n;++i){
        for(size_t j=i+1;j<n;++j){
            std::swap(m[i*n+j],m[j*n+i]);}}}

我使用了 C++11，但这可以很容易地翻译成其他语言。

9
投票

与其在内存中转置矩阵，为什么不将转置操作折叠到要对矩阵执行的下一个操作中？

7
投票

Steve Jessop 提到了一种缓存不经意矩阵转置算法。 作为记录，我想分享缓存不经意矩阵转置的可能实现。

public class Matrix {
    protected double data[];
    protected int rows, columns;

    public Matrix(int rows, int columns) {
        this.rows = rows;
        this.columns = columns;
        this.data = new double[rows * columns];
    }

    public Matrix transpose() {
        Matrix C = new Matrix(columns, rows);
        cachetranspose(0, rows, 0, columns, C);
        return C;
    }

    public void cachetranspose(int rb, int re, int cb, int ce, Matrix T) {
        int r = re - rb, c = ce - cb;
        if (r <= 16 && c <= 16) {
            for (int i = rb; i < re; i++) {
                for (int j = cb; j < ce; j++) {
                    T.data[j * rows + i] = data[i * columns + j];
                }
            }
        } else if (r >= c) {
            cachetranspose(rb, rb + (r / 2), cb, ce, T);
            cachetranspose(rb + (r / 2), re, cb, ce, T);
        } else {
            cachetranspose(rb, re, cb, cb + (c / 2), T);
            cachetranspose(rb, re, cb + (c / 2), ce, T);
        }
    }
}

有关缓存遗忘算法的更多详细信息可以在此处找到。

3
投票
我想到了

矩阵乘法，但是缓存问题更加明显，因为每个元素都被读取N次。

使用矩阵转置，您正在单次线性传递中进行读取，并且无法对其进行优化。但是您可以同时处理多行，以便写入多列，从而填充完整的缓存行。您只需要三个循环。

或者反过来，在线性书写的同时分栏阅读。

0
投票

对于一个大矩阵，可能是一个大稀疏矩阵，将其分解为更小的缓存友好块（例如 4x4 子矩阵）可能是一个想法。您还可以将子矩阵标记为标识，这将帮助您创建优化的代码路径。

0
投票

以下是用于缓存友好的并行矩阵转置的示例 C++ 代码：

constexpr auto range(unsigned max) 
    { return std::views::iota(0u, max); }
constexpr auto range(unsigned start, unsigned max) 
    { return std::views::iota(start, max); }

constexpr auto SZ = 1000u;
std::array<std::array<int, SZ>, SZ> mx;

void testtranspose(void)
{
    for (int i = 0; i < SZ; i++) 
        for (int j = 0; j < SZ; j++)
            mx[i][j] = j;

    // ***************
    //       cache friendly in place parallel transposition
    constexpr unsigned block_size = 10; assert(SZ % block_size == 0);
    std::for_each(std::execution::par_unseq,
        range(SZ / block_size).begin(), range(SZ / block_size).end(), 
        [&](const unsigned iib)
        {   //cache friendly transposition
            const unsigned ii = iib * block_size;
            //transpose the diagonal block
            for (unsigned i = ii; i < ii + block_size; i++)
                for (unsigned j = i + 1; j < ii + block_size; j++)
                    std::swap(mx[i][j], mx[j][i]);
            //transpose the rest of blocks
            for (unsigned jj = ii + block_size; jj < SZ; jj += block_size)
                for (unsigned i = ii; i < ii + block_size; i++)
                    for (unsigned j = jj; j < jj + block_size; j++)
                        std::swap(mx[i][j], mx[j][i]);
        });
    // ***************

    for (int i = 0; i < SZ; i++)
        for (int j = 0; j < SZ; j++)
            assert(mx[i][j] == i);
}
0
投票

在我的机器上，这段代码达到了 memcpy 速度的 85%（单线程 ofc）： https://godbolt.org/z/9Y3GzE1Tx

template <typename T>
void transpose2(const T* in, T* out, size_t nrow, size_t ncol) {
    static constexpr size_t kTileRowCount = 4;
    static constexpr size_t kTileColCount = 16;

    const auto transpose_block = [=](size_t col, size_t row) {
        double scratch[kTileRowCount * kTileColCount];

        for (size_t icol = 0; icol < kTileColCount; ++icol) {
            for (size_t irow = 0; irow < kTileRowCount; ++irow) {
                scratch[irow + icol * kTileRowCount] =
                    in[(row + irow) + (col + icol) * nrow];
            }
        }

        for (size_t irow = 0; irow < kTileRowCount; ++irow) {
            for (size_t icol = 0; icol < kTileColCount; ++icol) {
                out[(col + icol) + (row + irow) * ncol] =
                    scratch[irow + icol * kTileRowCount];
            }
        }
    };

    for (size_t iblockcol = 0; iblockcol < ncol; iblockcol += kTileColCount) {
        for (size_t iblockrow = 0; iblockrow < nrow; iblockrow += kTileRowCount) {
            transpose_block(iblockcol, iblockrow);
        }
    }
}
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.