Почему транспонирование матрицы 512x512 намного медленнее, чем транспонирование матрицы 513x513?

После проведения некоторых экспериментов с квадратными матрицами разных размеров возникла закономерность. Неизменно транспонирование матрицы размера 2^nпроисходит медленнее, чем транспонирование матрицы размера2^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 мс

Почему это происходит?

Объяснение исходит от Agner Fog в Оптимизации программного обеспечения на C ++ и сводится к тому, как данные доступны и хранятся в кеше.

Условия и подробную информацию смотрите в статье о кешировании в вики , здесь я ее сужу.

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

Для конкретного адреса памяти мы можем рассчитать, какой набор должен его зеркально отображать по формуле:

set = ( address / lineSize ) % numberOfsets

Такая формула в идеале дает равномерное распределение по наборам, потому что каждый адрес памяти с большей вероятностью будет прочитан (я сказал в идеале ).

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

Я постараюсь немного последовать примеру Агнера:

Предположим, что каждый набор имеет 4 строки, каждая из которых содержит 64 байта. Сначала мы пытаемся прочитать адрес 0x2710, который идет в комплекте 28. И тогда мы также попытаться прочитать адреса 0x2F00, 0x3700, 0x3F00и 0x4700. Все они принадлежат одному и тому же набору. Перед чтением 0x4700все строки в наборе были бы заняты. Чтение этой памяти высвобождает существующую строку в наборе, строку, которая изначально удерживала 0x2710. Проблема заключается в том, что мы читаем адреса, которые (для этого примера) 0x800разделены. Это критический шаг (опять же, для этого примера).

Критический шаг также может быть рассчитан:

criticalStride = numberOfSets * lineSize

Переменные, разнесенные criticalStrideили разделенные множеством, конкурируют за одни и те же строки кэша.

Это часть теории. Далее объяснение (также Агнер, я внимательно слежу за ним, чтобы не ошибиться):

Предположим, что матрица размером 64x64 (помните, эффекты варьируются в зависимости от кеша) с кешем 8 КБ, 4 строки в наборе * размер строки 64 байта. Каждая строка может содержать 8 элементов в матрице (64-битной int).

Критическим шагом будет 2048 байтов, что соответствует 4 строкам матрицы (которая непрерывна в памяти).

Предположим, что мы обрабатываем строку 28. Мы пытаемся взять элементы этой строки и поменять их местами с элементами из столбца 28. Первые 8 элементов строки составляют строку кэша, но они перейдут в 8 различных строки кэша в столбце 28. Помните, что критический шаг составляет 4 строки (4 последовательных элемента в столбце).

Когда в столбце достигнут элемент 16 (4 строки кэша в наборе и 4 строки друг от друга = проблема), элемент ex-0 будет удален из кэша. Когда мы дойдем до конца столбца, все предыдущие строки кэша будут потеряны и потребуется перезагрузка при доступе к следующему элементу (вся строка перезаписывается).

Имея размер, не кратный критическому шагу, портит этот идеальный сценарий катастрофы, так как мы больше не имеем дело с элементами, которые имеют критический шаг по вертикали, поэтому количество перезагрузок кэша значительно сокращается.

Еще один отказ от ответственности - я только обдумал объяснение и надеюсь, что прибил его, но я могу ошибаться. В любом случае, я жду ответа (или подтверждения) от Mysticial . :)

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

Ваш алгоритм в основном делает:

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

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

Можем ли мы сделать лучше, чем это? Да, мы можем: Общий подход к этой проблеме - алгоритмы, не обращающие внимания на кеш которые, как следует из названия, избегают зависимости от конкретных размеров кеша [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, testcode для 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- это в пределах погрешности, колебания находятся в области 100 мс; этот «эталон» не очень удобен, если мы хотим получить абсолютно точные значения])

[1] Источники для этого материала: Хорошо, если вы не можете получить лекцию от кого-то, кто работал с Лайзерсоном и соавторами по этому вопросу ... Я считаю их статьи хорошей отправной точкой. Эти алгоритмы все еще довольно редко описаны - 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];
            }
        }
    }
}

В качестве иллюстрации к объяснению в ответе Лучиана Григоре , вот как выглядит присутствие матричного кэша для двух случаев матриц 64x64 и 65x65 (подробности о числах см. По ссылке выше).

Цвета в анимации ниже означают следующее:

• - не в кеше,
• - в кеше,
• - попадание в кеш,
• - просто читать из оперативной памяти,
• Промах кеша.

Корпус 64х64:

Обратите внимание, что почти каждый доступ к новой строке приводит к отсутствию кэша. А теперь как выглядит обычный корпус, матрица 65х65:

Здесь вы можете видеть, что большинство обращений после начального прогрева являются попаданиями в кэш. Так работает кеш процессора в целом.

_{Код, сгенерировавший кадры для вышеуказанных анимаций, можно увидеть здесь .}