Открытая архитектура RISC-V активно развивается: в стандарт добавляются новые расширения и инструкции, разрабатываются новые ядра и SoC. Поскольку многие компании видят перспективы архитектуры и готовы использовать ее в продакшене, создается программный стек для высокопроизводительных вычислений — RISC-V HPC (High Performance Computing). Прогресс сопровождает формирование нового тренда — OpenHPC. Он заключается в технологической независимости от решений коммерческих компаний. Причем это относится не только к ПО, но и к железу. 

Чтобы концепция OpenHPC реализовывалась быстрее, нужно, чтобы к инициативе присоединилось как можно больше компаний, помогающих в развитии экосистемы решений для RISC-V HPC. Меня зовут Андрей Соколов, я инженер-программист в компании YADRO. В R&D-команде мы поставили перед собой задачу: изучить, как можно поддержать архитектуру RISC-V со стороны библиотек линейной алгебры BLAS и LAPACK. Тестирование одной из open source-библиотек привело нас к интересным открытиям, о которых я расскажу в статье. 

HPC на 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.

Почему OpenBLAS

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

• Она полностью реализует интерфейсы 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 и 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). Первые выполняют отдельную вычислительную задачу, в то время как вторые решают более общую проблему.

Например, решение системы линейных уравнений в виде квадратной матрицы с помощью функции-драйвера SGESV — это последовательный вызов двух вычислительных функций: SGETRF, вычисление LU-разложения, и SGETRS, решение СЛАУ с использованием LU-разложения.

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

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

Архитектура OpenBLAS

Архитектура библиотеки 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 сборок библиотек (три набора оптимизации в четырех сборках).

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

Для сборки библиотеки использовали два компилятора (точнее — кросс-компилятора) для 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 без реального железа!

Ранее мои коллеги из YADRO писали о том, как можно симулировать работу и железные части процессора, чтобы протестировать и отладить программное обеспечение для чипа, которого еще не существует. Один из вариантов — как раз эмулятор QEMU, однако с некоторыми доработками.

С железом на 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/CBLAS-тестов

Запуск тестов 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-тестах →

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

Тестирование LAPACK

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

Перед тем как показать результаты тестирования, немного расскажу про тестовую систему 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 —  процент, который определяет отношение пройденных тестов ко всем тестам.

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

Generic-сборка

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

Сборка RVV 0.7.1

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

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

Сборка RVV 1.0

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

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

Исправление проблем в LAPACK

Generic-сборка

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

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

 SHS -- Real Non-symmetric eigenvalue problem
 Matrix types (see xCHKHS for details): 

 Special Matrices:
  1=Zero matrix.                          5=Diagonal: geometr. spaced entries.
  2=Identity matrix.                      6=Diagonal: clustered entries.
  3=Transposed Jordan block.              7=Diagonal: large, evenly spaced.
  4=Diagonal: evenly spaced entries.      8=Diagonal: small, evenly spaced.
 Dense, Non-Symmetric Matrices:
  9=Well-cond., evenly spaced eigenvals. 14=Ill-cond., geomet. spaced eigenals.
 10=Well-cond., geom. spaced eigenvals.  15=Ill-conditioned, clustered e.vals.
 11=Well-conditioned, clustered e.vals.  16=Ill-cond., random complex pairs 
 12=Well-cond., random complex pairs     17=Ill-cond., large rand. complx prs.
 13=Ill-conditioned, evenly spaced.      18=Ill-cond., small rand. complx prs.
 19=Matrix with random O(1) entries.     21=Matrix with small random entries.
 20=Matrix with large random entries.   

 Tests performed:   (H is Hessenberg, T is Schur, U and Z are orthogonal,
                    '=transpose, W is a diagonal matrix of eigenvalues,
                    L and R are the left and right eigenvector matrices)
  1 = | A - U H U' | / ( |A| n ulp )           2 = | I - U U' | / ( n ulp )
  3 = | H - Z T Z' | / ( |H| n ulp )           4 = | I - Z Z' | / ( n ulp )
  5 = | A - UZ T (UZ)' | / ( |A| n ulp )       6 = | I - UZ (UZ)' | / ( n ulp )
  7 = | T(e.vects.) - T(no e.vects.) | / ( |T| ulp )
  8 = | W(e.vects.) - W(no e.vects.) | / ( |W| ulp )
  9 = | TR - RW | / ( |T| |R| ulp )      10 = | LT - WL | / ( |T| |L| ulp )
 11= |HX - XW| / (|H| |X| ulp)  (inv.it) 12= |YH - WY| / (|H| |Y| ulp)  (inv.it)
 Matrix order=    2, type= 3, seed=2332,1865,2657,1485, result  11 is 3.193E+05
 Matrix order=    2, type= 3, seed=2332,1865,2657,1485, result  12 is 1.789E+05
 Matrix order=    2, type= 3, seed=2332,1865,2657,1485, result  13 is 3.193E+05
 Matrix order=    2, type= 3, seed=2332,1865,2657,1485, result  14 is 1.789E+05
 SCHKHS: Selected Right Eigenvectors from STREVC do not match other eigenvectors          N=     2, JTYPE=     4, ISEED=( 2332, 1865, 2657, 1485)
 SCHKHS: Selected Left Eigenvectors from STREVC do not match other eigenvectors          N=     2, JTYPE=     4, ISEED=( 2332, 1865, 2657, 1485)
 SCHKHS: Right Eigenvectors from SHSEIN incorrectly normalized.
 Bits of error= 0.315E+07,         N=     2, JTYPE=     4, ISEED=( 2332, 1865, 2657, 1485)
 SCHKHS: Left Eigenvectors from SHSEIN incorrectly normalized.
 Bits of error= 0.202E+07,         N=     2, JTYPE=     4, ISEED=( 2332, 1865, 2657, 1485)
 Matrix order=    2, type= 4, seed=2332,1865,2657,1485, result  11 is 6.386E+05
 // 
 // ...
 // 
 Matrix order=   16, type=19, seed=4009,2844,2856,2805, result  14 is 1.425E+05
 SCHKHS: Selected Right Eigenvectors from STREVC do not match other eigenvectors          N=    16, JTYPE=    20, ISEED=( 3866,  318, 1530, 1781)
 SCHKHS: Selected Left Eigenvectors from STREVC do not match other eigenvectors          N=    16, JTYPE=    20, ISEED=( 3866,  318, 1530, 1781)
 SCHKHS: Right Eigenvectors from SHSEIN incorrectly normalized.
 Bits of error= 0.976E+05,         N=    16, JTYPE=    20, ISEED=( 3866,  318, 1530, 1781)
 SCHKHS: Left Eigenvectors from SHSEIN incorrectly normalized.
 Bits of error= 0.220E+06,         N=    16, JTYPE=    20, ISEED=( 3866,  318, 1530, 1781)
 Matrix order=   16, type=20, seed=3866, 318,1530,1781, result  11 is 1.160E+05
 Matrix order=   16, type=20, seed=3866, 318,1530,1781, result  12 is 1.291E+05
 Matrix order=   16, type=20, seed=3866, 318,1530,1781, result  13 is 9.886E+04
 Matrix order=   16, type=20, seed=3866, 318,1530,1781, result  14 is 1.253E+05
 SCHKHS: Selected Right Eigenvectors from STREVC do not match other eigenvectors          N=    16, JTYPE=    21, ISEED=( 3575, 1284, 3532,  757)
 SCHKHS: Selected Left Eigenvectors from STREVC do not match other eigenvectors          N=    16, JTYPE=    21, ISEED=( 3575, 1284, 3532,  757)
 SCHKHS: Right Eigenvectors from SHSEIN incorrectly normalized.
 Bits of error= 0.230E+06,         N=    16, JTYPE=    21, ISEED=( 3575, 1284, 3532,  757)
 SCHKHS: Left Eigenvectors from SHSEIN incorrectly normalized.
 Bits of error= 0.249E+06,         N=    16, JTYPE=    21, ISEED=( 3575, 1284, 3532,  757)
 Matrix order=   16, type=21, seed=3575,1284,3532, 757, result  11 is 1.217E+05
 Matrix order=   16, type=21, seed=3575,1284,3532, 757, result  12 is 1.231E+05
 Matrix order=   16, type=21, seed=3575,1284,3532, 757, result  13 is 1.045E+05
 Matrix order=   16, type=21, seed=3575,1284,3532, 757, result  14 is 9.467E+04
 SHS:  332 out of  1764 tests failed to pass the threshold

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

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

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

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

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

*>    (1)     | A - U H U**T | / ( |A| n ulp )
*>
*>    (2)     | I - UU**T | / ( n ulp )
*>
*>    (3)     | H - Z T Z**T | / ( |H| n ulp )
*>
*>    (4)     | I - ZZ**T | / ( n ulp )
*>
*>    (5)     | A - UZ H (UZ)**T | / ( |A| n ulp )
*>
*>    (6)     | I - UZ (UZ)**T | / ( n ulp )
*>
*>    (7)     | T(Z computed) - T(Z not computed) | / ( |T| ulp )
*>
*>    (8)     | W(Z computed) - W(Z not computed) | / ( |W| ulp )
*>
*>    (9)     | TR - RW | / ( |T| |R| ulp )
*>
*>    (10)    | L**H T - W**H L | / ( |T| |L| ulp )
*>
*>    (11)    | HX - XW | / ( |H| |X| ulp )
*>
*>    (12)    | Y**H H - W**H Y | / ( |H| |Y| ulp )
*>
*>    (13)    | AX - XW | / ( |A| |X| ulp )
*>
*>    (14)    | Y**H A - W**H Y | / ( |A| |Y| ulp )

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

Matrix order=    2, type= 3, seed=2332,1865,2657,1485, result  11 is 3.193E+05

Где:

• Matrix order — размер матрицы.

• type — способ заполнения тестируемой матрицы.

• seed — сид для случайной генерации матриц.

• result 11 ... — номер упавшего теста.

• ... is 3.193E+05 — погрешность вычисления, которая должна быть меньше порога, заданного во входном файле тестирования.

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

SHS: 332 out of 1764 tests failed to pass the threshold

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

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

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

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

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

STREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,

*                          LDVR, MM, M, WORK, INFO )

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

typedef blasint integer;

typedef unsigned int uinteger;
typedef char *address;
typedef short int shortint;
typedef float real;
typedef double doublereal;
// ...
typedef int logical;
// ...
int strevc_(char *side, char *howmny, logical *select, 
	integer *n, real *t, integer *ldt, real *vl, integer *ldvl, real *vr, 
	integer *ldvr, integer *mm, integer *m, real *work, integer *info)

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

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

// Реализация на Fortran

SUBROUTINE XERBLA( SRNAME, INFO )

// Вызов функции на C

xerbla_("CGBCON", &i__1);

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

// Вызов функции на C

xerbla_("CGBCON", &i__1, (ftnlen)6);

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

Сборка RVV 0.7.1

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

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

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

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

 	// Заполнение векторов нулями 
	v_res = VFMVVF_FLOAT_M1(0, gvl);
    v_z0 = VFMVVF_FLOAT_M1(0, gvl);
	// Определение количества элементов, которые умещаются в векторный регистр (или группу регистров)
	gvl = VSETVL(n);
	// Инициализация вектора индексов максимальных элементов нулями
	v_max_index = VMVVX_UINT(0, gvl);
	// Инициализация вектора максимальных элементов -1
	// Так как у нас максимальное значение по модулю, -1 всегда меньше любого рассматриваемого числа
	v_max = VFMVVF_FLOAT(-1, gvl);
	for(i=0, j=0; i < n/gvl; i++) {
		// Загружаем элементы вектора в векторный регистр
		vx = VLEV_FLOAT(&x[j], gvl);

		// Взятие модуля элементов
		// Поиск элементов вектора ниже нуля, запись индексов элементов в маску
		mask = VMFLTVF_FLOAT(vx, 0, gvl);
		// Вычитание скаляра от элементов вектора
		// vx[i] = 0 - vx[i] по маске, содержащую только отрицательные элементы 
		vx = VFRSUBVF_MASK_FLOAT(mask, vx, vx, 0, gvl);

		// Поиск элементов vx[i] > v_max[i], запись индексов элементов в маску
		mask = VMFLTVV_FLOAT(v_max, vx, gvl);
		// Запись индексов максимальных элементов из маски в векторный регистр по маске максимальных элементов  
		v_max_index = VIDV_MASK_UINT(mask, v_max_index, gvl);
		// Сложение записанных индексов максимальных элементов с отступом обрабатываемого блока по маске максимальных элементов 
		// Таким образом индексы в v_max_index указывают на элементы в массиве, а не векторном регистре
		v_max_index = VADDVX_MASK_UINT(mask, v_max_index, v_max_index, j, gvl);

		// Обновление вектора максимальных элементов
		v_max = VFMAXVV_FLOAT(v_max, vx, gvl);
		j += gvl;
	}
	// Запись максимального элемента v_max в первый элемент v_res 
	v_res = VFREDMAXVS_FLOAT(v_res, v_max, v_z0, gvl);
	// Запись максимального элемента в скаляр
	maxf = *((FLOAT*)&v_res);
	// Поиск позиции максимального элемента, запись в маску 
	mask = VMFGEVF_FLOAT(v_max, maxf, gvl);
	// Поиск индекса позиции максимального элемента в маске
	max_index = VMFIRSTM(mask, gvl);
	// Взятие фактического индекса максимального элемента в массиве
	max_index = *((UINT_T*)&v_max_index+max_index);

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

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

Если хотите ознакомиться с тем, как выглядит код оптимизированного ядра, вот ссылка на GitHub.

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

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

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

When multiple floating-point precisions are supported, then valid values of narrower n-bit types, n < FLEN, are represented in the lower n bits of an FLEN-bit NaN value, in a process termed NaN-boxing.

Также нашли ошибку в функции nrm2 — вычислении L2-нормы вектора, которую мы также исправили.

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

Сборка RVV 1.0

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

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

// Максимальное количество элементов, которые умещаются в векторный регистр (или группу регистров)
size_t vlmax = VSETVL_MAX;
// Заполнение векторов нулями
v_res = VFMVVF_FLOAT_M1(0, vlmax);
vr = VFMVVF_FLOAT(0, vlmax);
 
if(inc_x == 1) {
    for (size_t vl; n > 0; n -= vl, x += vl) {
        // Определение количества элементов, которые умещаются в векторный регистр (или группу регистров)
        vl = VSETVL(n);
        // Загрузка элементов вектора в векторный регистр
        v0 = VLEV_FLOAT(x, vl);
        // Перемножение элементов вектора и аккумулирование суммы квадратов в vr 
        vr = VFMACCVV_FLOAT(vr, v0, v0, vl);
    }
} else {
    // Реализация аналогична inc_x == 1, только загрузка с шагом inc_x
}

// Сохранение суммы в первый элемент v_res 
v_res = VFREDSUM_FLOAT(v_res, vr, v_res, vlmax);
// Сохранение суммы квадратов в скаляр
ssq = VFMVFS_FLOAT_M1(v_res);

return sqrt(ssq);

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

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

t=0 # scale factor
s=1 # scaled sum
for i = 1:n
  if |x_i| > t
      s = s*(t/x_i)^2 + 1
      t = |x_i|
  else
    s += (x_i / t)^2
return t*sqrt(s)

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

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

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

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

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

Заключение

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

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