О перемножении матриц и прочих архитектурных заморочках

Осваиваю тут новую NUMA-архитектуру (пока не скажу какую, хотя многие уже в курсе). Хочется ее адекватно сравнивать с возможностями PC-CPU, т.е. опять гонять тесты. Начал с умножения матриц только на CPU.

Если оставить в стороне продвинутые алгоритмы Штрассена и Копперсмита с Виноградом, ограничившись классическим подходом с умножением "строка на столбец", то умножение матриц вполне прямолинейно: C(i,j) += A(i,k)+B(k,j) С другой стороны, входные матрицы могут быть транспонированы (одна или обе). При этом меняется pattern доступа к памяти, за счет чего ожидается некоторая разница в производительности.

Рассматривая результаты Linpack, я ожидал получить разницу между «плохим» и «хорошим» паттерном в разы. Результат превзошел все ожидания, получилась разница на полтора порядка.

Оборудование

Все эксперименты делались на dual Opteron 275 (итого 4 ядра) под Windows XP x64, собирался и выполнялся 32-битный код (сравнение 32 и 64 бит буду делать отдельно). Памяти на машине достаточно. Компиляторам отдельно сообщалось, что оборудование умеет SSE2/SSE3.

Вычисления производились с 32-битными числами. Время измерялось системной GetTickCount(), разрешение таймера оценено как 16 миллисекунд. Впрочем, для размерностей матриц начиная с 512x512 этого разрешения более чем достаточно.

Код

Прямолинейный код умножения матриц выглядит просто и незамысловато. Для квадратных матриц размерности N (в псевдокоде a-la C) выражение C=A*B кодируется так: Если входные матрицы заранее транспонировать (одну или обе), то в строчке помеченной как **** можно поменять порядок индексов у A и B. Время транспонирования я в бенчмарки не считал, во-первых оно n², а во-вторых данные обычно можно готовить уже транспонированными.

Выполнение в один поток

Для однопоточных (без попыток распараллеливания) операций были получены такие вот времена исполнения (в секундах):

Размер задачи	C(i,j)*=A(k,i)*B(k,j)	C(i,j)*=A(i,k)*B(k,j)	C(i,j)*=A(i,k)*B(j,k)
Visual Studio 2005, SSE2, 1 thread
256x256	0.188	0.094	0.031
512x512	4.59	2.78	0.33
1024x1024	45.8	24.8	2.58
2048x2048	368	245	20
Intel C++ compiler 9.01, SSE3, 1 thread
256x256	0.172	0.093	0.016
512x512	4.98	2.62	0.23
1024x1024	44.9	24.7	1.81
2048x2048	392	291	14.6

По данным из таблички можно сделать сразу много интересных наблюдений:

• Даже когда данные помещаются в L-2кэш (256x256), скорость вычислений отличается как 6:3:1 (неоптимальное-стандартное-оптимальное расположение элементов в памяти).
• Как только в кэш перестаем помещаться (начиная с 1024x1024), мы имеем отношение времен порядка 20:10:1
• Для неоптимального расположения данных в памяти рост времени исполнения происходит быстрее чем n³.
• (ожидавшееся) SSE3 (Intel) лучше чем SSE2 (VS2005)

Параллельное исполнение

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

Получилось следущее:

Размер задачи	C(i,j)*=A(k,i)*B(k,j)	C(i,j)*=A(i,k)*B(k,j)	C(i,j)*=A(i,k)*B(j,k)
Visual Studio 2005, SSE2, OpenMP, 2 threads
1024x1024	31.1	16.1	1.31
2048x2048	262	156	10.5
Intel C++ compiler 9.01, SSE3, OpenMP, 2 threads
1024x1024	33.74	16.7	1.02
2048x2048	302	168	7.8
Visual Studio 2005, SSE2, OpenMP, 4 threads
1024x1024	25.2	12.9	0.77
2048x2048	231	120	6.9
Intel C++ compiler 9.01, SSE3, OpenMP, 4 threads
1024x1024	30.79	15.2	0.75
2048x2048	280	136	6.2

Результаты еще более любопытные:

• Код с неоптимальным расположением данных в памяти распараллеливается плохо, а с оптимальным - хорошо.
• Для оптимальной раскладки по памяти, VisualStudio дает практически удвоение производительности при удвоении числа тредов.
• Компилятор Intel линейного роста не дает, разница в быстродействии между 2-мя и 4-мя threads небольшая. Однако самый быстрый код - это 4 threads в исполнении Intel

Заключение

Все полученные цифры производительности прекрасно ложатся на очевидную идею - доступ к памяти должен быть cache-aware, полезно в деталях знать как оно все устроено.

Единственной неожиданностью была величина разницы: самое быстрое умножение матриц 2048x2048 на моей машине укладывается в 6.2 секунды, самое медленное - в 60 раз медленнее, а "классическое" - «всего лишь» в 20 раз медленнее оптимального при одинаковом числе потоков исполнения.