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

Осваиваю тут новую 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 кодируется так:
 for(i=0;i<N;I++)   {
    for(j=0;j<N;j++) {
       C(i,j)=0;
       for(k=0;k<N;k++) {
           C(i,j)+=A(i,k)*B(k,j); /*  **** */
       }
    }
 }
Если входные матрицы заранее транспонировать (одну или обе), то в строчке помеченной как **** можно поменять порядок индексов у A и B. Время транспонирования я в бенчмарки не считал, во-первых оно n2, а во-вторых данные обычно можно готовить уже транспонированными.

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

Для однопоточных (без попыток распараллеливания) операций были получены такие вот времена исполнения (в секундах):
Размер задачи 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
  • Для неоптимального расположения данных в памяти рост времени исполнения происходит быстрее чем n3.
  • (ожидавшееся) SSE3 (Intel) лучше чем SSE2 (VS2005)

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

Автоматическое распараллеливание, предлагаемое компилятором Intel, не дало никаких значимых результатов. Однако оба компилятора поддерживают OpenMP, а использвание его для рассматриваемой задачи тривиально:
register float sum;
#pragma omp parallel
#pragma omp for private(i, j, k, sum)
        for(i=0;i<Ne;i++){     
                for(j=0;j<psize;j++){
                        sum = 0.0;
                        for(k=0;k<psize;k++){
                                sum+=a(i,k)*b(k,j);
                        }
                        c(i,j)=sum;
                }
        }

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

Размер задачи 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 раз медленнее оптимального при одинаковом числе потоков исполнения.

Comments

о перемножении м-ц см. страницу, посвященную Core 2 Duo (IA32) (таблица 3):
http://www.thesa-store.com/products/

Тезисы - это прикольно, а как бы код пощупать

Add new comment