Вычислительные ресурсы процессоров NeuroMatrix с плавающей точкой в задачах обработки больших потоков данных
УДК 004.383.3 / 519.684.6
DOI: 10.22184/1993-8578.2018.82.110.118
Современные архитектуры процессоров обладают сложной распределенной системой памяти, могут иметь гетерогенную структуру с несколькими процессорами. Каждый из процессоров может включать несколько сложно-функциональных вычислительных узлов. Для эффективного программирования специализированных процессоров необходимо учитывать структуру и принципы взаимодействия отдельных узлов.
В данной статье будут рассматриваться вопросы программирования задач по обработке больших потоков данных с плавающей точкой на процессоре NM6407 производства НТЦ «Модуль». Будут описаны основные вычислительные узлы, режимы и особенности их использования в наиболее типовых задачах, а также принципы распараллеливания, как по вычислительной нагрузке, так и по взаимодействию с памятью.
АРХИТЕКТУРА ПРОЦЕССОРА
Рассмотрим процессор, начиная от обобщенной структуры и заканчивая уровнем элементарных базовых операций. На рис. 1 представлена структурная модель процессора: два независимых процессорных узла с плавающей и фиксированной точкой NMPU0 и NMPU1. Каждый узел имеет собственную накристальную однотактовую память, DMA-контроллер и DDR-интерфейс c внешней памятью. Доступ в память осуществляется по 64-разрядным шинам. Обмен данными между процессорными узлами может осуществляться через общую разделяемую память (Shared Memory 0, 1).
Процессорный узел NMPU0 с плавающей точкой состоит из управляющего RISC-ядра и 4-х вычислительных ячеек (FPU Cell0-Cell3). Целочисленный процессорный узел NMPU1 также включает свое RISC ядро и матрично-векторный умножитель [1].
На рис. 2 показана схема процессорного узла NMPU0. Внутренняя память разделена на 8 банков по 128 кБ псевдо-двухпортовой памяти, к которым может быть осуществлен одновременный внешний доступ (например, со стороны DMA или коммуникационных портов) и обращение со стороны процессорного ядра. Память и процессорные ячейки объединены между собой коммутатором (switch 9 → 11), обеспечивающим параллельное обращение сразу к нескольким банкам памяти. Как видно из рис. 2, каждая процессорная ячейка имеет две входных шины и одну выходную, что позволяет ячейке параллельно принимать два потока данных и отдавать один. Процессорные ячейки осуществляют доступ в память через шинный коммутатор 9 → 11. В общей сложности на все 4 ячейки (Cell0–Cell3) коммутатор обеспечивает максимум 4 одновременных доступа на чтение из памяти и 2 на запись.
Каждая ячейка FPU Cell0–FPU Cell3 состоит из операционного узла и восьми векторных регистров, вмещающих до 32-х 64-разрядных регистровых ячеек (см. рис. 3). В зависимости от режима работы операционного узла каждая регистровая ячейка (reg0, reg1–reg31) векторного регистра интерпретируется либо как 64-разрядное число с плавающей точкой двойной точности, либо как пара чисел одинарной точности, либо как комплексное число одинарной точности. На рис. 3 показано, как векторные регистры выступают в качестве входных и выходных операндов для векторных инструкций операционного узла (Operation Unit). Коммутатор 8 → 4 позволяет одновременно выбрать любые 3 регистра в качестве входных операндов, а из 4-го производить выгрузку в память. Коммутатор 3 → 8 обеспечивает на фоне вычислений параллельную загрузку двух регистров и одного регистра с выхода вычислителя по обратной связи новой порцией данных. Пересылки могут быть как между памятью и регистрами, так и между регистрами процессорных ячеек. Данные из ячеек (reg0, reg1. reg31) векторных регистров выбираются последовательно с темпом в один такт и обрабатываются заданной инструкцией.
Как показано на рис. 4, операционный узел работает с векторами трех типов: данные двойной точности, комплексные числа одинарной точности и числа одинарной точности. Операционный узел выполняет над векторными регистрами SIMD операции сложения, вычитания, умножения, взятие модуля, выделение знака, маскирование, а также операцию вида A · х + b. В режиме комплексных чисел (мнимая и действительная часть по 32 бита) или в режиме 64-разрядных чисел двойной точности A, x, b — это обычные числа. В случае одинарной точности A является матрицей 2 × 2, x и b-вектора из двух 32-разрядных элементов с плавающей точкой. Выход с операционного вычислителя замкнут на вход через коммутатор 3 → 8, что позволяет накапливать результат вычислений в векторных регистрах.
В каждом процессорном такте происходит последовательная выборка данных из ячеек (reg0, reg1. reg31) векторных регистров (vreg0.vreg8), между которыми производится заданное командой действие. Операции между векторными регистрами могут выполняться как поэлементно (например, при суммировании двух векторных регистров), так и выборочно, где в качестве второго операнда используется не весь векторный регистр, а только одна регистровая ячейка (например, при прибавлении к вектору числа).
В режиме двойной точности или в комплексном режиме в каждом такте производится операция над отдельными числами. В режиме одинарной точности операционный узел также может выполнять умножение с накоплением матрицы 2 × 2 на вектор из 2 элементов. Схематично процесс умножения двух чисел с накоплением показан на рис. 4а. На рис. 4б показана схема умножителя в случае режима одинарной точности, где в весовые коэффициенты матрицы 2 × 2 (B0, B1, B2, B3) загружаются числа из пары векторных регистров [vreg2, verg3].
Существует четыре способа загрузки коэффициентов в матричный умножитель (см. рис. 5):
вертикальная загрузка (коэффициенты [B0, B1] лежат смежно в vreg2, [B2, B3] лежат смежно в vreg3);
горизонтальная загрузка (коэффициенты [B0, B2] лежат смежно в vreg2, [B1, B3] лежат смежно в vreg3);
диагональная загрузка (коэффициенты [B1, B2] лежат смежно в vreg2);
загрузка для комплексного умножения (коэффициенты [B.im,B.re] лежат смежно в vreg2).
На рис. 6. показаны варианты загрузки в весовые коэффициенты матричного умножителя и схема умножения. В зависимости от режима вычислитель осуществляет различные математические операции. Вертикальная загрузка позволяет выполнять транспонирование или умножать матрицу на вектор-столбец. Горизонтальная загрузка — умножение вектор-строки на матрицу. При диагональной загрузке осуществляется поэлементное умножение. Специальный инвертирующий вход сумматора на рис. 6в позволяет осуществлять комплексное умножение.
УМНОЖЕНИЕ МАТРИЦЫ НА ВЕКТОР-СТОЛБЕЦ
Рассмотрев элементарную команду умножения, распространим процесс до полного умножения матрицы на вектор произвольного размера. Несмотря на то что процессор поддерживает SIMD инструкции с числом повторений до 32, для упрощения опишем процесс умножения, используя 4 векторных регистра глубиной только в одно 64-разрядное слово:
регистр, накапливающий результат (пара чисел [D00, D01] на рис. 7);
регистр, содержащий данные вектор-столбца А ([А20, A30] на рис. 7);
региcтровая пара для хранения матрицы 2 × 2 коэффициентов B ([B02, B03] и [B12, B13]).
На рис. 7 показаны матрица B размером 8 × 8 и вектор A из 8 элементов в памяти и порядок загрузки данных в векторные регистры (регистры обозначены цветными овальными прямоугольниками) с последующим умножением и накоплением в выходной регистр. Процесс вычисления пары чисел D00, D01 состоит из 4 шагов. На рис. 7 справа развернута схема умножения второго шага. Если повторить данную операцию аналогично еще 3 раза, получится полное умножение матрицы на вектор. Отметим, что загрузка множителей B происходит вертикально (обозначено пунктиром). Другими словами, матрица коэффициентов 2 × 2 загружается в умножитель транспонированной, что математически позволяет реализовать умножение матрицы на вектор-столбец.
УМНОЖЕНИЕ ВЕКТОР-СТРОКИ НА МАТРИЦУ
Для осуществления умножения вектор-строки на матрицу загрузка коэффициентов матрицы B производится горизонтально по две пары чисел: [B20, B21] и [B31, B30] (объединены пунктиром на рис. 8). Процесс вычисления пары чисел [D00, D01] осуществляется также за четыре шага и показан на рис. 8.
ОПТИМИЗАЦИЯ РАБОТЫ С ПАМЯТЬЮ
В приведенных примерах присутствуют три потока данных: два входных для чтения A и B и один выходной для записи результата D. Поскольку на каждом шаге умножения требуется две пары чисел матрицы B, то на подкачку данных уходит 2 такта. Если данные массивов A, B и D размещены по разным банкам памяти, то обращение к ним будет параллельно и умножение займет 2 такта. Однако, если матрицу B разбить на четные и нечетные строки и хранить их в разных банках, то загрузку можно ускорить до одного такта.
В данных примерах было рассмотрено умножение векторных регистров глубиной в одно слово. SIMD инструкции с малым числом повторения (<4) менее эффективны из-за задержек в конвейере, поэтому значения в тактах приведены условно, для наглядности. Однако указанные цифры соответствуют действительности для SIMD инструкций с большим числом повторений. При таком допущении справедливо оценивать суммарную производительность, принимая время исполнения векторной команды за один такт. Далее рассмотрим использование SIMD инструкций с числом повторений 8 на примере перемножения матриц 8 × 8.
ПЕРЕМНОЖЕНИЕ МАТРИЦ НА ОДНОЙ ЯЧЕЙКЕ
Умножение двух матриц строится на принципах умножения вектор-строки на матрицу, описанного выше.
Главным отличием является то, что теперь на первом шаге две пары чисел ([B00, B01] и [B10, B11]) матрицы B можно умножать сразу на два столбца ([Аx0, Ax1]) матрицы А, которые целиком помещаются в векторный регистр глубиной в 8 слов. На 2-м шаге матрица ([B20, B21] и [B31, B30]) умножается на следующие два столбца ([Аx2, Ax3]) с прибавлением к вектору-результату от первого шага и так далее. За 4 шага полностью осуществляется процесс умножения двух матриц. В отличие от векторно-матричного умножения, в данном случае не требуется непрерывная перезагрузка коэффициентов матричного умножителя. Так как она хранится в регистрах в течение 8 тактов умножения, ее подкачка осуществляется на фоне непрерывного чтения матрицы А и не оказывает влияния на производительность.
На рис. 10 показана временная диаграмма процессов умножения и обращения в память. Из рисунка видно, что если массивы расположены в разных банках, то доступ к ним осуществляется параллельно на фоне умножения. За счет так называемого «зацепления» производительность всей функции определяется лишь временем чтения матрицы А и временем разгона конвейера.
РАСПАРАЛЛЕЛИВАНИЕ МАТРИЧНОГО УМНОЖЕНИЯ НА 4 ЯЧЕЙКИ
Сопроцессор с плавающей точкой содержит 4 вычислительных ячейки, поэтому предыдущую схему умножения распространим на все 4 ячейки следующим образом: матрица B вертикально разбивается на 4 равных подматрицы (см. рис. 11). Каждое умножение подматрицы производится согласно выше описанной схеме. При этом, так как во всех четырех ячейках участвуют одни и те же данные матрицы А, обращение к памяти в массив A можно сократить в 4 раза за счет регистровых пересылок между ячейками (см. рис. 11). При матрично-матричном умножении не требуется непрерывная подкачка данных матрицы B. Во время умножения в каждой ячейке образуются временные окна, которые могут использоваться для чтения данных B другими ячейками. На рис. 12 показана временная диаграмма обращений в память каждой ячейкой. Видно, что благодаря «зацеплению», конфликтов по обращению в память не происходит. В результате, каждая ячейка совершает 4 умножения с накоплением (MAC) с темпом в один такт, а производительность всего процессора достигает 16 MAC или 32 FLOPs за такт. Следует отметить, что значение 16 MAC является пиковой производительностью. Реальные цифры меньше и зависят от размеров матриц, где сказываются потери, обусловленные расходами на вызов функции, инициализацию и разгон конвейера вычислений. При больших размерах матриц реальная производительность приближается к предельной и достигает 15,8 MAC за такт (31,6 FLOPS).
АНАЛИЗ ПРОИЗВОДИТЕЛЬНОСТИ ОТ РАЗМЕРА ДАННЫХ
Для процессора NM6407 были реализованы простейшие функции линейной алгебры из состава библиотеки BLAS. Производительность функций определяется тремя факторами:
1 — числом задействованных вычислительных ячеек в задаче;
2 — распределением данных по разным банкам;
3 — размером данных.
Рассмотрим простейшие функции линейной алгебры из состава библиотеки BLAS: умножение вектора на скаляр y͢ = ax͢ и поэлементное перемножение векторов z͢ = x͢y͢. Данные функции являются достаточно простыми — всего с одной операцией умножения на каждый элемент. Для умножения вектора на константу и поэлементного умножения чисел одинарной точности ячейка работает в режиме диагональной загрузки (рис. 6б) и способна осуществлять по 2 умножения за такт. В данном случае максимальная производительность ограничивается лишь скоростью обращения за данными и составляет 2 FLOPs за такт, остальные 3 ячейки из 4 остаются незадействованными. Для функций на рис. 13 приведен график зависимости усредненного числа операций за такт от длины векторов. Как видно из графика, при большом объеме данных (более 1000 элементов) реальная производительность приближается к 2 FLOPs за такт. При малых размерах (менее 1000) сказывается влияние потерь на вызов С-функции, чтение параметров и прочих подготовительных операций.
В случае скалярного произведения dot = (x͢ · y͢) ячейка загружается вертикально (пара чисел вектора x͢ загружается в правую половину и нули в левую, см. рис. 6а). Ячейка производит 2 умножения с накоплением, что соответствует 4 FLOPs за такт.
Рассмотрим более сложную функцию вида y͢ = aAx͢ + by͢, где A матрица размера h на L a, b скаляры, а y͢ и x͢ векторы длиной h. Данную функцию можно реализовать, задействовав уже три ячейки. Первая вычисляет Ax͢, результат с помощью регистровой пересылки поступает на вторую, где умножается на a, а третья ячейка вычисляет by͢ и прибавляет к результату от второй. Первая ячейка выполняет матрично-векторное умножение, которое является самой длинной операцией. 2-я и 3-я ячейки загружаются лишь по мере готовности вектора результата от первой. Производительность, главным образом, определяется временем вычисления Ax͢. Для больших размеров А производительность можно приблизительно принять 4 FLOPs за такт.
Максимально задействовать все 4 ячейки удается при матричном умножении, где достигается максимально возможная пиковая производительность 32 FLOPs за такт.
В табл. 1 приведены значения предельной (расчетной) и реальной производительностей некоторых функций линейной алгебры из состава библиотеки BLAS для данных одинарной точности. Значения производительности были замерены на максимальных объемах данных, помещающихся во внутренних банках памяти процессора, что соответствует размеру векторов 32000 или матриц 180x180.
Как видно из полученных результатов, реальная производительность достигает 95–98 % от теоретически возможной. Приведенные функции жестко поддерживают интерфейс стандартной библиотеки BLAS (Basic Linear Algebra Subprograms), что в некоторых случаях при их реализации ограничивает программиста в распределении вычислений на несколько ячеек, а также не позволяет максимально распараллелить потоки входных и выходных данных. Во второй колонке таблицы указано количество задействованных ячеек для разных функций, где видно, что оно не везде равно максимальному числу. Максимальная производительность достигает предельного пикового значения — (32 FLOPs/cycle) для функций матричного умножения, где участвуют все 4 ячейки.
Однако при реализации частной задачи, как правило, всегда есть возможность более эффективно распределить данные между банками памяти, распараллелив обмен с памятью. Также можно сгруппировать вычисления до более сложных выражений и задействовать таким образом большее число вычислительных ячеек. Так, в целях поддержки процессора математическими функциями, помимо библиотеки BLAS, разрабатываются аналогичные, либо более сложные функции с максимальной адаптацией под архитектуру процессора. Данные функции разрабатываются в составе библиотеки NMPP [2]. Интерфейс этих функций уже отличается от библиотеки BLAS, но в некоторых случаях позволяет значительно повысить производительность (в 2 и более раз).
РЕАЛИЗАЦИЯ БЫСТРОГО ПРЕОБРАЗОВАНИЯ ФУРЬЕ (БПФ)
В заключение рассмотрим реализацию алгоритмически сложных потоковых вычислений на примере БПФ-256. Стандартным подходом для реализации БПФ является использование схемы Кули — Тьюки [3]. Однако для процессора NM6407 она не совсем эффективна. Из-за сложного порядка выборки данных на первых слоях алгоритма не удается использовать SIMD инструкции с большим числом повторений. Кроме того, алгоритм необходимо адаптировать под 4 вычислительные ячейки и 8-банковую структуру внутренней памяти. Наиболее оптимальным вариантом получилась реализация комбинированной схемы, где первый слой реализован через последовательность дискретных преобразований Фурье по основанию 8 см (рис. 14). Последующие слои представляют набор классических «бабочек» с комплексным умножением, сложением и вычитанием. Подробно алгоритм БПФ описан в статье [4]. В данном разделе наибольший интерес представляет блок ДПФ-8, поэтому разберем только его реализацию.
Блок ДПФ представляет собой умножение матрицы 8 × 8 на вектор, описанный ранее, однако, поскольку таких блоков в алгоритме БПФ-256 — 32, появляется возможность более эффективно задействовать процессорные ячейки.
Как уже сказано, ДПФ-8 можно вычислить на NM6407 так же, как и умножение матрицы на вектор. На двух процессорных ячейках вычисляется один блок ДПФ-8. Соответственно 4 ячейки вычисляют сразу 2 блока ДПФ-8, а матрица коэффициентов W для каждого ДПФ-8 одна и та же и хранится в 2-х ячейках FPU Cell одновременно (рис. 15).
Рассмотрим более подробно вычисление блока ДПФ-8. Столбцы матрицы W — W(, 1), W(, 2), W(, 3) и W(, 4) — загружаются из памяти в векторные регистры: VR1, VR2, VR3 и VR4 ячейки CELL0 (рис. 16). Аналогично из памяти в регистр VR0 той же ячейки попадают первые 4 элемента вектора x͢. После чего происходит поэлементное комплексное умножение с накоплением каждого столбца матрицы, находящегося в векторном регистре процессорной ячейки, на один из первых 4-х элементов вектора x͢, лежащих в VR0. Накопление результата происходит в регистре VR5. Поочередный выбор действующего элемента регистра VR0 на каждом этапе вычислений показан розовым квадратом на рис. 16 и производится автоматически с помощью специальной команды.
Такие же операции, как в CELL0, выполняются и в CELL1. Причем 8-элементный вектор, полученный в ходе вычислений на CELL0, копируется в один из регистров CELL1 и затем суммируется с регистром VR5 ячейки CELL1. После этого в регистре VR5 находится 8-элементный вектор: ДПФ-8.
При вычислении серии блоков ДПФ-8 (от 100 блоков и больше) время, затраченное на вычисление одного такого блока, равно 83 тактам. Причем коэффициенты матрицы W хранятся в векторных регистрах ячеек и повторно используются при вычислении каждого блока ДПФ-8 из всей серии. При однократном вычислении такого блока число тактов увеличивается до 127. Так получается, что векторный сопроцессор с плавающей точкой эффективнее работает с большими потоками данных, в данном случае с серией блоков ДПФ-8.
В табл. 2 приведены значения производительности функций БПФ в тактах для процессоров NM6407, Texas Instruments C674x, Intel Pentium 4 и ARM Cortex A-15. На основе этих данных можно сделать вывод, что NM6407 обгоняет все процессоры из таблицы при вычислении БПФ-256. Остальные функции БПФ NM6407 вычисляет быстрее Pentium 4 в среднем на 22 %.
ЗАКЛЮЧЕНИЕ
При работе с большими массивами данных за счет распределенной структуры памяти на 8 банков и конвейерного выполнения команд процессор показывает высокую степень распараллеливания потоков данных, в которых реальная производительность достигает порядка 95–98 % от расчетной. При максимальной загрузке процессорных ячеек в операциях типа матричного умножения, производительность достигает 98 % от пиковой. Анализ и сравнение производительностей, представленных в табл. 1 и 2, показывают высокую эффективность процессора в задачах обработки больших потоков данных, однако для достижения максимально высоких показателей необходимо учитывать все особенности процессора и адаптировать вычисления под его архитектуру.
ЛИТЕРАТУРА
1. Косоруков Д. Е., Эйсымонт А. Л., Осипов В. Г., Панфилов А. П., Черников В. М., Виксне П. Е., Шелухин А. М., Насонов И. И. «СБИС на базе ядра NMC3 для программного приемника навигационных сигналов», Сборник докладов Международной конференции «Микроэлектроника 2015».
2. URL: https://github.com/RC-MODULE/nmpp.
3. Отнес Р., Эноксен Л. «Прикладной анализ временных рядов». — М.: Мир, 1982.
4. «Адаптация алгоритма БПФ для матрично-векторного сопроцессора nm6407 с плавающей точкой», 19-я Международная конференции «Цифровая обработка сигналов и ее применение».