Программирование экзафлопсных систем

Известно, что иногда распараллеливание программ приводит не к ускорению, а к замедлению. Известны также случаи, когда суперкомпьютеры используются менее чем на 1% своей производительности. Не угрожает ли это и первым экзафлопсным суперкомпьютерам, появление которых ожидается уже через 6–7 лет? Различные задачи максимально эффективны на различных вычислительных архитектурах [1], следовательно, на фиксированной архитектуре у каждого приложения будет свой потолок производительности. Все ли сделано, чтобы приблизиться к этому потолку? Попробуем ответить на этот вопрос на примере анализа обычной программы перемножения матриц, различные варианты которой сегодня «зашиты» во множество стандартных библиотек LAPACK, ATLAS, MKL, OpenBLAS и т. п., обрузующих основу почти любых приложений моделирования, анализа, математической физики (задача Дирихле), а также управления базами данных, предназначенных и для выполнения на будущем экзафлопсном суперкомпьютере:

do i = 1, n; do j = 1, n; do k = 1, n

{X( i,j ) = X( i,j )+A( i,k )*B( k,j )}

Запустим ее без оптимизации, в том числе и предоставляемой компилятором, на компьютере с двухъядерным или четырехъядерным процессором сначала последовательно, а потом с использованием стандартных средств параллельного программирования. Для размерности матрицы n = 30 явно видно преимущество параллельной программы по сравнению с последовательной. Дальнейшее увеличение размерности матрицы также приводит к ускорению, но при размере матриц n = 600 преимущество параллельной программы начинает уменьшаться, а при n > 1000 сходит на нет. А ведь на экзафлопсном суперкомпьютере будут решаться задачи с большими объемами данных — значит, возможности параллелизма в рамках одного узла весьма ограниченны? Странно делать недешевый экзафлопсный компьютер, программа для которого будет использовать менее 1% от его производительности — для суперкомпьютера рекордной производительности и программное обеспечение должно быть рекордным.

Пока матрицы имеют небольшую размерность, они целиком попадают в кэш-память. Время передачи данных между оперативной и кэш-памятью равно 3*n² * T_M, где 3*n² — количество чисел в трех матрицах, T_M — время чтения числа из оперативной памяти в кэш. Время (последовательного) выполнения арифметических операций равно n³ T_a , где T_a — время выполнения умножения и сложения. Соотношение между временем обращения к памяти и временем выполнения арифметических операций у разных компьютеров разное, поэтому примем правдоподобную величину T_M = 15*T_a. Если n — велико, но не настолько, чтобы все матрицы не помещались полностью в кэш, то время передачи данных между оперативной памятью и кэш-памятью несущественно — в этом случае распараллеливание может помочь ускорить выполнение программы.

Объем кэш-памяти процессоров настольных компьютеров — 8 Мбайт, и три матрицы восьмибайтных чисел размера n = 600 точно не поместятся в кэш, и тогда процессор постарается разместить в кэш данные, необходимые для выполнения двух самых вложенных циклов — циклов со счетчиками j и k. Это целиком матрица B, пусть целиком матрица X и строка матрицы A. В этом случае каждая строка матрицы A читается n раз. Следовательно, время обращения к памяти теперь (2*n² + n³ )* T_M. Время выполнения арифметических операций остается прежним n³ T_a и становится значительно меньше времени обращений к памяти. Поэтому уменьшение времени выполнения арифметических операций за счет распараллеливания относительно общего времени выполнения уже несущественно, а если предположить, что только одна из трех матриц помещается в кэш-памяти, то ускорения от распараллеливания вообще нет. Более того, затраты на синхронизацию параллельных процессов могут дать замедление. Производительность выполнения данной программы значительно ниже пиковой производительности процессора. Не окажется ли реальная производительность экзафлопсной системы на несколько порядков меньше пиковой?

Когда процессоры были большими

В течение нескольких десятилетий с каждым годом уменьшается отношение времени выполнения арифметической операции ко времени считывания ее аргументов из памяти. Пятьдесят лет назад, когда проектировались первые параллельные системы, умножение вещественных чисел было дорогостоящей операцией, а время обращения к памяти было незначительно. Сегодня — все наоборот: обращение к памяти осуществляется на порядки дольше, чем выполнение арифметических операций. Это привело к необходимости использования модулей вспомогательной памяти с быстрым доступом: кэш-память, управляемая процессором; локальная программируемая память (например, у процессора IBM Cell) и т. д. Так возникла иерархия памяти (рис.1), которую теперь необходимо учитывать при разработке быстрых программ. Но чтобы малая по объему быстрая память эффективно использовалась, необходимо, чтобы попадающие в нее данные многократно использовались алгоритмом, — только в этом случае параллельное программирование может принести хоть какое-то ускорение. А если используются модули быстрой памяти нескольких уровней, то данное правило надо учитывать на каждом уровне.

Рис. 1. Вычислительная система с тремя уровнями памяти. Задача дважды разбивается на подзадачи

Традиционно сложностью алгоритма считается количество вычислительных операций как функция от количества входных данных. Анализ сложности алгоритмов получил активное развитие в середине 60-х годов прошлого века после публикации известного алгоритма Штрассена перемножения матриц. Этот алгоритм оказался лучше стандартного по количеству операций при перемножении матриц и давал преимущества на матрицах размера меньше 100. Но сегодня сложность должна оцениваться обращениями к памяти.

Иерархия и процесс вычисления

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

• вычислительное ядро со своей локальной кэш-памятью;
• многоядерный процессор: множество ядер на кристалле с общей памятью на этом же кристалле;
• вычислительный узел: несколько кристаллов с одной общей оперативной памятью;
• кластер: несколько вычислительных узлов, соединенных коммуникационной сетью (возможно, с общей внешней памятью);
• сеть, объединяющая множество компьютеров.

Чем ближе по сети коммуникаций модуль памяти к вычислительному устройству, тем он меньше по объему, но быстрее по доступу к данным. Пересылки данных между модулями памяти тем медленнее, чем больше модули памяти.

Рис. 2. Вычислительная система со сложно организованной памятью и коммуникациями

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

Иерархия архитектуры предполагает адекватное ей разбиение задачи на подзадачи — на вычислительном устройстве каждого уровня архитектуры решается подзадача соответствующего уровня иерархии подзадач. Разбиение задач на подзадачи необходимо даже для простой (рис. 1) иерархической структуры памяти, состоящей из большой, но медленной и из быстрой, но малой по объему памяти. Адекватное разбиение задачи на подзадачи определяет эффективность всей высокопроизводительной системы. Компиляторы и процессоры выполняют такое разбиение автоматически, но практически без минимизации обращений к памяти — хорошую оптимизацию программист должен выполнять сам.

Таким образом, основная цель разбиения задач на подзадачи состоит в том, чтобы данные каждой подзадачи помещались в памяти соответствующего по иерархии вычислительного устройства, а самих подзадач было как можно меньше. Однако оба условия противоречат друг другу — например, в качестве подзадачи можно взять выполнение одной операции, и тогда будет выполняться второе условие, но подзадач будет столько же, сколько и операций. Если в качестве подзадачи взять саму исходную задачу, то ее данные могут не поместиться в кэш-памяти. Для достижения баланса между этими условиями требуется соблюдать следующие принципы разбиения задачи на подзадачи:

• подзадачи должны формироваться так, чтобы с быстрой памятью ограниченного объема можно было выполнять много операций;
• каждая подзадача должна быть максимально оптимизирована для своего уровня иерархии.

В работе [2] рекомендуется разбивать задачу на подзадачи, уже начиная от модулей памяти, наиболее близких к вычислительным устройствам, и передвигаясь к большим модулям памяти. Этот подход привел автора пакета GotoBLAS, японского программиста Казушиге Гото, к рекордной в свое время производительности программы перемножения матриц. Однако сегодня этот результат превзойден в пакете PLASMA Джека Донгарры, и достичь этого удалось за счет нестандартных размещений данных при соблюдении принципа разбиения задач на подзадачи от быстрой памяти к большой.

Стандартные распределения против производительности

Сегодня почти нет технологических препятствий для создания на одном кристалле тысяч вычислительных ядер, однако медленная скорость чтения данных из памяти создает «стену памяти» — вычислительные ядра не успевают получать данные на обработку, и с каждым годом эта проблема все обостряется. Для оптимизации применения иерархии памяти приходится использовать не только метод разбиения задач на подзадачи — на эффективность обращений к памяти влияет и способ размещения данных в этой памяти.

Из большой памяти в быструю данные обычно читаются порциями (кэш-линейками, страницами и т. п.), и при считывании нужной порции данных заодно считываются и соседние. Важно разместить данные в памяти так, чтобы и соседи оказались востребованы для обработки в следующий момент, а иначе налицо кэш-промах. Минимизация таких промахов — одна из задач разработчика эффективных программ. Хорошо известна техника выравнивания (alignment) — распределение данных в памяти с позиций, адреса которых кратны длине кэш-линейки, но и она оставляет еще много кэш-промахов. Существуют стандарты распределения данных в оперативной памяти для компиляторов. Например, компилятор языка Фортран размещает матрицы в оперативной памяти по столбцам, а Си — по строкам. Но при блочной обработке матриц желательно, чтобы элементы каждого блока (а не строки или столбца) оказывались в памяти рядом. В этом случае стандартное распределение массива может создавать кэш-промахи при обращении к нему.

На данный момент существует множество алгоритмов перемножения матриц, реализованных в стандартных пакетах LAPACK, ATLAS, MKL, OpenBLAS и т. п., но главным ограничением их всех является то, что они принимают на вход только матрицы, размещенные стандартным способом (по строкам либо по столбцам). В то же время при блочном размещении матриц можно существенно улучшить использование иерархии памяти — уменьшить количество промахов, увеличить эффективность неявной предвыборки данных из оперативной памяти. Пакет PLASMA (использующийся, в частности, для тестирования компьютерных систем на попадание в рейтинг Top500 и, естественно, оптимизированный для вычислительных систем с общей памятью) применяет именно блочное размещение матриц вместо стандартного. При таком размещении элементы каждого блока матрицы располагаются в линейной оперативной памяти рядом, что способствует минимизации кэш-промахов.

Основа рекордного быстродействия

Многие современные процессоры содержат несколько уровней кэш-памяти — из быстрой памяти данные по пути к регистрам попадают в очень быструю, и для организации эффективных вычислений на такой архитектуре иногда выгодно произвести разбиение подзадач на подзадачи второго уровня, что может привести к двойному блочному размещению данных в памяти. Алгоритм перемножения матриц, осуществленный в распараллеливающей системе OPS с дважды блочно размещенными матрицами, был реализован авторами и оптимизирован для процессоров семейства SandyBridge. Оказалось, что данный подход дает прирост производительности примерно 5–7% по сравнению с аналогом из пакета PLASMA (рис. 3), хотя и этот пакет уже обгоняет по производительности популярный пакет линейной алгебры Intel MKL. Учитывая то, что к перемножению матриц сводятся многие алгоритмы линейной алгебры, такое увеличение производительности может сильно повысить скорость выполнения программы. Для экзафлопсного суперкомпьютера 5% производительности означает снижение затрат на ресурсы на 5% (процессоров, коммуникаций и т. п.), а если компьютер стоит миллиард, то и экономия составит 50 млн, что уж точно окупит работу программиста.

Рис. 3. График сравнения производительности алгоритма OPS и пакета PLASMA

В предложенном алгоритме учитывается иерархическая структура кэша данных (первый, второй и третий уровни), поэтому исходный стандартный алгоритм перемножения матриц был заменен на эквивалентный блочный. Размеры блоков подбираются таким образом, чтобы все элементы матриц, используемые при перемножении блоков, полностью поместились в кэш-память соответствующего уровня. Далее все три матрицы нужно разместить блочно в оперативной памяти, что позволяет уменьшить количество кэш-промахов и увеличивает эффективность неявной предвыборки данных из оперативной памяти.

Внутренний цикл, в котором производится перемножение элементов и накапливается результат, можно векторизовать. Самый эффективный алгоритм векторизации внутреннего цикла в задаче перемножения матриц был опубликован в [2] и реализован в библиотеке GotoBLAS. Алгоритм использует методику register packing, при которой данные внутри блоков каждой из матриц размещаются полосами. Ширина полос подбирается исходя из размеров векторных регистров целевой архитектуры. При перемножении матриц C=A*B для AVX (Advanced Vector Extensions — набор векторных команд, оперирующих с 256-битными регистрами) ширина полосы для матрицы-результата (C) и матрицы, перемножаемой слева (A), равна 4, а для матрицы, перемножаемой справа (B), — 8.

Для того чтобы ускорить результирующую программу, алгоритм перемножения полос был запрограммирован на языке Ассемблера, использовались стандартные методы низкоуровневой оптимизации, такие как программная предвыборка данных, оптимизация индексных выражений массивов, увеличение дистанции между операциями, имеющими зависимости по данным. Рекордная производительность без использования низкоуровневых программ недостижима.

Размещения с перекрытиями в распределенной памяти

Экзафлопсный суперкомпьютер обязательно будет иметь распределенную память. Есть ли здесь резервы повышения производительности?

В системах с распределенной памятью самая медленная операция — обращение к локальной памяти другого процессора, и для минимизации межпроцессорных пересылок в некоторых случаях можно использовать размещения данных с перекрытиями, суть которых в том, чтобы, потратив память на дублирование некоторых данных (до 7%), выиграть в быстродействии (несколько десятков процентов). В качестве примера рассмотрим задачу итерационного умножения ленточной матрицы на вектор.

В стандартном параллельном алгоритме такая матрица разрезается на равные полосы, вектор разделяется на равные части и на каждой итерации происходят пересылки элементов вектора между узлами. При размещении данных с перекрытием алгоритм видоизменяется: пересылки координат вектора происходят не на каждой итерации, а матрица разрезается так, что некоторые ее элементы могут принадлежать нескольким узлам одновременно.

Если перед нечетной итерацией процесса данные распределены с перекрытием, то после выполнения умножения координаты вектора окажутся распределены так, что четную итерацию можно выполнить без пересылки. Численные эксперименты размещения данных с перекрытием проводились на кластере INFINI ЮГИНФО ЮФУ для итерационного умножения матрицы на вектор и для метода Якоби решения задачи Дирихле для уравнения Лапласа.

Для итерационного умножения ленточной матрицы на вектор (итераций: 10, размер вектора: 65536, ширина ленты: 513) размещение данных с перекрытием показывает выигрыш быстродействия на 22–24% по сравнению со стандартным. Для метода Якоби решения трехмерной задачи Дирихле для уравнения Лапласа (размерность сетки: 500*500*500, количество итераций: 2000) размещение данных с перекрытием показывает выигрыш быстродействия в два раза.

Цена вопроса

Рекордную производительность на конкретной архитектуре достигнуть непросто — например, рекордная программа перемножения матриц занимает порядка 1000 строк кода, большинство которых завязаны на особенности конкретной архитектуры процессора. Это примерно в 100 раз больше объема кода наивной программы, представленной в начале статьи.

Программы с оптимизированным размещением данных (с перекрытиями) в распределенной памяти также требуют существенно большего объема кода, чем параллельные программы, реализующие стандартные размещения (например, размещение полосами). Для решения задачи Дирихле параллельная программа, использующая размещения данных с перекрытиями, имеет объем кода 33 Кбайт, что в два раза больше обычной параллельной программы.

Кроме того, для для разработки таких программ требуются программисты очень высокой квалификации. Сокращение времени разработки может быть достигнуто путем автоматической оптимизации, однако разработка таких средств сегодня находится в зачаточном состоянии.

***

Известны случаи закупок суперкомпьютеров без специального программного обеспечения или даже инженерных инфраструктур (систем охлаждения), но сколько еще можно наступать на одни и те же грабли? Разработка программного обеспечения для зкзафлопсного суперкомпьютера — сложная задача, требующая многолетних усилий программистов высокой квалификации.

Создание экзафлопсного суперкомпьютера бессмысленно без анализа специфики его будущих приложений и особенностей их отображения на архитектуру. При этом оценивать следует не только количество вычислительных операций, но и количество обращений к памяти на всех уровнях иерархий — объемы и скорости доступа всех модулей памяти должны быть достаточны, чтобы обеспечить данными вычислительные устройства на всех этапах алгоритма. Уместно моделирование времени выполнения различных версий будущих алгоритмов на прототипах экзафлопсных систем, причем такое моделирование следует начинать уже сейчас. Некоторые принципы разработки эффективного программного обеспечения для перспективных вычислительных систем можно сформулировать уже сегодня — описание иерархии подзадач должно соответствовать иерархии архитектуры суперкомпьютера, а структуры хранения данных и инструменты автоматизации разработки быстрых программ должны быть ориентированы на минимизацию обращений программы к памяти.

Литература

1. Штейнберг Б.Я. Зависимость оптимального распределения площади кристалла процессора между памятью и вычислительными ядрами от алгоритма // Труды международной конференции PACO’2012. Параллельные вычисления и задачи управления. 24-26 октября. М.: ИПУ РАН, 2012. с. 99-108.
2. Kazushige Goto, Robert A. van de Geijn. Anatomy of High-Performance Matrix Multiplication. ACM Trans. Math. Softw., Vol. 34, No. 3. (May 2008), pp. 1-25.

Лев Гервич (lgervith@gmail.com), Борис Штейнберг (borsteinb@mail.ru), Михаил Юрушкин (m.yurushkin@gmail.com) — Южный федеральный университет (Ростов-на-Дону). Статья подготовлена на основе материалов доклада, представленного авторами на IV Московский суперкомпьютерный форум (МСКФ-2013, грант РФФИ 13-07-06046г).