Почему транспонирование матрицы 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 мс
Почему это происходит?
3 ответа
Объяснение исходит от 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.:)
В качестве иллюстрации к объяснению в ответе Лучиана Григоре, вот как выглядит присутствие матричного кэша для двух случаев матриц 64x64 и 65x65 (подробности о числах см. По ссылке выше).
Цвета в анимации ниже означают следующее:
-
- не в кеше, -
- в кеше, -
- попадание в кеш, -
- просто читать из оперативной памяти, -
Промах кеша.
Корпус 64x64:
https://media.giphy.com/media/3o752kgqGGmfcti8KY/giphy.gif
Обратите внимание, что почти каждый доступ к новой строке приводит к отсутствию кэша. А теперь как выглядит обычный корпус, матрица 65х65:
https://media.giphy.com/media/d1G6dyzcFW8CvL3y/giphy.gif
Здесь вы можете видеть, что большинство обращений после начального прогрева являются попаданиями в кэш. Так работает кеш процессора в целом.
Лучиан объясняет, почему происходит такое поведение, но я подумал, что было бы неплохо показать одно возможное решение этой проблемы и в то же время немного рассказать о алгоритмах, забывающих о кеше.
Ваш алгоритм в основном делает:
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];
}
}
}
}