Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM)21.12.2015 15:48

Сопроцессоры Intel Xeon Phi^™ представляют собой PCI Express устройство и имеют x86 архитектуру, обеспечивая высокую пиковую производительности — до 1,2 терафлопс (триллион операций с плавающей запятой в секунду) двойной точности на сопроцессор. Xeon Phi^™ может обеспечивать одновременную работу до 244 потоков, и это нужно учитывать при программировании для достижения максимальной эффективности.

Недавно мы вместе с компанией Intel проводили небольшое исследование эффективности реализации алгоритма Штрассена для сопроцессора Intel Xeon Phi^™. Кому интересны тонкости работы с этим устройством и просто любящих параллельное программирование, прошу под кат.

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

Алгоритм Штрассена (1969), который относится к быстрым алгоритмам умножения матриц требует:

Таким образом, алгоритм Штрассена имеет меньшую асимптотическую сложность и потенциально быстрее на больших размерах матриц. Но рекурсивная природа данного алгоритма представляет определенную сложность для эффективного распараллеливания на современных вычислительных системах и использования данных из памяти.

В данной работе рассматривается комбинированный подход, в котором при достижении порогового значения из алгоритма Штрассена вызывается стандартный алгоритм умножения матриц (в качестве стандартного алгоритма мы использовали Intel MKL DGEMM). Это связано с тем, что при маленьких размерах матриц (десятки или сотни, в зависимости от архитектуры процессора) алгоритм Штрассена начинает проигрывать стандартному алгоритму и по количеству операций, и за счет большего числа кэш-промахов. Теоретические оценки для порогового значения размера матриц для перехода на стандартный алгоритм (без учета кэширования и конвейерной обработки) дает различные оценки — 8 у Хайема и 12 у Хасс-Ледермана, на практике пороговое значение зависит от архитектуры вычислительной системы и должно оцениваться экспериментально. Для нашей системы с Intel Xeon Phi^™ 7120D наиболее эффективным оказалось значение 1024.

Результаты вычислительных экспериментов сравнивались с функцией DGEMM из библиотеки Intel MKL. Рассматривалось умножение квадратных матриц размера , где и . Под M здесь понимается пороговое значение размера матриц, после которого вызывается стандартный алгоритм. Умножение двух матриц записывается как , где A, B и C матрицы размера . Метод Штрассена для умножения матриц основан на рекурсивном делении каждой перемножаемой матрицы на 4 подматрицы и выполнении операций над ними. Требование к размеру матриц ( и ) необходимо, чтобы была возможность разбить каждую матрицу на 4 подматрицы на требуемую глубину рекурсии.

Метод Штрассена описывается следующей формулой:

где

Ниже показано, как можно выделить 4 группы операций, в каждой из них все операции могут быть выполнены параллельно.

В реализации распараллеливания нет ничего сложного. Использовалось OpenMP и механизм задач (omp task), распределение нагрузки между задачами показано на рисунке. Первые группы операций были совмещены в одну, так получается меньше точек синхронизации, так как это очень дорогая операция.

Для замера времени использовался стандартный подход с «разогревающим» запуском и усреднением времени нескольких запусков. В качестве таймера используется функция Intel MKL dsecnd.

double run_method(MatrixMultiplicationMethod mm_method, int n, double *A, double *B, double *C)
{
                int runs = 5;
 
                mm_method(n, A, B, C);
 
                double start_time = dsecnd();
 
                for (int i = 0; i < runs; ++i)
                               mm_method(n, A, B, C);
 
                double end_time = dsecnd();
 
                double elapsed_time = (end_time - start_time) / runs;
 
                return elapsed_time;
}

После реализации распараллеливания мы провели ряд тестов. Тесты проводились в нативном режиме на Intel Xeon Phi^™ 7120D, 16 Гб.

какие Xeon Phi^™ и режимы использования ещё бывают

Наименование	TDP (WATTS)	Число ядер	Частота (GHz)	Пиковая производительность (GFLOP)	Пиковая пропускная способность (GB/s)	Объем памяти (GB)
3120A	300	57	1.1	1003	240	6
3120P	300	57	1.1	1003	240	6
5110P	225	60	1.053	1011	320	8
5120D	245	60	1.053	1011	352	8
7120P	300	61	1.238	1208	352	16
7120X	300	61	1.238	1208	352	16
7120A	300	61	1.238	1208	352	16
7120D	270	61	1.238	1208	352	16

Подходит для высоко-параллельного и векторизуемого кода, последовательный код выполняется медленно	Подходит для последовательного кода с большими параллельными регионами, потенциальная проблема — передача данных по PCIe	Подходит для высоко-параллельного кода, эффективно выполняющегося на обеих платформах, потенциальная проблема — дисбаланс нагрузки

Чтобы не усложнять статью, все тесты будут проводиться на размере матриц 8192×8192 (данный размер является почти предельным по потреблению памяти для распараллеленного алгоритма Штрассена — около 10GB — и является достаточно большим, чтобы ощутить эффект меньшей асимптотической сложности алгоритма).

На небольшом количестве потоков алгоритм Штрассена оказался быстрее Intel MKL DGEMM. Также было замечено, что на большом количестве потоков наблюдается падение производительности (задача почти не масштабируется больше 60 потоков). Для эффективного использования ресурсов Intel Xeon Phi^™ в многопоточном приложении, рекомендуют использовать число потоков, кратное NP-1, где NP — количество процессоров в устройстве (в нашем случае 61). Подробнее можно почитать здесь.

Возникла мысль, что помочь дальнейшему масштабированию может использование параллелизма внутри вызываемого из Штрассена DGEMM. Для управления количеством потоков в MKL есть несколько функций: mkl_set_num_threads, mkl_set_num_threads_local, mkl_domain_set_num_threads. При попытке использовать mkl_set_num_threads мы обнаружили, что данная функция не влияет на количество потоков в MKL DGEMM, вызываемой из алгоритма Штрассена (на количество потоков в MKL за пределами Штрассена данная функция влияла). Вложенный параллелизм при этом был включен (OMP_NESTED=true).

Как выяснилось, MKL активно сопротивляется вложенному параллелизму. MKL использует переменную окружения MKL_DYNAMIC, которая по умолчанию равна true. Данная переменная позволяет MKL уменьшать количество потоков, которое задает пользователь, в частности при вызове функций MKL из параллельной области будет принудительно установлено количество потоков, равное 1. Чтобы разрешить параллелизм в MKL DGEMM, мы использовали функцию mkl_set_dynamic (0).

mkl_set_dynamic(0);    
 
omp_set_nested(1);     
 
set_num_threads(num_threads);
 
mkl_set_num_threads(mkl_num_threads);
 
strassen( … );
…
mkl_set_num_threads(num_threads * mkl_num_threads);
 
dgemm( … );

Результаты алгоритма Штрассена для 120 потоков при использовании многопоточного MKL DGEMM стали немного лучше, но серьезного выигрыша мы не получили.

Следующим нашим шагом было изучение привязки программных потоков OpenMP к физическим и логическим ядрам (binding). На Xeon Phi для управления привязкой потоков к ядрам служат переменные окружения KMP_AFFINITY и KMP_PLACE_THREADS. Хорошее описание нашли здесь.

Самым важным среди параметров KMP_AFFINITY является type, который управляет очередностью назначения потоков на ядра (см. рис. ниже).

Были использованы следующие значения KMP_AFFINITY:

KMP_AFFINITY = granularity=fine, balanced

KMP_AFFINITY = granularity=fine, compact

По умолчанию, переменная KMP_AFFINITY= granularity=core, balanced. При тестировании выяснилось, что лучшие параметры для умножения матриц KMP_AFFINITY= granularity=fine, balanced, причем на результаты MKL DGEMM данные параметр влияет не так сильно, как на алгоритм Штрассена, где по сравнению с KMP_AFFINITY= granularity=fine, compact наблюдается двукратное сокращение времени работы. Также заметим, что изменение переменной окружения KMP_AFFINITY с его значения по умолчанию (KMP_AFFINITY= granularity=core, balanced) сократило время работы MKL DGEMM с минимальных 1.90 с до 1.15 с (примерно в 1,5 раза). Результаты MKL DGEMM с compact и balanced отличаются предсказуемым образом, например при 120 потоках вариант с compact использует 30 ядер по 4 потока, а balanced — 60 по 2, за счет большего количества ядер и получается почти двукратное ускорение в варианте balanced.

Еще мы попробовали профилировать программу на Xeon Phi, как это сделать можно почитать здесь. Чтобы профилировать только алгоритм умножения мы использовали возможности VTune API.

В итоге мы не смогли догнать MKL DGEMM на максимальном количестве потоков, но немного больше узнали про программирование такого мощного вычислителя, как Intel Xeon Phi^™.

© Habrahabr.ru