CC BY
22
ББК 3973:3972.1 ГРНТИ 50.09.31, 50.33.04 УДК 004.222+004.272.25
К. С. Исупов, В. С. Князьков
Матрично-векторное умножение многократной точности на графическом процессоре
Аннотлция. Мы рассматриваем параллельную реализацию матрично-вектор-ного умножения (GEMV, уровень 2 BLAS) для графических процессоров (GPU) с использованием арифметики многократной точности на основе системы остаточных классов. В нашей реализации GEMV покомпонентные операции с многоразрядными векторами и матрицами разбиваются на части, каждая из которых выполняется отдельным CUDA ядром. Это исключает ветвление логики исполнения и позволяет добиться более полного использования ресурсов GPU. Эффективная структура данных для хранения многоразрядных массивов обеспечивает объединение доступов параллельных потоков к глобальной памяти GPU в транзакции. Для предложенной реализации GEMV выполнен анализ ошибок округления и получены оценки точности. Представлены экспериментальные результаты, показывающие высокую эффективность разработанной реализации по сравнению с существующими программными пакетами многократной точности для GPU.
Ключевые слова и фразы: вычисления высокой точности, BLAS, GEMV, параллельные
алгоритмы, CUDA, GPU, система остаточных классов.
Введение
Операции с плавающей точкой имеют ошибки округления, возникающие непосредственно во время вычислений. Такие ошибки являются естественными из-за ограниченной длины мантиссы в стандартных форматах IEEE 754 одинарной и двойной точности (binary32 и binary64, соответственно). Для многих задач эти ошибки не препятствуют получению корректного результата вычислений. Более того, для ряда приложений, таких как машинное обучение, лучшим вариантом является использование форматов пониженной (половинной) точности [1]. Вместе с тем, в настоящее время возникает все больше критичных
Исследование выполнено за счет гранта Российского научного фонда (проект № 18-71-00063).
© К. С. Исупов( , В. С. Князьков«, 2020 © Вятский государственный университет« , 2020 © Пензенский государственный университет« , 2020
© Программные системы: теория и приложения (дизайн), 2020
к ошибкам округления задач, для которых стандартных форматов IEEE 754 оказывается недостаточно [2-7]. Такие задачи решаются с использованием арифметических библиотек многократной точности, позволяющих выполнять операции с числами, разрядность которых превышает стандартные форматы, а в общем случае ограничена лишь объемом доступной памяти вычислительной системы.
К числу широко известных библиотек многократной точности относятся GMP, MPFR, MPIR и MPDECIMAL. Эти и многие другие высокоточные библиотеки ориентированы на системы с центральными процессорами (CPU). Вместе с тем приложения, требующие вычислений повышенной точности, зачастую являются крайне ресурсоемкими и могут быть ускорены за счет использования аппаратных платформ с массивно-параллельной архитектурой, таких как графические процессоры видеокарт (GPU).
Современный графический процессор представляет собой массив потоковых мультипроцессоров, каждый из которых содержит в своем составе множество потоковых процессоров. Массивный параллелизм GPU обеспечивается репликацией общей архитектуры мультипроцессора, причем каждый потоковый процессор в один момент времени может выполнять одну и ту же инструкцию над множеством данных. В 2007 году компания NVIDIA выпустила архитектуру параллельных вычислений CUDA (Compute Unified Device Architecture). Используя CUDA, программисты получили возможность разработки параллельных алгоритмов для GPU на специальном диалекте языка C, без необходимости сложного графического программирования [8]. Производительность последних версий GPU, таких как NVIDIA Tesla V100, сравнима с производительностью 100 CPU, что делает GPU эффективным инструментом для ресурсоемких вычислений.
В статье мы рассматриваем параллельную реализацию матрично-векторного умножения (GEMV) многократной точности для CUDA-совместимых GPU. GEMV входит в состав уровня 2 базовых подпрограмм линейной алгебры (BLAS) [9] и выполняет одну из следующих операций:
(1) y ^ aAx + ßy или y ^ aATx + ßy.
где A —плотная M x N матрица, x и y — векторы, а и ß — скаляры. Наряду с другими BLAS операциями, GEMV является важным строительным блоком для многих алгоритмов линейной алгебры, таких как решатели линейных систем и задач на собственные значения.
Одной из основных стратегий оптимизации вычислений линейной алгебры на GPU является блочное разбиение (blocking), в соответствии с которым матрица или вектор разбивается на небольшие части, которые однократно загружаются в разделяемую память GPU и затем многократно используются в арифметических операциях [10]. Целью является снижение числа обращений к глобальной памяти. Однако при работе с многоразрядными числами интенсивное использование разделяемой памяти (равно как и регистров) часто оказывается неэффективным, так как при большом размере каждого числа снижается масштабируемость и быстродействие CUDA ядер. Поэтому в нашей реализации GEMV используется другой подход, предложенный в [11]. Он основан на разбиении арифметических операций многократной точности на несколько этапов, каждый из которых выполняется отдельным CUDA ядром со своей конфигурацией, причем все цифры многоразрядных чисел вычисляются параллельно. Такой подход приводит к некоторому увеличению числа обращений к глобальной памяти, но вместе с тем обеспечивает высокую производительность и хорошую масштабируемость вычислений многократной точности на GPU по сравнению с традиционной парадигмой, в соответствии с которой каждая арифметическая операция с многоразрядными числами выполняется одним потоком. Для реализации этого подхода мы используем систему остаточных классов (СОК) [12].
В СОК число представляется своими остатками от деления на заданный набор модулей, называемый базисом. Разрядность базиса определяется суммой разрядностей всех модулей [13]. Остатки являются взаимно независимыми, а для операций сложения, вычитания и умножения вычисления с каждым остатком выполняются в кольце целых чисел по соответствующему модулю, так что цепочки переносов между остатками исключаются, как показано на рисунке 1. Это позволяет вычислять все остатки параллельно.
Многоразрядные операнды
Преобра-
зование
2СС
СОК
Вычисления по модулю mi (+,-*)
Вычисления по модулю т2 (+,-,*)
Вычисления по модулю т„ (+,-*)
Преобра-
зование
—>
СОК
2СС
>
Многоразрядный результат
J
Рисунок 1. Параллельные арифметические свойства СОК
Оставшаяся часть статьи построена следующим образом. В разделе 1 приведен краткий обзор связанных работ. В разделе 2 описан используемый формат чисел многократной точности. В разделе 3 приводится структура данных для хранения многоразрядных массивов и используемая схема размещения матрицы в памяти GPU. В разделе 4 рассматриваются алгоритмы и основные особенности предложенной реализации GEMV, а в разделе 5 оценивается ее точность. Результаты экспериментов представлены в разделе 6, а последний раздел содержит выводы и направления дальнейших исследований.
1. Высокоточные вычисления и BLAS для GPU
Операции умножения плотной матрицы на вектор (GEMV, SYMV, TRMV) часто возникают в приложениях научных вычислений, и в литературе рассматриваются эффективные реализации этих операций на GPU с использованием арифметики IEEE 754 [10,14-16]. В то же время существуют работы, посвященные созданию BLAS реализаций повышенной точности с поддержкой GPU [11,17-19].
Распространенным в настоящее время способом повышения точности вычислений является использование разложений с плавающей точкой (floating-point expansions): число расширенной точности представляется в виде невычисленной суммы нескольких стандартных чисел с плавающей точкой. Примером таких разложений является известный формат double-double, в котором каждое число представляется в виде невычисленной суммы двух чисел двойной точности (binary64), что соответствует четверной точности (106 бит мантиссы). В свою очередь, формат quad-double обеспечивает восьмерную точность (212 бит мантиссы) за счет использования четырех чисел в формате binary64 для представления каждого числа расширенной точности [20]. Алгоритмы вычисления разложений с плавающей точкой называются безошибочными преобразованиями (error-free transformations) [21,22]. Арифметика double-double используется в нескольких программных пакетах, например в QPBLAS-GPU [17].
В работе [18] представлен ExBLAS — пакет оптимизированных операций линейной алгебры, обеспечивающий воспроизводимость результатов вычислений. Под воспроизводимостью понимается возможность получения идентичных результатов при многократных запусках программы с одними и теми же входными данными. В ExBLAS для обеспечения воспроизводимости используются безошибочные
преобразования и специальные аккумуляторы — представления с фиксированной точкой большой длины, способные хранить каждый бит информации стандартного формата с плавающей точкой (binary64). Использование аккумуляторов позволяет заменить неассоциативные операции с плавающей точкой операциями с фиксированной точкой, являющимися ассоциативными. ExBLAS предоставляет реализации ряда воспроизводимых операций линейной алгебры для CPU, сопроцессоров Intel Xeon Phi, а также GPU компаний NVIDIA и AMD.
В недавнем исследовании [19] представлены оптимизированные CUDA реализации операций DOT, GEMV, GEMM и SpMV, составляющие пакет BLAS-DOT2. В этих реализациях входные данные хранятся в стандартных форматах binary32 и binary64, а для внутренних вычислений используется двукратная точность: для входных данных в формате binary32 вычисления выполняются в формате binary64; в свою очередь, для входных данных в формате binary64 вычисления выполняются с использованием алгоритма Dot2 [23], основанного на безошибочных преобразованиях.
Также существует ряд арифметических библиотек расширенной или многократной точности для GPU. Поддержка форматов double-double и quad-double реализована в библиотеке GQD [24], которая, помимо основных арифметических операций, позволяет вычислять ряд часто используемых математических функций расширенной точности, таких как квадратный корень, логарифм, экспонента и тригонометрические функции. Алгоритмы GQD по большей части основаны на алгоритмах, реализованных в известной CPU-библиотеке QD1. Для хранения чисел расширенной точности в GQD используются предоставляемые CUDA векторные типы double2 и double4.
В библиотеке CAMPARY [25] используются n-компонентные разложения с плавающей точкой. Данная библиотека предоставляет гибкие CPU и GPU реализации основных арифметических операций многократной точности. В качестве базовых типов для отдельных компонент разложения поддерживаются форматы binary32 и binary64, а точность вычислений (количество компонент в разложении) задается шаблонными параметрами функций. Каждая операция сложения и умножения n-компонентных разложений в CAMPARY в общем случае требует выполнения, соответственно, 3n2 + 10n — 4 и 2n3 + 2n2 + 6n — 4 стандартных операций с плавающей точкой, а для double-double арифметики (при использовании двухкомпонентных разложений) применяются оптимизированные алгоритмы.
1https://www.davidhbailey.com/dhbsoftware
Библиотеки GARPREC [24] и CUMP [26] поддерживают произвольную точность на GPU с использованием формата "multi-digit". Число в таком формате представляется в виде последовательности цифр (целых чисел машинной точности) с одной общей экспонентой. Алгоритмы GARPREC аналогичны алгоритмам, реализованным в известном пакете ARPREC1 для CPU. CUMP, в свою очередь, основана на библиотеке GMP (GNU MP Bignum Library)2. Большинство функций CUMP имеют интерфейс, подобный интерфейсу GMP. В GARPREC и CUMP каждая высокоточная операция реализована в виде одного потока, а для объединения доступов к памяти GPU в транзакции используется схема интервальной индексации.
В [11] авторы данной статьи предложили новые алгоритмы арифметики многократной точности на основе СОК, а также оптимизированные алгоритмы для параллельной обработки векторов многократной точности на GPU, которые в дальнейшем были расширены и адаптированы для обработки матриц. На основе этих алгоритмов разработан новый высокоточный пакет базовых подпрограмм линейной алгебры для GPU—MPRES-BLAS. Проведенные эксперименты показали, что во многих случаях MPRES-BLAS обладает более высоким быстродействием по сравнению с реализациями, основанными на существующих позиционных библиотеках многократной точности для CPU и GPU. Рассматриваемая в данной статье реализация GEMV является частью MPRES-BLAS.
Таблица 1 содержит сводную информацию по рассмотренному программному обеспечению со ссылками на исходный код.
Таблица 1. Высокоточное программное обеспечение с поддержкой GPU
Пакет Ссылка Исходный код
QPBLAS-GPU [17] https://ccse.jaea.go.jp/software/QPBLAS-GPU
ExBLAS [18] https://github.com/riakymch/exblas
BLAS-DOT2 [19] http://www.math.twcu.ac.jp/ogita/post-k/results.html
GQD [24] https://code.google.com/archive/pZgpuprec
CAMPARY [25] http://homepages.laas.fr/mmjoldes/campary
GARPREC [24] https://code.google.com/archive/p/gpuprec
CUMP [26] https://github.com/skystar0227/CUMP
MPRES-BLAS [11] https://github.com/kisupov/mpres-blas
2https://gmplib.org
2. Представление многоразрядных чисел с плавающей точкой с использованием СОК
Система остаточных классов задается набором из n попарно взаимно простых целых чисел (модулей) {mi,m2,...,mn}. Динамический диапазон СОК определяется произведением модулей3
(2) M = mi • m-2.....mn.
Целое число X G [0, M — 1] однозначным образом представляется в СОК остатками (xi, X2,..., xn), где Xj = |X|mi соответствует операции X mod mj. Преобразование из СОК в двоичную систему выполняется с использованием китайской теоремы об остатках (КТО)[12]:
(3) X
м
У^ Мг |хгШг|т
¿=1
где Мг и тг — константы, такие что Мг = М/тг и тг — мультипликативная инверсия4 Мг по модулю тг.
В отличие от операций сложения, вычитания и умножения, которые выполняются параллельно над всеми остатками хг, операции сравнения (за исключением случая сравнения на равенство), определения переполнения, вычисления знака, масштабирования и деления, называемые "немодульными", являются трудоемкими для СОК. Классический метод выполнения немодульных операций, основанный на КТО, становится неэффективным при больших динамических диапазонах, состоящих из сотен и тысяч бит. Поэтому мы используем альтернативный метод, основанный на интервальной оценке дробного представления числа в СОК [27]. Для X, заданного остатками (ж1,ж2,... , хп), дробное представление (относительная величина) вычисляется следующим образом:
(4) X
£
|XjWj|r
Интервальная оценка для Х/М обозначается I(Х/М) = [X/М, X/М] и представляет собой такой интервал, что X/М < Х/М < X/М.
3Обычно для обозначения произведения модулей СОК используют символ М. В данной работе, чтобы избежать коллизии с количеством строк матрицы, мы используем символ М для обозначения произведения модулей СОК.
4Если х является ненулевым целым числом, то у называется мультипликативной инверсией х по модулю т, если |х X у|т = 1.
m
i
Границами этого интервала, X/М и X/М, являются числа с плавающей точкой машинной точности, и для их вычисления требуются только стандартные арифметические операции, вне зависимости от размера набора модулей СОК и величины динамического диапазона. С использованием интервальных оценок разработаны эффективные алгоритмы выполнения ряда немодульных операций в СОК, включая сравнение по величине и общее деление [27].
Мы используем следующий способ представления многоразрядных чисел с плавающей точкой:
(5) x =(-1)s
Mi \xiWi\r
i=i
x 2e
M
где s € {0,1} — знак числа, e — экспонента (порядок) и xi, x2,... ,xn — цифры (остатки) мантиссы X, представленной в СОК, которая принимает значения в диапазоне от 0 до M — 1. Каждая цифра xi = X mod mi является целым машинным числом со знаком.
В качестве дополнительной информации в числовой формат включена интервальная оценка мантиссы, I(X/M), которая позволяет реализовать эффективное сравнение, вычисление знака результата, выравнивание экспонент чисел и оценку необходимости округления.
В предыдущих работах авторов рассматривались различные способы выполнения арифметических операций с числами вида
(5). В работе [11] предложены улучшенные алгоритмы сложения и умножения многократной точности, а также получена следующая граница относительной ошибки вычислений. Если fl(xoy) — округленный результат операции o € {+, —, х} с числами вида (5), то
(6) fl(x o y) = (x o y)(1 + S), |S| < u, u < 4/VM.
В разделе 5 мы используем (6) для оценки точности разработанной операции матрично-векторного умножения.
3. Размещение данных
В нашей реализации GEMV используется column major формат хранения матрицы, который является стандартным для библиотек BLAS и LAPACK [14]. Для того, чтобы обеспечить объединение доступов к глобальной памяти GPU в транзакции (coalesced access), массив многоразрядных чисел (вектор или матрица) хранится в виде структуры mp_array_t, описание которой представлено в таблице 2.
Таблица 2. Спецификация структуры тр_аггау_ для хранения массива многоразрядных чисел; здесь п — число модулей СОК, Ь — размер массива
Поле Описание
Массив размера п х Ь, который хранит цифры многоразрядных мантисс. Все цифры одного числа хранятся последовательно в памяти.
digits
len
sign Целочисленный массив размера L, который хранит знаки чисел. о
Е
exp Целочисленный массив размера L, который хранит экспоненты чисел. s
о
ж
Массив чисел с плавающей точкой с расширенной экспонентой размера 2L, н
eval который хранит интервальные оценки мантисс. Первые L элементов хранят Е
нижние границы интервальных оценок, а вторые L — верхние. н
ог
Массив (буфер) длины L с элементами векторного типа int4 из стандартных Р
buf заголовков CUDA C/C++. Используется в операциях сложения и вычитания для т
передачи вспомогательных переменных между вычислительными ядрами. °
й
т
Количество элементов в многоразрядном массиве (L). Для вектора длины ч
N массив должен содержать не менее (1 + (N — 1) х |incx|) элементов, где с
incx определяет расстояние между двумя соседними элементами вектора; и
для матрицы размера M х N, массив должен содержать не менее LDA х N а
элементов, где LDA определяет ведущую размерность матрицы; значение LDA j-,
должно быть не менее max(1,M). U
На рисунке 2 приведен пример размещения многоразрядной матрицы 3 х 4 в памяти СРи. В этом примере п = 4, т.е. мантисса каждого многоразрядного элемента а^ состоит из четырех цифр: X = (хх, Х2, хз, х4). Символ "." используется для доступа к отдельным частям элементов. Символы 1о и ир обозначают нижнюю и верхнюю границы интервальной оценки мантиссы, соответственно, т.е. 1о := X/М и ир :=
X /M.
Матрица А
Логическое представление (column major)
ft Зц a i2 а 1з 314
321 322 323 324
a3i а32 азз Эз4
1
«- N -s
> Я
au 321 331 3l2 322 а32
3l3 323 азз а 14 324 Эз4
о
\7
\7
Многоразрядный массив (тр array t) digits:
an.xi ЭЦ.Х2 Эц.Хз Эц.Х4 a2i.Xi a2i.x2 a2i.x3 Э21-Х4 a3i.xi a3i-x2 a3i-X3 Э31-Х4
a12.X! ai2-Xz ai2.x3 a12-X4 a22.xi a22.x2 a22.x3 a22.X4 a32.xi a32.x2 a32.x3 a32.X4
a13.Xi ai3-x2 au.x3 al3.X4 a23.Xi a23.x2 a23.x3 a23.X4 a33.Xi a33.x2 a33.x3 a33.X4
ai4.Xi ai4-x2 ai4.x3 ai4.X4 a24.Xi a24.x2 a24.x3 Э24-Х4 a34-Xi a34-x2 a34-x3 a34-X4
sign:
au.s a2i.s a3i.s a12.s a22.s a32.s ai3.s a23.s a33.s ai4.s a24.s a34.s
exp:
ац.е a2i.e a3i-e ац.е a22.e a32.e ац.е a23.e a33e ац.е a24.e a34e
eval:
au.lo 321-10 a31.lo au.lo 322-10 a32.lo a13.lo a23.lo a33.lo ai4.lo a24- lo a34. lo
au.up a21.up a31.up a12.up a22.up a32.up a13.up a23.up a33.up a14.up a24.up a34.up
buf:
ац.Ь | a2i.b a3i.b a12.b a22.b a32.b ац.Ь a23.b азз-Ь a14.b a24.b аз4-Ь
len: H
Рисунок 2. Размещение многоразрядной матрицы 3 х 4 (LDA = 4) в памяти GPU.
Многоразрядные элементы
/ \ {+,-,*} ч
{+-,*}
{+,-*}
{+,-,*}
{+,-*}
|—j |—11—j j—j Вычисление знаков,
4 Kernel2 0 экспоненти
| интервальных оценок
□ □□□ _
□ □□□
II II II I Вычисление цифр
□жН.....-........-■■
мантисс в системе остаточных классов
□□□□ □гигп[]
1 Kernel 3 hd
Округление результата
I_II_II_I
□□□□
Рисунок 3. Выполнение покомпонентных операций с массивами многократной точности на GPU.
4. Алгоритмы реализации GEMV на GPU
В работе [11] авторы данной статьи предложили алгоритмы для вычислений с векторами многократной точности, ориентированные на GPU. В предложенных алгоритмах арифметические операции разбиваются на следующие части, каждая из которых выполняется как отдельное CUDA ядро:
• вычисление знаков, экспонент и интервальных оценок;
• вычисление цифр (остатков) многоразрядных мантисс в СОК;
• округление.
Такая декомпозиция (см. рисунок 3) исключает ветвления при выполнении последовательных этапов арифметических операций с многоразрядными числами. Кроме этого, для каждого вычислительного ядра задается своя конфигурация запуска (количество блоков потоков и размер блока), что позволяет эффективно задействовать все имеющиеся ресурсы GPU. Данный подход лежит в основе нашей реализации GEMV, описание которой представлено в алгоритме 1. Каждый шаг алгоритма, кроме шага 4, состоит из запуска трех CUDA ядер. Шаг 4 выполняется одним ядром, так как каждая отдельная операция сложения многократной точности на этом шаге вычисляется последовательно (без распараллеливания по модулям СОК).
Алгоритм 1. Операция GEMV многократной точности 1: Вектор х умножается на скаляр a: d ^ ах. Для хранения результатного
вектора d используется буфер в глобальной памяти GPU. 2: Вектор у умножается на скаляр в: у ^ ву.
3: Вектор d рассматривается как главная диагональ матрицы D:
где L = N для нетранспонированной матрицы и L = M для транспонированной матрицы. На данном шаге матрица A масштабируется с левой или с правой стороны (в зависимости от типа выполняемой операции) диагональной матрицей D. Масштабирование справа (B ^ AD) заключается в умножении каждого i-го столбца матрицы A на соответствующий диагональный элемент di. Масштабирование слева (B ^ DA) заключается в умножении каждой i-й строки матрицы A на соответствующий диагональный элемент di. В результате этой операции формируется промежуточная матрица B размера M х N, для хранения которой используется буфер в глобальной памяти GPU.
4: Элементы матрицы B суммируются по строкам или по столбцам, в зависимости от выполняемой операции, и вычисленный вектор сумм складывается с ву.
Таким образом, выполнение операции GEMV многократной точности на GPU заключается в запуске в общей сложности десяти CUDA ядер. Стоит отметить, что накладные расходы, связанные с этими запусками, не вносят существенного вклада в общее время вычислений даже при небольших матрицах, поскольку операции многократной точности сами по себе достаточно затратные.
Умножение векторов на скаляры на шагах 1 и 2 алгоритма 1 может быть выполнено с использованием функции SCAL (уровень 1 BLAS), параллельная CUDA реализация которой с использованием чисел вида
(5) рассмотрена в [11]. Поэтому мы более подробно остановимся лишь на шагах 3 и 4. Так как для нетранспонированной и транспонированной матрицы используются различные схемы вычислений, далее эти случаи рассматриваются по отдельности.
di 0 0 0 d2 0
(7) d =(di,d2,ds ...,db) ^ D = 0 0 d3
0 0 0
000
4.1. Нетранспонированная матрица
Схема вычисления у ^ аАх + ву представлена на рисунке 4. Каждый столбец матрицы А умножается на один элемент вектора ! = ах (т.е. выполняется масштабирование диагональной матрицей справа), в результате чего формируется соответствующий столбец матрицы В. Затем производится суммирование элементов матрицы В по строкам с прибавлением г-го элемента вектора ву к г-й вычисленной сумме.
Рисунок 4. Схема выполнения GEMV для нетранспониро-ванной матрицы A
Вычисление матрицы B
Как указано выше, вычисление матрицы B производится тремя ядрами CUDA: Kernel 1 (вычисление знаков, экспонент и интервальных оценок), Kernel 2 (вычисление цифр многоразрядных мантисс) и Kernel 3 (округление элементов вычисленной матрицы). Далее рассматриваются принципы работы каждого ядра.
Kernel 1. Вычисления выполняются на двумерной сетке блоков потоков. Размер сетки — gridDiml х gridDiml; размер каждого блока равен blockDiml. Параметры gridDiml и blockDiml являются настраиваемыми. Координата blockldx.y определяет смещение по столбцам матрицы, т.е. все блоки для которых blockIdx.y = i (индексация начинается с нуля) участвуют в формировании i-го столбца матрицы B, выполняя умножение i-го столбца матрицы A и i-го элемента вектора d = ax, который предварительно загружается в регистры или локальную память потока. Если gridDiml < N, то блоки для которых blockIdx.y = i вычисляют в цикле столбцы матрицы с индексами i, (i + gridDiml), (i + 2 х gridDiml) и т.д.
В данном CUDA ядре знак, экспонента и границы интервальной оценки одного многоразрядного числа вычисляются одним потоком. Благодаря использованию структуры mp_array_t для хранения многоразрядных массивов, доступы потоков к глобальной памяти GPU объединяются в транзакции.
Kernel 2. Вычисления выполняются на двумерной gridDim2 х gridDim2 сетке блоков потоков, где gridDim2 — настраиваемый параметр. Назначение блоков столбцам матрицы аналогично описанному выше для Kernel 1, однако в данном случае каждый поток ассоциирован со своим модулем СОК, и таким образом все n цифр многоразрядной мантиссы вычисляются одновременно n потоками в рамках одного блока, причем в каждом блоке потоков могут вычисляться одновременно несколько многоразрядных мантисс, в зависимости от соотношения между величиной набора модулей системы остаточных классов n и минимальным размером блока, необходимым для полной загрузки потокового мультипроцессора GPU.
Для данного вычислительного ядра размер блока потоков рассчитывается автоматически. Минимальный размер блока зависит от версии аппаратной совместимости GPU (Compute Capability) и определяется следующим образом [11]:
(8) MinBS = MaxThreads/MaxBlocks,
где MaxThreads и MaxBlocks — максимальное количество резидентных потоков и блоков в мультипроцессоре, соответственно. Например, если MinBS = 64 и n = 4, то в каждом блоке потоков будут одновременно вычисляться цифры для 16 многоразрядных мантисс.
Kernel 3. Данное вычислительное ядро выполняется на одномерной сетке одномерных блоков потоков, причем матрица B рассматривается как вектор длины M х N. Используемый алгоритм округления и его CUDA реализация рассмотрены в [11]. При округлении числа вида (5) необходимо заново вычислять интервальную оценку мантиссы. Актуальный алгоритм вычисления интервальной оценки представлен в [27].
Редукция матрицы B
Суммирование элементов каждой строки матрицы B выполняется в соответствии с Алгоритмом 6 из [11]. Для этого запускается M блоков потоков, и в каждом г-м блоке элементы г-й строки матрицы загружаются в разделяемую память, после чего выполняется редукция по разделяемой памяти. Максимальный размер каждого блока зависит от точности вычислений (размера модулей СОК) и ограничен доступным объемом разделяемой памяти.
Каждая операция сложения многократной точности выполняется одним потоком. Для хранения многоразрядных чисел в разделяемой памяти используется массив структур mp_float_t. В отличие от структуры mp_array_t, которая представляет вектор многоразрядных чисел и может быть инициализирована только на стороне хоста, mp_float_t представляет одно многоразрядное число и может быть инициализирована на стороне GPU. Преобразование между mp_array_t и mp_float_t осуществляется непосредственно при выполнении операции сложения и не приводит к накладным расходам.
В результате основного цикла редукции в разделяемой памяти г-го блока хранится вычисленная сумма элементов г-й строки матрицы B, к которой затем прибавляется соответствующий элемент вектора у, предварительно умноженный на скаляр в. Так как каждый блок потоков ассоциирован со своим элементом вектора у, то синхронизация между блоками не требуется.
4.2. Транспонированная матрица
Схема вычисления у ^ aATx + ву представлена на рисунке 5. В данном случае каждый столбец матрицы A покомпонентно умножается на весь вектор d = ax (иначе говоря, выполняется масштабирование
матрицы А диагональной матрицей Б слева), в результате чего формируется соответствующий столбец матрицы В. Затем производится суммирование элементов матрицы В по столбцам с прибавлением г-го элемента вектора вУ к г-й вычисленной сумме.
Blocks # (0,0), (1, 0), (2, 0)... Blocks # (0,1), (1,1), (2,1)... Blocks # (0, 2), (1,2), (2, 2)... Blocks # (0, 3), (1,3), (2, 3)...
Матрица А
Block #0 Block # 1 Block # 2 Block # 3
Матрица В
Л/
\ ! \ ! \ ! \ / \ !
Ж
Рисунок 5. Схема выполнения GEMV для транспонированной матрицы A
Также как и в первом случае, для вычисления матрицы B производится запуск трех вычислительных ядер — Kernel 1, Kernel 2 и Kernel 3, причем Kernel 1 и Kernel 2 выполняются на двумерных сетках из блоков потоков, а все блоки с индексом blockldx.y = i участвуют в формировании i-го столбца матрицы B. Схема обращения потоков к матрице остается прежней. Отличие заключается в том, что каждый поток оперирует соответствующим элементом вектора d. После того как матрица B сформирована, вычисление i-го элемента вектора результата сводится к суммированию элементов i-го столбца матрицы с прибавлением к вычисленной сумме i-го элемента вектора у, предварительно умноженного на скаляр в.
5. Оценка точности
Полагая М = N, получим оценки точности для представленной реализации СЕМУ. Пусть выполняется операция у ^ аАх + ву. Обозначим у* = (у**, у***,... ) — точный результат и у = (уг,у2,... ,уN) — результат применения операции СЕМУ. Будем оценивать границы абсолютной и относительной ошибки, ||у* — у У и ||у* — у ||/||у*||, соответственно (где || || — ^1-норма), используя стандартную модель арифметики с плавающей точкой [28, с. 40].
Для чисел вида (5) стандартной модели (в ограниченной форме) соответствует соотношение (6). Для упрощения промежуточных выражений будем использовать нотацию вп, которая определяется следующим образом [28, с. 63]:
п
(9) Ш + = 1 + °п, К .
х х 1 — пи
¿=1
Здесь 5ъ — относительная ошибка, вносимая г-й операцией с плавающей точкой из последовательности операций, входящих в исследуемое выражение, причем \5ъ\ < и. Целью нашего анализа является получение границы итоговой ошибки СЕМУ и нет необходимости различать вклад отдельных арифметических операций. Поэтому примем далее 5г = 5, \5\ < и. В соответствии с (6) имеем и < 4/%/М.
Вначале вычисляется вектор с! = ах, элементы которого в соответствии со стандартной моделью выражаются следующим образом:
(10) ! = й(ахъ) = ахъ(1 + 5) = ахъ(1 + в1), 1 < г < N.
С учетом используемых обозначений, вычисленные элементы матрицы В могут быть записаны в следующей форме:
(11) Ъщ = ¿у) = а^¿у (1 + 5) = аа^х^ (1 + в^), 1 < г,3 < N.
Далее вычисляется сумма элементов каждой строки матрицы В. Для г-й строки обозначим эту сумму зъ. Хорошо известно, что точность суммирования зависит от выбранного алгоритма. Как отмечено выше, зъ вычисляется одним блоком потоков. И хотя используемый Алгоритм 6 из [11] реализует попарное суммирование, максимальный размер блока (количество параллельных потоков, выполняющих суммирование) зависит от величины набора модулей СОК, т.е. от разрядности вычислений.
Когда размер блока равен единице, вг вычисляется последовательно, что приводит к классическому рекурсивному алгоритму: вг ^ Ъц
з = 2 N ёо яг ^ й(вг + Ъ^) епё Рог
Для этого алгоритма вычисленное значение вг может быть представлено в следующем виде [29]:
(12) 8г = (Ъг! + Ъг2)(1 + 0N -1) + Ъ«(1 + 0N-2) + Ъ<4 (1 + 0N _з) + • • •
+ Ъг№ (1+ 01),
или то же самое:
вг = (аацХ1 + ааг2Ж2)(1 + 0^+1) + ««¿3X3(1 + 0N)
+ ааг4 Х4(1 + 0N-1) +----+ аагм XN (1 + 0з)
Сложение вг и вуг(1 + 01) дает
(13) Уг = й(вг + вУг(1 + 01)) = вг(1 + 01) + вУг (1 + 02),
и после подстановки получаем финальное выражение для г-го элемента вычисленного вектора:
N
Уг = вУг + а.а,цХ^ + вУг02 + аа^X10N+2 + ааi2X20N+2 3=1
+ аaгзXз0N +1 + аaг4X40N + • • • + аaiN XN 04.
Для получения верхней границы ошибки следует исходить из того, что отдельные компоненты 0г взаимно не компенсируются. Будем считать все 0г > 0. Тогда приведенное выше выражение примет вид
N
Уг = вУг + 53 аацXj + 02 |вУг | + 0N+2 ^1X1 | + 0N+2 ^2X2 | 3=1
+ 0N+l|aaгзxз| + 0N ^«¿4X41 + • • • + 04|aaгN XN |.
Для дальнейшего упрощения перейдем к записи 7„ [28, с. 63], которая определяется следующим образом:
Пи 1
(14) 7„ := -- при пи < 1.
1 — пи
Одним из свойств этой записи является 7п < ^п+1 (см. [28, с. 67]). Учитывая это свойство, представленное выше выражение для уъ можно переписать в следующем виде:
N
(15) уъ = вуъ + аащх^ + EN,
3 = 1
где Еп < YN+2(\вуъ \ + \аазх^ \).
Поскольку г-й элемент точного вектора-результата у* выражается в форме у* = вуъ + £^=1 аа^х3-, имеем следующую верхнюю границу абсолютной ошибки:
(16) ||у* — у || < 1 {-(++и)и • £ (\в*\ + £ \аац хз
и верхнюю границу относительной ошибки:
\\у* - у\\ ^ (N + 2)u Ef=i(Ißyil + ЕN=1 laaijXj|)
N in , v^N
(17) \\y ~y" <
||у*| " 1 — ^ + 2)и £N=1 \вуг + Е;=1 аазхз \
где и < 4/%/М и М — произведение модулей СОК. То же самое в более лаконичной форме:
(18) ||у* — у|| < 1 ^^ • ||\аА\ • \х\ +
(19) ||у* — у|| < ^ + 2)и ||\аА\ • \х\ +
||у*\\ -1 - (N + 2)u \\aAx +
При M = N полученные оценки справедливы также и для операции с транспонированной матрицей, у ^ aATx + ßy.
6. Оценка производительности
В данном разделе представлены выполненные эксперименты и полученные результаты. Все эксперименты выполнялись на видеокарте NVIDIA GeForce GTX 1080 GPU (8 GB GDDR5X, 2560 CUDA cores, Compute Capability 6.1). Хост имел следующую конфигурацию: Intel Core i5 7500 (3.40 GHz) / 16 GB DDR4 RAM / Ubuntu 18.04.4 LTS / CUDA 10.2 / NVIDIA Driver 440.33.01. Исходный код компилировался с использованием компилятора nvcc с опциями -O3 -use_fast_math -Xcompiler=-O3,-fopenmp,-ffast-math.
6.1. Эффективность отдельных CUDA ядер
В данном эксперименте с использованием средства профилирования NVIDIA Visual Profiler мы исследовали эффективность использования ресурсов GPU отдельными вычислительными ядрами (глобальными функциями), запускаемыми в рассматриваемой реализации GEMV. Эксперимент выполнялся с нетранспонированной матрицей размера M = N = LDA = 1000 с 424-битной точностью вычислений. Для достижения указанной точности был использован набор из 32 модулей СОК, каждый 32 бита, с динамическим диапазоном M = 850 бит. Точность вычислений с числами вида (5) составляет [log^ %/Mj — 1 бит [11]. Необходимо отметить, что хотя множественные округления могут оказать существенное влияние на общую производительность нашей реализации, мы всегда можем снизить их количество, используя набор модулей СОК большего размера. В связи с этим исходные данные для эксперимента были такими, что множественные округления не требовались.
В таблице 3 представлен вклад отдельных CUDA ядер в общее время выполнения операции. При указанных параметрах эксперимента наиболее трудоемким является суммирование элементов матрицы B по строкам, на которое приходится 71.5% от общего времени вычислений. Также существенный вклад вносит вычисление цифр многоразрядных мантисс при формировании B.
Таблица 3. Вклад отдельных CUDA ядер в общее время выполнения GEMV многократной точности
Функция Описание %
Умножение x на a, y на в: вычисление
mp vec2scal mul esi знаков, экспонент, интервальных оценок 0.2
mp vec2scal mul v Умножение x на a, y на в: вычисление digits , '" г - ь цифр многоразрядных мантисс 0.1
mp mat2vec right Формирование матрицы B: вычисле-scal esi ние знаков, экспонент, интервальных оценок 5.6
mp mat2vec right scal di ^Формирование матрицы B: вычисле-— g ние цифр многоразрядных мантисс 20.8
mp vector round Округление промежуточных результатов 1.7
mp_ . ■ Построчная редукция матрицы B и matrix row sum ^ a — — сложение с вектором py 71.5
В таблице 4 представлены значения фактической загрузки доступных ресурсов GPU (occupancy), эффективной пропускной способности глобальной памяти (bandwidth) и ветвления логики исполнения (divergence) для двух наиболее трудоемких CUDA ядер. Разбиение операций многократной точности на части (mp—mat2vec—right—scal—digits) приводит к существенно более эффективным реализациям по сравнению с традиционным подходом, где каждая операция многократной точности выполняется одним потоком GPU (mp—matriX—roW—sum).
Таблица 4. Показатели эффективности наиболее затратных вычислительных ядер
Ядро Occupancy, % Divergence, % Bandwidth, ГБ/с
mp mat2vec right scal digits 95.5 0 163.3
mp matrix row sum 26.8 66.8 49.9
Отметим, что NVIDIA 1080 GPU использует память GDDR5X с рабочей частотой 1251 МГц и шириной шины 256 бит, следовательно теоретическая пиковая пропускная способность памяти составляет
(20) 1251 х 106 х (256/8) х 8/109 = 320.3 ГБ/с.
Для ядра mp_mat2vec_right_scal_digits эффективная пропускная способность составляет 51% от пиковой в связи с тем, что вычисление каждой цифры мантиссы сопровождается нахождением остатка от деления по соответствующему модулю (mod), но, как известно, операция mod является одной из наиболее трудоемких арифметических операций. Поэтому дальнейшее повышение быстродействия вычислений с цифрами многоразрядных мантисс в СОК возможно за счет ускорения операции mod. Этого можно достичь, например, используя редукцию Барретта. Для ядра mp_matrix_row_sum фактическая пропускная способность составляет 16% от пиковой, так как все промежуточные результаты хранятся в разделяемой памяти, а в глобальную память записывается лишь финальный результат работы каждого блока потоков. Однако разделяемая память в данном случае является ограничивающим фактором для полного использования доступных ресурсов GPU.
6.2. Сравнение с другими реализациями
В ходе экспериментов быстродействие представленной реализации GEMV сравнивалось с реализациями, основанными на CUDA библиотеках многократной точности GARPREC, CAMPARY и CUMP, которые рассматриваются в разделе 1. Все эти библиотеки основаны на позиционной системе счисления. Для того, чтобы оценить практический эффект от использования структуры mp_array_t и разбиения операций многократной точности на несколько CUDA ядер, мы реализовали классический алгоритм GEMV, в котором каждая арифметическая операция с многоразрядными числами вида (5) выполняется одним потоком GPU, а для хранения массива многоразрядных чисел используется массив структур mp_float_t. В дальнейшем эта версия GEMV называется "базовая реализация".
Эксперименты выполнялись при фиксированном размере матрицы M = N = LDA = 1000. Рассматривались различные уровни арифметической точности, от 106 до 1696 бит. Исходные наборы данных были представлены псевдослучайными многоразрядными числами с плавающей точкой, равномерно распределенными в интервале [— 1,1]. Чтобы уменьшить влияние шума, оцениваемые реализации GEMV повторялись в цикле несколько раз, измерялось общее время выполнения всех итераций в миллисекундах, после чего рассчитывалось среднее время выполнения одной итерации.
Таблица 5. Время выполнения (в миллисекундах) GEMV многократной точности на GeForce GTX 1080.
Точность в битах
106 212 424 848 1696
Предложенная реализация 3.1 5.6 8.5 12.9 24.6
Предложенная реализация (трансп.) 2.9 4.6 7.1 11.2 19.4
Базовая реализация 12.1 22.8 42.4 75.5 150.6
CUMP 20.1 20.8 26.8 56.8 154.0
CAMPARY 0.7 6.0 26.3 124.4 2499.9
GARPREC 26.5 37.4 70.2 206.0 689.1
Результаты экспериментов, представленные в таблице 5, показывают, что во многих случаях быстродействие предложенной реализации существенно выше, чем у аналогов, а также всегда выше чем у базовой
реализации. С ростом точности время предложенной реализации возрастает линейно. Операция с транспонированной матрицей оказалась незначительно быстрее операции с нетранспонированной матрицей, что обусловлено лучшим паттерном чтения матрицы B при ее редукции по столбцам.
При точности вычислений 106 бит (четверная точность) библиотека CAMPARY оказывается быстрее нашей реализации, однако при увеличении точности время CAMPARY существенно возрастает, что соответствует приведенным в разделе 1 оценкам сложности арифметических операций для данной библиотеки.
7. Выводы
В этой статье мы представили параллельную реализацию операции GEMV многократной точности для вычислительных систем с CUDA-совместимыми графическими процессорами. В основе нашей реализации лежит двухуровневая схема распараллеливания вычислений: на верхнем уровне осуществляется параллельная обработка элементов вектора/матрицы, а на нижнем уровне параллельно обрабатываются отдельные цифры многоразрядных мантисс элементов. Последнее стало возможным за счет перехода от традиционной двоичной арифметики к системе остаточных классов (СОК).
Однако простое использование СОК на GPU не дает значительного эффекта, так как операции с плавающей точкой многократной точности содержат последовательные и параллельные секции, что приводит к ветвлению путей выполнения. Для преодоления этого недостатка в предложенной реализации GEMV использована оригинальная схема вычислений, основанная на разбиении арифметических операций многократной точности на составные части, каждая из которых выполняется отдельным CUDA ядром. Такая схема увеличивает количество обращений к глобальной памяти устройства, так как все промежуточные результаты вычислений должны храниться в глобальной памяти. Однако это компенсируется тем, что ветвления внутри вычислительных ядер исключаются и обеспечивается более полная загрузка GPU. Эффективность такой схемы подтверждается экспериментами.
Разработанная реализация GEMV входит в состав новой GPU-библиотеки многократной точности MPRES-BLAS, которая доступна
на GitHub. Отметим, что функциональность MPRES-BLAS не ограничена базовыми операциями линейной алгебры. Фактически MPRES-BLAS является также арифметической библиотекой общего назначения, поскольку предоставляет CPU и CUDA реализации основных арифметических операций многократной точности. Кроме этого, в MPRES-BLAS реализованы оптимизированные алгоритмы для вычислений в СОК, такие как сравнение по величине и масштабирование степенью двойки. Также библиотека поддерживает арифметику расширенного диапазона (extended-range) для CPU и GPU, которая позволяет исключить переполнение и потерю значимости в процессе вычислений с разномасштабными величинами.
Направлениями дальнейших исследований является автоматический подбор оптимальной конфигурации вычислительных ядер в зависимости от точности вычислений и размера задачи. Также мы планируем адаптировать использованные подходы для реализации разреженного матрично-векторного умножения (SpMV) многократной точности на GPU.
Список литературы
[1] M. Courbariaux, Y. Bengio, J. David. Training deep neural networks with low precision multiplications, 2014. arXivJyJ 1412.7024 33
[2] D. H. Bailey, J. M. Borwein. "High-precision arithmetic in mathematical physics", Mathematics, 3:2 (2015), pp. 337-367. I t34
[3] J. Danek, J. PospiSil. "Numerical aspects of integration in semi-closed option pricing formulas for stochastic volatility jump diffusion models", International Journal of Computer Mathematics, 97:6 (2020), pp. 1268-1292.
d ' 34
[4] Y. Feng, J. Chen, W. Wu. "The PSLQ algorithm for empirical data", Math. Сотр., 88:317 (2019), pp. 1479-1501. 34
[5] S. Leweke, E. von Lieres. "Fast arbitrary order moments and arbitrary precision solution of the general rate model of column liquid chromatography with linear isotherm", Comput. Chem. Eng., 84 (2016), pp. 350-362. 34
[6] M. Kyung, E. Sacks, V. Milenkovic. "Robust polyhedral Minkowski sums with GPU implementation", Comput. Aided Des., 67—68 (2015), pp. 48-57.
d ' 34
[7] B. Pan, Y. Wang, S. Tian. "A high-precision single shooting method for solving hypersensitive optimal control problems", Mathematical Problems in Engineering, 2018 (2018), 7908378, 11 pp. 34
[8] Y. Xuan, D. Li, W. Han. "Efficient optimization approach for fast GPU computation of Zernike moments", Journal of Parallel and Distributed Computing, 111 (2018), pp. 104-114. 34
[9] C. L. Lawson, R. J. Hanson, D. R. Kincaid, F. T. Krogh. "Basic linear algebra subprograms for Fortran usage", ACM Trans. Math. Softw., 5:3 (1979), pp. 308—323. t34
[10] R. Nath, S. Tomov, T. Tim Dong, J. Dongarra. "Optimizing symmetric dense matrix-vector multiplication on GPUs", SC '11: Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis, ACM, New York, NY, USA, 2011, pp. 1-10. 36 36
[11] K. Isupov, V. Knyazkov, A. Kuvaev. "Design and implementation of multiple-precision BLAS Level 1 functions for graphics processing units", Journal of Parallel and Distributed Computing, 140 (2020), pp. 25-36.
^35,36,38,40,43,44,46,47,49,52
[12] A. Omondi, B. Premkumar. Residue number systems: theory and implementation, Imperial College Press, London, UK, 2007. 35 39
[13] K. Bigou, A. Tisserand. "Single base modular multiplication for efficient hardware RNS implementations of ECC", Cryptographic hardware and embedded systems - CHES 2015, eds. T. Giineysu, H. Handschuh, Springer Berlin Heidelberg, Berlin, Heidelberg, 2015, pp. 123-140. 35
[14] A. Abdelfattah, D. Keyes, H. Ltaief. "KBLAS: an optimized library for dense matrix-vector multiplication on GPU accelerators", ACM Trans. Math. Softiu., 42:3 (2016), 18. 36 40
[15] G. He, J. Gao, J. Wang. "Efficient dense matrix-vector multiplication on GPU", Concurrency and Computation: Practice and Experience, 30:19 (2018), e4705. 36
[16] T. Inoue, H. Tokura, K. Nakano, Y. Ito. "Efficient triangular matrix vector multiplication on the GPU", PPAM 2019: Parallel Processing and Applied Mathematics, Lecture Notes in Computer Science, vol. 12043, eds. R. Wyrzykowski, E. Deelman, J. Dongarra, K. Karczewski, Springer International Publishing, Cham, 2020, pp. 493-504. 136
[17] Quadruple precision BLAS routines for GPU: QPBLAS-GPU ver.1.0. User's manual, 2013 (accessed 19 May 2019), 58 pp. .urn. 36 38
[18] R. Iakymchuk, S. Collange, D. Defour, S. Graillat. "ExBLAS: reproducible and accurate BLAS library", Numerical Reproducibility at Exascale (NRE2015) workshop held as part of the Supercomputing Conference (SC15) (November 20, 2015, Austin, TX, USA), url 36 38
[19] D. Mukunoki, T. Ogita. "Performance and energy consumption of accurate and mixed-precision linear algebra kernels on GPUs", Journal of Computational and Applied Mathematics, 372 (2020), 112701. I t36 37 38
[20] Y. Hida, X.S. Li, D.H. Bailey. "Algorithms for quad-double precision floating point arithmetic", Proceedings 15th IEEE Symposium on Computer Arithmetic, ARITH-15 (11-13 June 2001, Vail, CO, USA), pp. 155-162.
^36
[21] D.E. Knuth. The art of computer programming. V. 2: Seminumerical algorithms, 3rd ed., Addison-Wesley Longman Publishing Co., Inc., USA, 1997, ISBN 978-0201896848. 36
[22] J.R. Shewchuk. "Adaptive precision floating-point arithmetic and fast robust geometric predicates", Discrete & Computational Geometry, 18:3 (1997), pp. 305-363. 36
[23] T. Ogita, S. M. Rump, S. Oishi. "Accurate sum and dot product", SI AM J. Sei. Comput., 26:6 (2005), pp. 1955—1988. 37
[24] M. Lu, B. He, Q. Luo. "Supporting extended precision on graphics processors", DaMoN'10: Proceedings of the Sixth International Workshop on Data Management on New Hardware (2010, Indianapolis, Indiana, USA), pp. 19-26. d t
37 38
[25] M. Joldes, J. Muller, V. Popescu. "Implementation and performance evaluation of an extended precision floating-point arithmetic library for high-accuracy semidefinite programming", 2017 IEEE 24th Symposium on Computer Arithmetic (ARITH) (24-26 July 2017, London, UK), pp. 27-34.
d t
37 38
[26] T. Nakayama, D. Takahashi. "Implementation of multiple-precision floatingpoint arithmetic library for GPU computing", The 23rd IASTED International Conference on Parallel and Distributed Computing and Systems PDCS 2011 (December 14-16 2011, Dallas, USA), pp. 343-349. 38
[27] K. Isupov. "Using floating-point intervals for non-modular computations in residue number system", IEEE Access, 8 (2020), pp. 58603-58619.
^39,40,47
[28] N. J. Higham. Accuracy and stability of numerical algorithms, 2nd, SIAM, Philadelphia, PA, USA, 2002, ISBN 978-0-89871-521-7, xxvii+663 pp.
^49,50,51
[29] J. Muller, N. Brunie, F. de Dinechin, C. Jeannerod, M. Joldes, V. Lefevre, G. Melquiond, N. Revol, S. Torres. Handbook of floating-Point arithmetic, 2, Birkhäuser, Basel, 2018, ISBN 978-3-319-76525-9. 60
Поступила в редакцию 29.04.2020 Переработана 24.07.2020
Опубликована 20.08.2020
Рекомендовал к публикации
к.ф.-м.н. С. А. Романенко
Пример ссылки на эту публикацию:
К. С. Исупов, В. С. Князьков. «Матрично-векторное умножение многократной точности на графическом процессоре». Программные системы: теория и приложения, 2020, 11:3(46), с. 33-59.
10.25209/2079-3316-2020-11-3-33-59 url, http: //psta.psiras. ru/read/psta2020_3_33-59 .pdf
Эта же статья по-английски: 10.25209/2079-3316-2020-11-3-61-84
Об авторах:
Константин Сергеевич Исупов
Кандидат технических наук, доцент кафедры электронных вычислительных машин Вятского государственного университета. Область научных интересов: высокоточные вычисления, система остаточных классов, компьютерная арифметика, параллельные алгоритмы, программирование для графических процессоров, GPGPU.
ПИ 0000-0003-0239-0404 e-mail: [email protected]
Владимир Сергеевич Князьков
Доктор технических наук, профессор, главный научный сотрудник научно-исследовательского института фундаментальных и прикладных исследований, Пензенский государственный университет. Область научных интересов: параллельные вычислительные архитектуры, высокопроизводительные вычисления, реконфигурируемые вычислительные системы.
ПИ 0000-0003-3820-6541 e-mail: [email protected]
Sample citation of this publication:
Konstantin S. Isupov, Vladimir S. Knyazkov. "Multiple-precision matrix-vector multiplication on graphics processing units". Program Systems: Theory and Applications, 2020, 11:3(46), pp. 33-59. (In Russian). 10.25209/2079-3316-2020-11-3-33-59
url http://pi3ta.pi3irai3.ru/read/pi3ta2020_3_33-59.pdf
The same article in English:
10.25209/2079-3316-2020-11-3-61-84
d
