Умножение матриц: небольшая разница в размере матрицы, большая разница во времени
У меня есть код умножения матрицы, который выглядит следующим образом:
for(i = 0; i < dimension; i++)
for(j = 0; j < dimension; j++)
for(k = 0; k < dimension; k++)
C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
Здесь размер матрицы представлен dimension, Теперь, если размер матриц равен 2000, запуск этого фрагмента кода занимает 147 секунд, тогда как если размер матриц равен 2048, то это занимает 447 секунд. Так что пока разницы нет. умножения (2048*2048*2048)/(2000*2000*2000) = 1,073, разница во времени составляет 447/147 = 3. Может кто-нибудь объяснить, почему это происходит? Я ожидал, что он будет линейно масштабироваться, чего не происходит. Я не пытаюсь сделать самый быстрый матричный код, просто пытаюсь понять, почему это происходит.
Технические характеристики: двухъядерный узел AMD Opteron (2,2 ГГц), 2 ГБ ОЗУ, gcc v 4.5.0
Программа составлена как gcc -O3 simple.c
Я также запустил это на ICC-компиляторе Intel и увидел похожие результаты.
РЕДАКТИРОВАТЬ:
Как предложено в комментариях / ответах, я запустил код с измерением =2060, и это занимает 145 секунд.
Вот полная программа:
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
/* change dimension size as needed */
const int dimension = 2048;
struct timeval tv;
double timestamp()
{
double t;
gettimeofday(&tv, NULL);
t = tv.tv_sec + (tv.tv_usec/1000000.0);
return t;
}
int main(int argc, char *argv[])
{
int i, j, k;
double *A, *B, *C, start, end;
A = (double*)malloc(dimension*dimension*sizeof(double));
B = (double*)malloc(dimension*dimension*sizeof(double));
C = (double*)malloc(dimension*dimension*sizeof(double));
srand(292);
for(i = 0; i < dimension; i++)
for(j = 0; j < dimension; j++)
{
A[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
B[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
C[dimension*i+j] = 0.0;
}
start = timestamp();
for(i = 0; i < dimension; i++)
for(j = 0; j < dimension; j++)
for(k = 0; k < dimension; k++)
C[dimension*i+j] += A[dimension*i+k] *
B[dimension*k+j];
end = timestamp();
printf("\nsecs:%f\n", end-start);
free(A);
free(B);
free(C);
return 0;
}
5 ответов
Вот мое дикое предположение: кеш
Это может быть, что вы можете соответствовать 2 рядам 2000 doubleс в кеш. Который чуть меньше, чем кеш L1 32 КБ. (оставляя в комнате другие необходимые вещи)
Но когда вы поднимаете его до 2048 года, он использует весь кеш (и вы проливаете часть, потому что вам нужно место для других вещей)
Предполагая, что политика кэширования - это LRU, разлив крошечный крошечный бит приведет к многократному сбросу всей строки и ее повторной загрузке в кэш L1.
Другая возможность - ассоциативность кэша из-за степени двойки. Хотя я думаю, что процессор является двухсторонним L1-ассоциативным, поэтому я не думаю, что это имеет значение в этом случае. (но я все равно брошу эту идею)
Возможное объяснение 2: Отсутствует конфликт кэша из-за супер-выравнивания в кэше L2.
Ваш B массив итерируется по столбцу. Таким образом, доступ облегчен. Ваш общий размер данных 2k x 2k что составляет около 32 МБ на матрицу. Это намного больше, чем ваш кэш L2.
Когда данные не выровнены идеально, вы будете иметь приличную пространственную локальность на B. Несмотря на то, что вы перемещаетесь по строкам и используете только один элемент на строку кэширования, линия кэширования остается в кэше L2 для повторного использования на следующей итерации среднего цикла.
Однако, когда данные выровнены идеально (2048), все эти скачки окажутся на одном и том же "пути кеширования" и намного превысят вашу ассоциативность кеша L2. Таким образом, доступ к строкам кэша B не будет оставаться в кэше для следующей итерации. Вместо этого их нужно будет тянуть полностью от барана.
Вы определенно получаете то, что я называю кеш- резонансом. Это похоже на псевдонимы, но не совсем то же самое. Позволь мне объяснить.
Кэши - это аппаратные структуры данных, которые извлекают одну часть адреса и используют ее в качестве индекса в таблице, в отличие от программного массива. (На самом деле мы называем их массивами в аппаратном обеспечении.) Массив кэша содержит строки данных кэша и теги - иногда одна такая запись на индекс в массиве (прямое отображение), иногда несколько таких (ассоциативность N-way set). Вторая часть адреса извлекается и сравнивается с тегом, хранящимся в массиве. Вместе индекс и тег однозначно идентифицируют адрес памяти строки кэша. Наконец, оставшиеся биты адреса идентифицируют, какие байты в строке кэша адресованы, а также размер доступа.
Обычно индекс и тег являются простыми битовыми полями. Таким образом, адрес памяти выглядит
...Tag... | ...Index... | Offset_within_Cache_Line
(Иногда индекс и тег - это хэши, например, несколько XOR других битов в битах среднего диапазона, которые являются индексом. Гораздо реже, иногда индекс, и реже тег, это такие вещи, как получение адреса строки кэша по модулю a простое число. Эти более сложные вычисления индекса являются попытками борьбы с проблемой резонанса, которую я объясняю здесь. Все страдают той или иной формой резонанса, но самые простые схемы извлечения битового поля страдают резонансом на общих схемах доступа, как вы обнаружили.)
Итак, типичные значения... есть много разных моделей "Opteron Dual Core", и я не вижу здесь ничего, что бы указывало, какая у вас есть. Наугад выбирая одно, самое последнее руководство, которое я вижу на веб-сайте AMD, в Руководстве разработчика по Bios и Kernel (BKDG) для 15-часовых моделей семейства AMD 00h-0Fh, 12 марта 2012 г.
(Семейство 15h = семейство Bulldozer, самый последний высокопроизводительный процессор - BKDG упоминает о двухъядерном процессоре, хотя я не знаю номер продукта, который в точности соответствует тому, что вы описываете. Но, в любом случае, одна и та же идея резонанса применима ко всем процессорам, просто такие параметры, как размер кэша и ассоциативность, могут немного отличаться.)
Из с.33:
Процессор AMD Family 15h содержит 16-килобайтный 4-сторонний прогнозируемый кэш данных L1 с двумя 128-битными портами. Это сквозной кэш, который поддерживает до двух загрузок по 128 байт за цикл. Он разделен на 16 банков, каждый по 16 байт. [...] Только одна загрузка может быть выполнена из данного банка кэша L1 за один цикл.
Подводить итоги:
64-байтовая строка кэша => 6 смещенных бит в строке кэша
16KB/4-way => резонанс 4KB.
Т.е. биты адреса 0-5 являются смещением строки кэша.
16 КБ / 64B строк кэша => 2^14/2^6 = 2^8=256 строк кэша в кэше.
(Исправление: я изначально просчитал это как 128. я исправил все зависимости.)4 way ассоциативно => 256/4 = 64 индекса в массиве кеша. Я (Intel) называю эти "наборы".
то есть вы можете считать кеш массивом из 32 записей или наборов, каждая из которых содержит 4 строки кеша и свои теги. (Это сложнее, чем это, но это нормально).
(Кстати, термины "набор" и "способ" имеют различные определения.)
в самой простой схеме есть 6 битов индекса, биты 6-11.
Это означает, что любые строки кэша, которые имеют абсолютно одинаковые значения в битах индекса, биты 6-11, будут отображаться в один и тот же набор кэша.
Теперь посмотрите на вашу программу.
C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
Петля k является самой внутренней петлей. Базовый тип - двойной, 8 байтов. Если размерность =2048, то есть 2K, то последовательные элементы B[dimension*k+j] доступ к циклу будет на расстоянии 2048 * 8 = 16 Кбайт. Все они будут отображаться в один и тот же набор кеша L1 - у них у всех будет одинаковый индекс в кеше. Это означает, что вместо 256 кеш-строк в кеше, доступных для использования, будет только 4 - "4-сторонняя ассоциативность" кеша.
Т.е. вы, вероятно, будете получать пропуск кеша каждые 4 итерации в этом цикле. Нехорошо.
(На самом деле, все немного сложнее. Но вышеизложенное является хорошим первым пониманием. Адреса записей B, упомянутых выше, являются виртуальными адресами. Так что могут быть немного другие физические адреса. Кроме того, в Bulldozer есть способ прогнозирующего кэша, возможно, с использованием битов виртуальных адресов, чтобы не пришлось ждать преобразования виртуальных адресов в физические. Но, в любом случае: ваш код имеет "резонанс" 16K. Кэш данных L1 имеет резонанс 16K. Не хорошо.)]
Если вы немного измените размерность, например, на 2048+1, то адреса массива B будут распределены по всем наборам кэша. И вы получите значительно меньше промахов кеша.
Это довольно распространенная оптимизация для заполнения ваших массивов, например, чтобы изменить 2048 на 2049, чтобы избежать этого резонанса. Но "блокировка кеша - еще более важная оптимизация. http://suif.stanford.edu/papers/lam-asplos91.pdf
Помимо резонанса строки кэша, здесь происходят и другие вещи. Например, кэш L1 имеет 16 банков, каждый по 16 байт. При измерении = 2048 последовательные доступы B во внутреннем цикле всегда будут идти в один и тот же банк. Таким образом, они не могут идти параллельно - и если доступ А случится в одном и том же банке, вы проиграете.
Я не думаю, глядя на это, что это так же сильно, как резонанс кеша.
И, да, возможно, может появиться псевдоним. Например, буферы STLF (Store To Load Forwarding) могут сравнивать только с использованием небольшого битового поля и получать ложные совпадения.
(На самом деле, если задуматься, резонанс в кеше подобен псевдониму, связанному с использованием битовых полей. Резонанс вызван тем, что несколько строк кэша отображают один и тот же набор, а не разбросаны по областям. Алисинг вызван сопоставлением на основе неполного адреса бит.)
В целом моя рекомендация по тюнингу:
Попробуйте заблокировать кеш без дальнейшего анализа. Я говорю это потому, что блокировка кэша проста, и вполне вероятно, что это все, что вам нужно сделать.
После этого используйте VTune или OProf. Или Кашегринд. Или же...
Еще лучше, используйте хорошо настроенную библиотечную процедуру для умножения матрицы.
Есть несколько возможных объяснений. Одним из вероятных объяснений является то, что предлагает Mysticial: исчерпание ограниченного ресурса (кеш или TLB). Другая вероятная возможность - ложное блокирование псевдонимов, которое может произойти, когда последовательные обращения к памяти разделены кратным некоторой степени двойки (часто 4 КБ).
Вы можете начать сужать то, что работает, нанося график времени / измерения ^3 для диапазона значений. Если вы взорвали кэш или исчерпали TLB, вы увидите более или менее плоский участок, за которым следует резкий рост между 2000 и 2048 годами, а затем еще один плоский участок. Если вы видите киоски, связанные с псевдонимами, вы увидите более-менее плоский график с узким всплеском вверх в 2048 году.
Конечно, это имеет диагностическую силу, но это не является окончательным. Если вы хотите точно знать, что является источником замедления, вам нужно узнать о счетчиках производительности, которые могут дать окончательный ответ на этот вопрос.
Я знаю, что это слишком долго, но я откушу. Это (как уже было сказано) проблема с кешем, которая вызывает замедление примерно в два раза. Но есть еще одна проблема: это слишком медленно. Если вы посмотрите на свой цикл вычислений.
for(i = 0; i < dimension; i++)
for(j = 0; j < dimension; j++)
for(k = 0; k < dimension; k++)
C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
Самый внутренний цикл меняет k на 1 на каждой итерации, что означает, что вы получаете доступ только на 1 двойное расстояние от последнего использованного вами элемента A, но целое "измерение" удваивается от последнего элемента B. Это не дает никаких преимуществ кеширование элементов Б.
Если вы измените это на:
for(i = 0; i < dimension; i++)
for(j = 0; j < dimension; j++)
for(k = 0; k < dimension; k++)
C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];
Вы получаете точно такие же результаты (по модулю ошибок двойного сложения), но они намного более удобны для кэша (локально). Я попробовал это, и это дает существенные улучшения. Это можно суммировать как
Не умножайте матрицы по определению, а скорее по строкам
Пример ускорения (я изменил ваш код, чтобы взять измерение в качестве аргумента)
$ diff a.c b.c
42c42
< C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
---
> C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];
$ make a
cc a.c -o a
$ make b
cc b.c -o b
$ ./a 1024
secs:88.732918
$ ./b 1024
secs:12.116630
В качестве бонуса (и что делает это связанным с этим вопросом) является то, что этот цикл не страдает от предыдущей проблемы.
Если вы уже знали все это, то я прошу прощения!
В нескольких ответах упоминаются проблемы с кэшем L2.
Вы можете проверить это с помощью имитации кэша. Инструмент Valgrind cachegrind может сделать это.
valgrind --tool=cachegrind --cache-sim=yes your_executable
Установите параметры командной строки, чтобы они соответствовали параметрам L2 вашего ЦП.
Протестируйте его с матрицами разных размеров, и вы, вероятно, увидите внезапное увеличение коэффициента промаха L2.