Умножение матриц, серия 2: MKL против компилятора, single/double и int

Продолжаем умножать матрицы. Для начала смоделируем sgemm/dgemm: C=alpha*A*B+beta*C

Нас интересует, естественно, самый быстрый способ из изученных ранее, а вопрос заключается в разнице в скорости между float и double и разницы в скорости между простым кодом, написанным вручную, и библиотечной реализацией.

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

Все тесты делались на Dual Opteron 275 под управлением Windows XP x64. Размер кэша составляет 1 мегабайт на ядро, памяти достаточно.

Ручное написание кода

Перемножающий код практически тот же самый, добавляется половинка строчки: C(i,j)=alpha*sum+beta*D(i,j).

Самодельный код компилировался Intel C++ Compiler 9.1.028 с включенной генерацией SSE3-кода и без распараллеливания.

В таблице времена исполнения в секундах для разных типов данных.

Тип данных	размер задачи
	256	512	1024	2048	4096
	SSE3
int (32 bit)	0.031	0.25	2.25	20.22	152.4
float (32 bit)	0.015	0.078	1.89	14.78	123.1
double (64 bit)	0.015	0.391	3.55	28.0	225
	Без SSE
float (32 bit)	0.047	0.281	2.58	20.2	163.9
double (64 bit)	0.031	0.453	4.03	31.55	257.6

SSE3 рулит, умножение single float даже быстрее, чем целых. Хочется еще сделать вывод о том, что все лимитировано скоростью памяти, однако полученные для библиотечной реализации результаты не позволяют.

Intel Math Kernel Library

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

Мгновенно выяснилось, что задачи с размером 128-512 бенчмаркать нет смысла, не хватает разрешения таймера Windows. Столь же мгновенно выяснилось, что транспонирование не имеет значения, отклонения в пределах ошибки таймера. Для задач размером 1024-4096 времена исполнения (в секундах) приведены в таблице:

Тип данных	размер задачи
	1024	2048	4096
single (32 bit)	0.375	2.89	23.25
double (64 bit)	0.688	5.56	49.375

Задача 1024x1024 single precision - это 8 мегабайт активных данных и еще 4 мегабайта матрицы C которая проходится один раз. Другими словами, в кэш помещается 1/8 активных данных. Для задач бОльшего размера достаточность кэша меньше в 4 и 16 раз соответственно.

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

Параллельная версия

MKL позволяет использовать распараллеливание через OpenMP без перекомпиляции, достаточно установить переменную окружения OMP_NUM_THREADS в положительное значение.

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

Количество потоков	размер задачи
	1024	2048	4096
	single precision
2 threads	0.187	1.453	11.797
4 threads	0.094	0.734	5.95
	double precision
2 threads	0.344	2.891	25.609
4 threads	0.187	1.438	12.609

Распараллеливание происходит линейно до 4 процессоров, что даже лучше чем очень хороший результат компилятора Intel с OpenMP-директивами. Хочется еще напомнить, что лучший результат для single precision с OpenMP, SSE3-оптимизацией и так далее, полученный написанием простого кода в 8 (прописью: восемь) раз медленнее, чем вариант MKL: 6.2 cекунды против 0.73 для задачи 2048x2048.

Выводы

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