Эффективная программа транспонирования матрицы в кеше?

Question

Итак, очевидный способ транспонировать матрицу - это использовать:

  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.

Answer 1

Возможно, вам понадобится четыре цикла - два для итерации по блокам, а затем еще два для выполнения транспонирования-копии одного блока.Предполагая для простоты размер блока, который делит размер матрицы, я думаю, что-то вроде этого, хотя я бы хотел нарисовать несколько картинок на обороте конвертов, чтобы быть уверенным:

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,, который использует именно эту проблему в качестве примера).Неформальное определение «не обращающего внимания на кэш» заключается в том, что вам не нужно экспериментировать с настройкой каких-либо параметров (в данном случае размера блоков), чтобы достичь хорошей / оптимальной производительности кэша.Решением в этом случае является транспонирование путем рекурсивного деления матрицы пополам и перемещения половин в их правильное положение в месте назначения.

Каким бы ни был размер кэша, эта рекурсия использует его в своих интересах.Я ожидаю, что по сравнению с вашей стратегией возникнут некоторые дополнительные издержки на управление, заключающиеся в том, чтобы использовать эксперименты с производительностью, чтобы, по сути, перейти прямо к той точке рекурсии, в которой кэш действительно срабатывает, и не идти дальше.С другой стороны, ваши эксперименты с производительностью могут дать вам ответ, который работает на вашей машине, но не на машинах ваших клиентов.

Answer 2

У меня вчера была точно такая же проблема.Я получил следующее решение:

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 но это можно легко перевести на другие языки.

Answer 3

Вместо транспонирования матрицы в памяти, почему бы не свернуть операцию транспонирования в следующую операцию, которую вы собираетесь выполнить с матрицей?

Answer 4

Стив Джессоп упомянул алгоритм транспонирования матрицы из кеша. Напомню, что я хочу поделиться возможной реализацией транспонированной матрицы кеша.

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);
        }
    }
}

Подробнее о алгоритмах кеширования можно найти здесь здесь .

Answer 5

Умножение матриц приходит на ум, но проблема с кешем там гораздо более выраженная, потому что каждый элемент читается N раз.

С помощью транспонирования матрицы вычитаем за один линейный проход, и нет способа оптимизировать это.Но вы можете одновременно обрабатывать несколько строк, чтобы написать несколько столбцов и заполнить все строки кэша.Вам понадобятся только три цикла.

Или сделайте это наоборот и читайте в столбцах при линейной записи.

Answer 6

С большой матрицей, возможно, большой разреженной матрицей, может быть идея разложить ее на более мелкие кеш-блоки (скажем, подматрицы 4x4). Вы также можете пометить подматрицы как идентификаторы, которые помогут вам в создании оптимизированных путей кода.