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

в 8:45, , рубрики: c++, Intel(R) Xeon Phi(TM), алгоритм Штрассена, Алгоритмы, Блог компании Intel, Блог компании Singularis, матричные операции, параллельное программирование, метки: , ,

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

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 1

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 2

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 3

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

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

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

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 11

где

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 12

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 13

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 14

Для замера времени использовался стандартный подход с «разогревающим» запуском и усреднением времени нескольких запусков. В качестве таймера используется функция 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(TM) 7120D, 16 Гб.

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

Наименование 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
Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 15 Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 16 Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 17
Подходит для высоко-параллельного и векторизуемого кода, последовательный код выполняется медленно Подходит для последовательного кода с большими параллельными регионами, потенциальная проблема – передача данных по PCIe Подходит для высоко-параллельного кода, эффективно выполняющегося на обеих платформах, потенциальная проблема – дисбаланс нагрузки

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

Количество потоков 1 4 8 16 60 120 180 240
Штрассен, с 114.89 29.9 15.75 8.85 3.68 3.58 4.22 8.17
MKL DGEMM, c 131.79 34.38 17.27 9 2.61 1.90 2.45 2.54

На небольшом количестве потоков алгоритм Штрассена оказался быстрее Intel MKL DGEMM. Также было замечено, что на большом количестве потоков наблюдается падение производительности (задача почти не масштабируется больше 60 потоков). Для эффективного использования ресурсов Intel Xeon Phi(TM) в многопоточном приложении, рекомендуют использовать число потоков, кратное 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( … );

Общее количество потоков Штрассен, с MKL DGEMM, с
Количество потоков в MKL DGEMM
1 2 4
60 3.68 3.89 5.19 2.61
120 3.58 3.50 3.82 1.90
240 8.17 4.23 3.59 2.54

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

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

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

Распараллеливание алгоритма Штрассена на Intel® Xeon Phi(TM) - 18

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

KMP_AFFINITY = granularity=fine,balanced

Общее количество потоков Штрассен, с MKL DGEMM, с
Количество потоков в MKL DGEMM
1 2 4
60 3.67 4.07 5.53 2.64
120 3.54 3.51 3.95 1.51
240 7.11 4.33 3.41 1.15

KMP_AFFINITY = granularity=fine,compact

Общее количество потоков Штрассен, с MKL DGEMM, с
Количество потоков в MKL DGEMM
1 2 4
60 9.29 9.96 10.27 4.58
120 6.23 6.79 6.04 2.31
240 9.32 5.21 4.21 1.15

По умолчанию, переменная 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(TM).

Автор: Intel

Источник


* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js