Как мы протестировали оптимизацию библиотеки OpenBLAS под RISC-V, исправив более миллиона тестов

В рамках инициативы EuroHPC европейские HPC-центры, такие как Barcelona Supercomputing Center (BSC) и Edinburgh Parallel Computing Centre (EPCC), открыли центры компетенции RISC-V. При этом решения на x86- и ARM-архитектурах не признаются частью европейской инициативы по развитию собственной технологической независимости.

На горизонте появляются первые RISC-V HPC CPU. Уже представлены процессоры SG2042 с 64 ядрами и Veyron V1 с чиплетами до 16 ядер. Первые RISC-V HPC-системы ожидаются в 2025—2026 годах в рамках проекта европейского суперкомпьютера BSC MareNostrum 6. Также реализован первый GPU на RISC-V, работающий на FPGA и поддерживающий OpenCL.

Как я уже писал выше, для успешного развития нужно больше компаний, которые бы помогали в развитии экосистемы решений для HPC. YADRO — одна из них. Мы подробно изучили OpenBLAS — известную open source-библиотеку, первую в мире из HPC-сегмента, портированную и оптимизированную на RISC-V.

У библиотеки много преимуществ:

• Она полностью реализует интерфейсы BLAS и LAPACK, а также предоставляет дополнительное множество BLAS-like функций для расширения функциональности BLAS.
• BLAS-часть библиотеки содержит в себе набор базовых операций линейной алгебры, которые чаще всего используются в прикладных программах. Все функции оптимизированы как алгоритмически, так и низкоуровнево, а на их основе строятся реализации более сложных алгоритмов, входящих в LAPACK.
• Поддерживается большим сообществом, в составе которого не только независимые разработчики, но и крупные компании.
• Распространяется под 3-Clause BSD License, которая допускает свободное использование в коммерческих целях.

Но главным преимуществом, конечно, стало наличие оптимизаций для RISC-V. OpenBLAS содержит оптимизации с использованием инструкций из RISC-V Vector Extension как для версии 0.7.1, так и для 1.0. Это позволяет не останавливаться на оптимизации, а сразу перейти к модификации уже имеющихся ядер и детальной настройке алгоритмов для требуемого железа.

Перед использованием лучше убедиться в функциональной корректности существующих оптимизаций, для чего в OpenBLAS есть множество тестов. До начала тестирования библиотеки немного расскажу про BLAS и LAPACK.

BLAS (c англ. Basic Linear Algebra Subprograms — базовые подпрограммы линейной алгебры) является де-факто стандартом интерфейса базовых операций линейной алгебры. BLAS API используется в многих высокоуровневых библиотек и языках программирования, а без линейной алгебры было бы невозможно развитие data science, машинного обучения и искусственного интеллекта, современных поисковых систем, робототехники и других перспективных направлений.

Функции BLAS делятся на три группы, называемые «уровнями»:

1. Операции над векторами — например, скалярное произведение или умножение вектора на скаляр.
2. Векторно-матричные операции — например, умножение матрицы на вектор, причем вид матрицы учитывается при оптимизации алгоритма.
3. Матрично-матричные операции — например, умножение матриц или решение системы линейных алгебраических уравнений (СЛАУ) с треугольной матрицей.

Также уровни BLAS можно интерпретировать как количество циклов в алгоритме или степень переиспользования данных.

BLAS был реализован на Fortran 77 в 80-х годах изначально для интеграции в высокоуровневые библиотеки, его референсную реализацию можно найти в Netlib-BLAS.

Также был разработан схожий API для языка С, который называется CBLAS. В целом, BLAS и CBLAS эквивалентны друг другу и используют одну и ту же реализацию алгоритмов.

LAPACK (Linear Algebra PACKage), которая написана на языке Fortran с использованием BLAS, содержит более тысячи функций, реализующих численные методы линейной алгебры. Их можно разделить на три группы:

• Решение систем линейных уравнений, включающие в себя задачи факторизации и уравновешивания матриц, поиска обратных матриц, оценки чисел обусловленности и уточнения решений СЛАУ для различных типов матриц.
• Решение задач наименьших квадратов, задач на собственные и сингулярные значения, а также уравнений Сильвестра.
• Вспомогательные (Auxiliary) и служебные (Utility) функции, используемые для выполнения отдельных вычислительных подзадач.

Первые две группы разделяются на вычислительные функции (Computational Routines) и функции-драйверы (Driver Routines). Первые выполняют отдельную вычислительную задачу, в то время как вторые решают более общую проблему.

У LAPACK, как и у BLAS, есть интерфейс LAPACKE для использования в программах на С. Он представляет собой функции-обертки для LAPACK-функций. Малая часть функциональности LAPACK реализована непосредственно в OpenBLAS на языке С и содержит дополнительные оптимизации для исполнения в многопоточной среде.

Для предоставления полного множества функций LAPACK и LAPACKE в библиотеке OpenBLAS используется реализация из LAPACK-Netlib, эталонной реализации библиотеки, которая написана на языке Fortran. На случай отсутствия Fortran-компилятора в OpenBLAS предусмотрен вариант сборки с использованием сконвертированных в язык С исходных Fortran-файлов LAPACK.

Архитектура библиотеки OpenBLAS состоит из трех уровней:

• На первом уровне с внешним API реализуется BLAS/CBLAS-интерфейсы функции. Интерфейсы BLAS для языков Fortran (BLAS) и С (CBLAS) отличаются только набором аргументов и их обработкой. В зависимости от сложности реализации функции вызывается средний или нижний уровень.
• На втором вычисления разделяются на потоки, выполняется часть вычислений, которые вынесены за оптимизированное ядро. В итоге данный уровень вызывает оптимизированные ядра из первого уровня.
• На третьем уровне оптимизированных ядер находится оптимизация алгоритмов, написанная на языке ассемблера или интринсиках.

Когда мы начали работать с OpenBLAS, ожидали, что средний и верхний уровни библиотеки, которые актуальны для всех платформ, в том числе x86, хорошо протестированы. К нашему удивлению, ошибки при тестировании наблюдались не только на нижнем уровне, где находилась оптимизация под RISC-V, но и на уровнях выше.

Но обо всем по порядку: перед тем как поделиться результатами тестов, опишу условия проведенного тестирования.

Для теста мы использовали все конфигурации OpenBLAS, актуальные под RISC-V. Всего их три:

• Generic — сборка без оптимизаций под RISC-V, ядра написаны на чистом C.
• RVV 0.7.1 — cборка с оптимизациями для ядер с RISC-V Vector Extension 0.7.1.
• RVV 1.0 — cборка с оптимизациями для ядер с RISC-V Vector Extension 1.0.

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

Каждый набор оптимизаций может собираться по-разному. Мы выбрали самые значимые для нас сборки:

• С Integer 32-bit API — собирается по умолчанию.
• С Integer 64-bit API — ключ INTERFACE64=1.
• LAPACK с использованием Fortran исходных файлов — собирается по умолчанию.
• LAPACK с использованием файлов, сконвертированных из Fortran в C (f2c) — ключ NOFORTRAN=1.

В итоге нам нужно протестировать 12 сборок библиотек (три набора оптимизации в четырех сборках).

Для сборки библиотеки использовали два компилятора (точнее — кросс-компилятора) для RISC-V:

• Для Generic и RVV 1.0 — Clang из пакета средств разработки от компании Syntacore, выложенный в открытый доступ.
• Для RVV 0.7.1 — GCC от компании T-Head, дочерней компании Alibaba Group, тулчейны выложен на GitHub.

Также эти два пакета содержат gfortran для RISC-V, необходимый нам для работы с OpenBLAS.

Есть два способа запустить исполняемый файл для RISC-V: на эмуляторах RISC-V и платах.

В качестве эмулятора мы использовали QEMU в двух вариантах. Аналогично компиляторам для сборок Generic и RVV 1.0 использовали QEMU, поставляемый в пакете для разработки от Syntacore. Для сборки с RVV 0.7.1 — QEMU от T-Head, который также выложен в открытый доступ на GitHub. Эмулятор QEMU поддерживает отладку программ — можно не только запускать, но и отлаживать локально программы для архитектуры RISC-V без реального железа!

С железом на RISC-V ситуация немного сложнее. На момент работы с OpenBLAS плат с поддержкой RVV 1.0 в открытой продаже не было — пришлось работать с эмулятором QEMU от Syntacore. В случае c Generic и RVV 0.7.1 нас снова выручила компания T-Head со своими ядрами XuanTie C906 и C910. Они установлены в платы MangoPi (одно ядро C906) и Lichee Pi 4A (четыре ядра C910).

Запуск тестов BLAS — первый этап в проверке состояния библиотеки для RISC-V. Тестировать LAPACK перед ним бессмысленно, поскольку библиотека, как говорилось ранее, полностью опирается на BLAS.

В OpenBLAS есть четыре набора тестов для BLAS/CBLAS API:

• Тесты для BLAS API. Перешли из Netlib-BLAS по наследству. Написаны на Fortran.
• Тесты для CBLAS API. Аналогичны BLAS-тестам, только CBLAS. Изначально написаны на Fortran, но также есть версия, сконвертированная в С.
• Utest. OpenBLAS-фреймворк для тестирования. Очень похож на GoogleTest, но с минимальной функциональностью. В фреймворке можно найти специфические сценарии тестирования — например, тесты на специальные значения аргументов или репродьюсеры багов.
• BLAS-Tester. Отдельный проект для тестирования BLAS-алгоритмов. В нем берется эталонная реализация алгоритмов и сравнивается с openBLAS. Тестовый исполняемый файл предусматривает возможность передачи параметров тестирования (размеры, инкременты, шаг перебора тестируемых параметров и т. д.), что позволяет значительно увеличить покрытие исходного кода тестами.

Мы будем использовать их все.

Переходим к результатам тестирования.

Generic-сборка. Все тесты в порядке. Реализации протестированы множество раз и не отличаются от других платформ.

Сборка RVV 0.7.1. Начали тестировать конфигурацию и сразу столкнулись с проблемами:

Вывод ошибки в тестах utest:

Из логов становится ясно, что проблемы в трех функциях:

• dsdot () — функция скалярного произведения, отличающаяся от привычных sdot и ddot типом данных: она принимает векторы из float, а возвращает результат в double. Погрешность совсем небольшая, что говорит о проблемах с точностью.
• *swap () и *rot () — две функции с одним тестируемым случаем: когда инкремент векторов (inc) равен нулю. Ожидается, что в этом случае будет обработан только первый элемент вектора.

Разобрались в причинах падений. Оказалось, что для dsdot () просто не было реализации, и вместо нее использовалась sdot (). Пришлось добавить. Для *swap () и *rot (), в свою очередь, нужно было добавить обработку краевого случая, приводящего к падению.

BLAS-тесты выявили проблему с функциями i[s/d]max () и i[s/d]min (), считающими индекс максимального или минимального элемента в векторе.

Вывод ошибки в тестах BLAS:

Здесь проблема была в оптимизациях, где следствием ошибки в приведении типов могло стать переполнение. Также исправили.

Тестирование библиотеки с помощью BLAS-Tester также выдало ошибку. Здесь и вовсе вылезла ошибка сегментации — SEGFAULT.

Вывод ошибки в тестах BLAS-Tester:

Исправили ошибку и взялись за следующую конфигурацию.

Сборка RVV 1.0. Здесь тесты также выявили ряд проблем. Так, в BLAS-тестах есть падения на SYMV — функции умножения симметричной матрицы на вектор.

Вывод ошибки в BLAS-тестах:

Очередная ошибка в вычислениях в оптимизированном ядре, когда в цикле аккумулятор «занулялся» на каждой итерации и накопление не происходило. Эту ошибку также исправили.

После исправления реализации BLAS-части мы взялись за LAPACK. Ожидали, что при пройденных и поправленных BLAS-тестах с этой библиотекой проблем быть не должно, ведь она полностью основана на BLAS, а исходники LAPACK хорошо протестированы.

Перед тем как показать результаты тестирования, немного расскажу про тестовую систему LAPACK. Про нее не так много информации в интернете, поэтому этот блок может быть полезен.

В тестовой системе функции делятся на 17 групп, соответствующих одному исполняемому файлу для одного из четырех типов данных: Float, Double, Complex Float и Complex Double.

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

Высокоуровневой документации для тестовой системы LAPACK нет, но исходные файлы содержат множество комментариев к тестовым функциям с подробным описанием аргументов и сценариев работы. Это очень помогает в процессе отладки тестов, так как тестовая система LAPACK основывается на математических свойствах и не использует никакие «эталонные» реализации. Часто входные данные тестируемого сценария могут зависеть от предыдущего теста, и при ошибке в начале тестирования последующие тесты могут быть заведомо ошибочными. Причем тестовая система это не отобразит, и причину ошибок нужно будет анализировать «вручную».

Чтобы понять, как работает тестовая система LAPACK, рассмотрим процесс тестирования функций, связанных с решением системы линейных уравнений с помощью LU-разложения:

1. Генерация исходной матрицы A.
2. Тестирование SGETRF () — вычисление LU-разложения с помощью функции LU.
   — Проверка выражения A = L*U.
3. Тестирование SGETRI () — вычисление обратной матрицы A^-1 с использованием LU-разложения.
   — Проверка выражения A*A^-1= I, где I — единичная матрица.
4. Генерация векторов решений X_act и вычисление векторов B = A*X.
5. Тестирование SGETRS () — решение СЛАУ A*X_comp = B.
   — Проверка X_comp = X_act.
6. Тестирование SGERFS () — итеративное уточнение X_comp.
   — Проверка X_comp = X_act.
   — Проверка численной устойчивости.
7. Тестирование SGECON () — оценка обратной величины числа обусловленности матрицы A, используя LU-разложение.
   — Сравнение результата с заранее вычисленным.

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

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

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

Далее в тексте я покажу результаты запуска тестов в табличном виде. Они представлены в формате х / y (z %), где:

• x — количество упавших тестов,
• y — количество тестов, которые были запущены; в таблице выше указаны референсные значения. Если значение y сильно ниже этого показателя, значит, многие тексты просто не запустились — и это говорит о проблеме.
• z — процент, который определяет отношение пройденных тестов ко всем тестам.

Первые запуски на выбранных конфигурациях

Тестирование первых трех сборок прошло достаточно успешно — если падения и есть, то число незначительно. Есть очевидные проблемы со сборкой NoFortran с интерфейсом Int64. Проблема не столько в упавших тестах, сколько в их количестве. На миллион тестов в итоге запустилось меньше проверок, чем должно было.

Аналогично запустим тесты для конфигураций с оптимизациями под RVV 0.7.1.

Здесь результаты запуска тестов еще хуже: около 75% не запускаются, а все тесты при использовании сконвертированных файлов для Double зависают.

Что ж, по сравнению с этой конфигурацией остальные запуски были еще вполне успешными. Половина тестов просто не дала данных, а те, что в итоге реализовались, показывают большой процент ошибок. Мы даже не смогли подвести саммари по этим тестам.

Итог: о функциональной корректности LAPACK с RVV-оптимизациями в текущем виде не может быть и речи. Необходимо разбираться в причинах падений.

Логично начать с исправления базовой сборки. Статистика падения тестов показала, что проблемы — в C-сконвертированных файлах c шириной Int 64 байта.

Для наглядности работы тестовой системы LAPACK рассмотрим сценарий с выводом результатов для тестирования функций поиска собственных чисел и векторов в случае действительных несимметричных матриц.

Пример вывода тестов:

Давайте разберем этот вывод.

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

Детальнее о способах можно почитать в комментариях тестовой функции. Также выводятся тестовые сценарии:

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

Также в выводе есть список тестов:

Далее идут сообщения о падениях тестов. Например:

Где:

• Matrix order — размер матрицы.
• type — способ заполнения тестируемой матрицы.
• seed — сид для случайной генерации матриц.
• result 11 … — номер упавшего теста.
• … is 3.193E+05 — погрешность вычисления, которая должна быть меньше порога, заданного во входном файле тестирования.

В завершение идет общая статистика по упавшим тестам. Общее число тестов: 1764 = 21 (типов матриц на входе) х 14 (типов тестов) х 6 (количество размеров матрицы).

Такой набор тестов запускается пять раз для каждого набора параметров окружения LAPACK — например, оптимальный размер блоков для блочных алгоритмов. В итоге умножаем 1764 на 5. Итого получается 8 820 тестов.

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

После длительных часов отладки тестов выяснилось, что сконвертированные С-файлы содержат две проблемы, связанные с взаимодействием тестовой системы на Fortran и сконвертированной реализации LAPACK на C.

Первая проблема состоит в том, что размер типа logical в C-файле равен 32 битам. В то время как при сборке Fortran-кода используется ключ -fdefault-integer-8, который изменяет размер не только INTEGER типа, но и LOGICAL.

Неудивительно, что не было ошибки компиляции. У Fortran все аргументы — это указатели, даже на скалярные значения. Например, API функции STREVC для вычисления правых и/или левых собственных векторов вещественной верхней квазитреугольной матрицы на Fortran выглядит так:

После конвертации реализации на язык С она выглядит следующим образом:

Размер blasint зависит от параметров сборки. В нашем случае он равен 8 байтам. Соответственно и INTEGER, и LOGICAL должны быть равны 8 байтам.

После исправления первой проблемы остались ошибки с неверными статусами для некоторых функций. Например, проблема в ошибке конвертации файла, в котором функция xerbla, реализованная на Fortran, в тестовой системе вызывается с неверным количеством аргументов:

Проблема не повсеместная — встречается всего пару раз. Возникает очевидный вопрос: «Почему неверное количество аргументов?». Дело в том, что по GCC Fortran ABI последним аргументом неявно передается размер строки. Аналогично первым аргументом передается указатель this в методах C++ классов. В случае вызова из Fortran-кода нам бы не пришлось об этом думать, но в нашем случае это необходимо указывать явно.

После всех исправлений посмотрим на результаты повторного запуска тестов (Упавшие/Запущенные).

После исправления Generic-сборки запустим тесты еще раз и посмотрим, с чем нам предстоит иметь дело (Упавшие/Запущенные):

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

В отличии от Generic-версии, тут причины падений — исключительно в оптимизациях BLAS-функций.

Первая из них — SEGFAULT в функции поиска индекса максимального по модулю элемента в векторе idamax (). Рассмотрим часть реализации алгоритма с помощью интринсиков с обработкой основной части вектора (есть еще обработка хвоста, длина которого меньше длины векторного регистра) при incx = 1. Этот файл оптимизации используется для сборки оптимизации для Float и Double, поэтому вместо интринсиков тут макросы.

Если не обращать внимание на оптимальность кода, может показаться, что все хорошо. Но есть один крайний случай, который не смогли выявить BLAS-тесты, но он проявился в тестировании LAPACK. Инструкция VMFIRSTM (для Float разворачивается в __riscv_vfirst_m_b4) может возвращать -1 в случае, если маска состоит из нулей. Что и происходит в случае передачи вектора в idamax, содержащий одни нули.

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

Проблема — в появлении статусов NaN во время выполнения LAPACK-тестов, что приводит к их зависаниям.

На самом деле исправление проблемы выглядит не так страшно, как может показаться. В самом начале ассемблерной вставки выполняется инструкция инициализации регистра нулем fmv.w.x ft11, zero. Генерация NaN происходит именно в этой строчке, так как инструкция fmv.w.x рассчитана на float-элементы (а у нас функция DGEMM). После ее замены на fmv.d.x, NaN пропадают.

Причина этого странного, на первый взгляд, поведения кроется в спецификации архитектуры RISC-V. Вместо того, чтобы оставить старшие биты регистра нетронутыми, она специально заполняет ее единицами.

После всех изменений посмотрим на окончательный результат (Упавшие/Запущенные):

Аналогичным образом — с помощью отладки и терпения — исправим зависания для конфигурации RVV 1.0.

Главная проблема оптимизации в этом случае — реализация уже известной нам функции nrm2. На первый взгляд, ничего сложного в ней нет: нужно сложить квадраты всех элементов вектора и взять из них корень. Это и предлагает нам оптимизация для RVV 1.0:

К сожалению, такая реализация работать не будет. Если в векторе содержится множество слишком малых или больших чисел, то при вычислении промежуточного результата может возникать overflow или underflow. Это приведет к потере точности вычислений или к Inf в результате.

Алгоритм с динамическим масштабированием суммы квадратов вектора решает эту проблему. В виде псевдокода он выглядит так:

В результате в сравнении с наивной реализацией существенно увеличивается количество операций и появляются операции деления чисел с плавающей точкой, которых при оптимизации все стараются избегать. Но мы избавляемся от проблемы с underflow и overflow, не возводя в квадрат самое большое число.

Например, для вектора из трех элементов порядок вычислений будет следующим:

Подробнее про вычисление нормы вектора можно прочитать здесь.

Также была проблема в функции [s/d]scal (умножение вектора на скаляр) — при умножении вектора с NaN на 0 в результате содержался NaN. От BLAS API ожидается, что результирующий вектор будет нулевым без NaN. Но тут математика и правило умножения на ноль стоит выше особенности вычислений с плавающей точкой. Проблема решилась правками в оптимизации.

Финальная статистика после всех исправлений (Упавшие/Запущенные):

Это была довольно поучительная история о том, что, даже если у open source-инструмента серьезное и обширное комьюнити, не будет лишним сначала протестировать его в своих задачах. Вас могут ожидать сюрпризы.

Над улучшением качества библиотеки OpenBLAS работали три человека, включая меня. Все исправления, которые мы сделали, выложили в open source — надеемся, это привлечет дополнительное внимание к оптимизации библиотек под RISC-V.

Ссылки на предложенные исправления

• https://github.com/OpenMathLib/OpenBLAS/pull/4452
• https://github.com/OpenMathLib/OpenBLAS/pull/4601
• https://github.com/OpenMathLib/OpenBLAS/pull/4454
• https://github.com/OpenMathLib/OpenBLAS/pull/4456
• https://github.com/OpenMathLib/OpenBLAS/pull/4499

Внимательный читатель заметит, что некоторых озвученных в тексте исправлений нет, — действительно, некоторые проблемы успели исправить до наших пул-реквестов. Помимо исправлений, про которые я рассказал в статье, мы также работали над повышением тестового покрытия OpenBLAS. Это позволило нам найти еще больше проблем, которые мы также выложили в open source.