CC BY
0
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ / _MATHEMATICAL MODELING, NUMERICAL METHODS AND PROGRAM COMPLEXES_
DOI: https://doi.org/10.60797/IRJ.2024.143.115
РЕАЛИЗАЦИЯ МЕТОДА СОПРЯЖЕННЫХ ГРАДИЕНТОВ НА ГРАФИЧЕСКОМ ПРОЦЕССОРЕ С ПРИМЕНЕНИЕМ МАТРИЧНЫХ МЕТОДОВ РЕШЕНИЯ СЛАУ
Научная статья
Устюжанин Д.А.1, Пицхелаури С.С.2' *, Некрасов К.А.3
2 ORCID : 0000-0003-2192-1880;
3 ORCID : 0000-0002-1863-2597;
1 2' 3 Уральский федеральный университет имени первого Президента России Б.Н. Ельцина, Екатеринбург, Российская
Федерация
* Корреспондирующий автор (lauri2011[at]mail.ru)
Аннотация
В работе представлен алгоритм решения систем линейных алгебраических уравнений (СЛАУ) методом сопряженных градиентов на графическом процессоре. В качестве основного инструмента реализации предложен класс-контейнер для работы с матрицами на базе технологии NVIDIA CUDA. Производительность решения СЛАУ методами Крамера и сопряженных градиентов была сопоставлена на центральном и графическом процессорах. Результаты исследования показали, что распараллеленный метод сопряженных градиентов, выполненный на графическом процессоре, обладает наибольшей эффективностью при обработке СЛАУ с симметричной положительно определенной основной матрицей. Рассматриваемый подход к параллельной обработке матричных операций имеет потенциал для применения в различных областях, где требуется решение крупных систем уравнений, таких как в науке, инженерии и финансах. В целом, данная работа представляет практическую значимость в области оптимизации производительности вычислений и закладывает фундамент для решения многих математических задач на графическом процессоре.
Ключевые слова: распараллеливание алгоритмов, решение СЛАУ на графическом процессоре, матричные операции.
REALIZATION OF THE METHOD OF CONJUGATE GRADIENTS ON A GRAPHIC PROCESSOR WITH APPLICATION OF MATRIX METHODS OF SLAU SOLUTION
Research article
Ustyuzhanin D.A.1, Pitskhelaury S.S.2' *, Nekrasov K.A.3
2 ORCID : 0000-0003-2192-1880;
3 ORCID : 0000-0002-1863-2597;
1 2' 3 Ural Federal University, Ekaterinburg, Russian Federation
* Corresponding author (lauri2011[at]mail.ru)
Abstract
The paper presents an algorithm for solving SLAEs by the method of conjugate gradients on a graphics processor. A matrix container class based on NVIDIA CUDA technology is proposed as the main implementation tool. The performance of SLAE solution by Cramer and conjugate gradient methods was compared on the central and graphic processors. The results of the study have shown that the parallelized method of conjugate gradients performed on a graphics processor has the highest efficiency when processing SLAEs with a symmetric positively defined principal matrix. The considered approach for parallel processing of matrix operations has potential for application in various areas where large systems of equations need to be solved, such as in science, engineering and finance. Overall, this work is of practical relevance in the field of computational performance optimization and lays the foundation for solving many mathematical problems on a GPU.
Keywords: parallelization of algorithms, solving SLAEs on a graphics processor, matrix operations.
Введение
Системы линейных алгебраических уравнений (СЛАУ) являются фундаментом задач для многих областей естествознания, включая механику, электротехнику, экономику и компьютерную графику. Способы решения СЛАУ имеют важное значение, особенно при обработке большого количества уравнений. Методы решения СЛАУ делятся на точные и приближенные, включая итерационные методы. Алгоритмы для вычислительных машин, основанные на точных методах, таких как метод Гаусса и метод Крамера, подходят для систем порядка не выше 20, в то время как итерационные методы могут обрабатывать системы более высокого порядка.
Одним из подходов к решению СЛАУ является параллельное программирование, включая использование графических процессоров (GPU). В отличие от центрального процессора (CPU), GPU имеет сотни или тысячи ядер, способных обрабатывать вычисления одновременно. Это делает GPU более эффективным для параллельных задач, таких как обработка графики, научные вычисления и машинное обучение. Кроме того, на GPU реализуются различные методы решения СЛАУ.
Метод сопряженных градиентов относится к итерационным методам и основывается на умножении матриц и вычислении скалярных произведений векторов. Умножение матриц строится на основе нахождения суммы
произведении соответствующих элементов исходных матриц, т. е. состоит из множества несвязных вычислении. Поэтому метод сопряженных градиентов подходит для распараллеливания на графическом процессоре.
Работа с матрицами при использовании технологии CUDA
Архитектура CUDA представляет собоИ специальный набор инструментов и библиотек, предназначенный для программирования графических процессоров. В качестве языков программирования CUDA поддерживает Python, C, C++ и т. д. Распараллеливание вычислений на GPU достигается за счет запуска большого количества потоков. Потоки в архитектуре CUDA представляются в виде сложной структуры. На верхнем уровне располагается сетка, которая является дискретным трехмерным пространством блоков, каждый блок в свою очередь трехмерным пространством потоков.
Матрицы заданного размера, состоящие из действительных чисел, образуют следующее линейное пространство:
^ Пе\г\ ... Пг\___ 1 ^
Mr
«00
^ ап-10
ЯОт-1
Яи-lm-l /
Ч
(1)
В настоящей работе матрицы представляли как одномерные динамические массивы. Доступ к элементу матрицы осуществляли по следующему соотношению: + { * М + ] / = 0 N — 1 = 0 М - 1) ' гДе N - число строк в
матрице, М - число столбцов в матрице, р1г - указатель на первый элемент массива. В качестве «обертки» для динамического массива использовали класс-контейнер devMatгix.
Таблица 1 - Основные методы класса devMatrix DOI: https://doi.Org/10.60797/IRJ.2024.143.115.1
Сигнатура метода Назначение метода
devMatrix () ; Конструктор по умолчанию.
devMatrix(size t N, size t M); Конструктор, выделяет память, но не инициализирует данные.
devMatrix(size t N, size t M, double a); Конструктор, инициализирует данные переданным значением.
devMatrix(devMatrix&& A); Конструктор, захватывает данные rvalue объекта.
-"devMatrix() ; Деструктор, освобождает память GPU.
devMatrix operator+(devMatrix const& A); Операция сложения матриц.
devMatrix operator*(devMatrix const& A); Операция умножения матриц.
Класс devMatrix включает все стандартные специальные методы (см. таб. 1), в том числе перемещающий конструктор и оператор присваивания. Класс Matrix является прообразом класса devMatrix и предназначен для работы с матрицами на CPU. Для удобства devMatrix имеет конструктор от Matrix.
Сложение матриц на GPU
В качестве иллюстрации внутреннего устройства класса devMatrix и работы CUDA, рассмотрим перегрузку оператора «+». Данный оператор осуществляет сложение двух матриц. Под сложением матриц будем понимать следующую операцию:
+ ВПт —
яоо
Я0т-1
* / ¿00 +
\ Яц-10 ■■ ■ ( аоо + boo
(in- ю + bn- ю
Яц-lm-l /
bn-10 ЯОт-1 + bom-
ûn-lm-l + Ьп-im-i )
bom-l \
bn-lm-1 /
(2)
При использовании CUDA, сумма каждой пары элементов (ац, Ьц) вычисляется в отдельном потоке. В данном случае использовали двумерное представление сетки и блока, при этом число блоков в одном из измерений вычисляли как ближайшее целое, большее или равное числу Ц . Блоки имели размер 32 х 32 потока. Каждый поток
обрабатывал одну сумму (см. рис. 1).
Рисунок 1 - Блок-схема суммы матриц DOI: https://doi.org/10.60797/IRJ.2024.143.115.2
Рисунок 2 - Сравнение производительности сложения матриц на CPU и GPU DOI: https://doi.org/10.60797/IRJ.2024.143.115.3
Примечание: оранжевая и синяя линии отображают зависимость времени от размерности матрицы для операции сложения в логарифмических координатах на CPU и GPU соответственно
Распараллеливание на GPU, при больших размерностях матриц, значительно сокращает время обработки операции сложения (см. рис. 2). При этом стоит отметить, что возрастание времени не связано на прямую с процессом суммирования элементов на GPU. Возрастание времени в данном случае определяется тем, что процесс выделения непрерывного участка памяти требует дополнительных вычислений. Также в данном случае размерность матрицы сильно ограничена объемом глобальной памяти, например квадратная матрица порядка 25000, состоящая из объектов типа double, занимает:
250002-8 10243
4.66 ГБ
(3)
Для корректной работы, стоит учитывать, что в какой-то момент времени будет как минимум два объекта типа devMatrix, при условии, что поддерживается перемещающий конструктор. Идеология перемещающих конструктора и оператора играют важную роль в ОПП реализации матричных операций. При манипуляции большими объемами данных стоит избегать их дублирования, которое возникает при инициализации или переопределении в результате возвращения из функции rvalue объектов. В противном случае минимальное количество объектов при обработке операции возрастает до трех, что соответствует (в случае квадратной матрицы порядка 25000) 13.98 ГБ. Если отбросить ограничения, связанные с памятью, на современных видеокартах на размерность сетки накладывается ограничение в максимальное число блоков, что соответствует максимальной размерности первого измерения. Максимальная размерность блока составляет 1024 потока (подробные данные см. гайд ссылка на гайд).
Таким образом максимальное количество строк/столбцов в квадратной матрице определяется следующим образом: Т к 655362 > 231 — 1 ' N ~ 32 • — 1 с учетом что в нашей реализации К должно быть кратно 32; ДГ = 1482880, что во много раз превышает лимит памяти.
В результате верхнюю границу размеров матрицы определяет именно объем глобальной память GPU.
Умножение матриц на GPU
Умножение матриц в алгебраической форме записывается следующим образом:
Anxm " BmXk — С„хЬ cij — aid ' ^dj
(4)
Классический распараллеленный алгоритм произведения матрицы схож с суммированием, но для каждого элемента отдельно запускается цикл, в котором рассчитывается сумма произведений (см. рис. 3). Важно, чтобы промежуточные значения суммы хранились в быстрой регистровой памяти GPU, иначе произойдет падение производительности программы.
Рисунок 3 - Блок-схема произведения матриц DOI: https://doi.org/10.60797/IRJ.2024.143.115.4
Рисунок 4 - Сравнение производительности умножения матриц на CPU и GPU DOI: https://doi.Org/10.60797/IRJ.2024.143.115.5
Примечание: оранжевая и синяя линии отображают зависимость времени от размерности матрицы для операции умножения в логарифмических координатах на CPU и GPU соответственно
На данную реализацию произведения матриц накладываются те же ограничения, что на суммирование. Разница во времени между перемножением на CPU и GPU, значительно превосходит разницу в случае суммы (см. рис. 4).
Представленный алгоритм является самой простой реализацией параллельного произведения. Существуют альтернативные алгоритмы:
1. Распараллеливание не только вычисления отдельного элемента результирующей матрицы, но и отдельных произведений соответствующей суммы, за счет третьего измерения z в сетке. У данного алгоритма есть две глобальные проблемы. Первая заключается в гонке потоков на запись произведений, появляется она из-за того, что мы вынуждены обращаться не к регистровой памяти, а к глобальной. Причем в подобной реализации в сочетание с нашим представлением класса devMatrix за суммирование произведений одного элемента не всегда отвечает одни блок, поэтому накладывается ограничение на использование разделяемой памяти. Вторая проблема связана с резким уменьшением максимальной размерности матриц из-за уменьшения числа блоков в x, y измерении из-за добавления z измерения: _
N » 32 • ^231 - 1 т. к. Ц должно быть кратно 32: N = 40960.
2. Рассматривать сетку в двух измерениях соответствуя элементам результирующей матрицы. В каждом блоке будет происходить распараллеленное вычисление суммы произведений с использованием разделяемой памяти. Данный алгоритм позволит вычислять произведения матриц с самой высокой скоростью, но если будет происходить перемножение матриц _АпХт • ВтХ^ , то максимальное значение к (число столбцов матрицы A, число строк матрицы B) равно 1024, из-за ограничения на число потоков внутри блока.
По сравнению с указанными альтернативами, классический (простой) алгоритм предпочтителен, он обеспечивает хорошую скорость около граничной области размерности матриц по памяти.
Метод сопряженных градиентов
Пусть система из п линейных алгебраических уравнений фоРме: А-Х = В,
' «11 Я12 • • • а\п 4
«21 Й22 • ■ • а2п
c n неизвестными представлена в матричной
где А =
\ ап\ аП2
основная матрица системы,
X -
I XX *2
\ лп /
*пп /
- столбец неизвестных, jj —
( h \ Ъг
Ьп
- столбец свободных членов.
Метод сопряженных градиентов, подразумевает, что основная матрица системы А - симметрична и положительно определена, т. е. \/хфв(=Яп:(А-х,х)>0 .
Для решения такой системы методом сопряженных градиентов можно использовать следующее рекуррентное соотношение: Шаг 0:
= А-Хо -ДА) = ДХ! = ft • Г0) ОД = 0
Х0 =
' \
Х2
ГО
- произвольный вектор.
(5)
\ хп / _
Шаг индукции /(i = l,B-l)
X,- = X,- + ДХг, где ДХ/ = at-\ ■ ДХг_1 + $-1 • r/_i(i Ф 1)
п = А ■ Xi - В,
_ (п,г,) (А-п,АХ1)-(г,,АХ,) (А-п,г,) 1 (ЛДХьДХг) ■(A-ri,rl)-{A-ri,KXi)2
а _ (П.АХ,-) • (Ат,-,АХ,-)-(п,п) • (Л-АХ,-,АХ,-) Pl (A-AXi,AXl)-(A-rl,ri)-(.A-rl,AXl)2 '
(6)
(7)
(8) (9)
Отметим, что метод сопряженных градиентов всегда сходится к решению исходной СЛАУ за £ < п итераций ( п - порядок системы).
return X;
Рисунок 5 - Вид итогового алгоритма с учетом построенных инструментов для работы с матрицами
DOI: https://doi.Org/10.60797/IRJ.2024.143.115.6
Расчет коэффициентов происходит в специальной функции, которая распараллеливает скалярное
произведение. Ядро реализовали с использованием редукции (см. рис. 6). Разбиение скалярных произведений происходит за счет добавления измерения у с индексацией ^ _ о""^
Рисунок 6 - Блок-схема реализации скалярного произведения DOI: https://doi.Org/10.60797/IRJ.2024.143.115.7
Рисунок 7 - Метод Крамера, спряжённых градиентов на CPU и GPU DOI: https://doi.Org/10.60797/IRJ.2024.143.115.8
Примечание: зеленая, оранжевая и синяя линии отображают зависимость времени от размерности основной матрицы системы при решении СЛАУ в логарифмических координатах на CPU методом Крамера и МСГ и МСГ на GPU соответственно
О 200 400 600 800 1000 1200 1400 Размерность матрицы
Рисунок 8 - Метод спряжённых градиентов на GPU DOI: https://doi.Org/10.60797/IRJ.2024.143.115.9
Примечание: зависимость времени от размерности основной матрицы системы при решении СЛАУ методом МСГ на GPU
Итоговое сравнение производительности различных методов решения СЛАУ, показывает, что распараллеливание на GPU сильно уменьшает время (см. рис. 7, см рис. 8).
Все тесты проводились на следующей конфигурации:
Процессор: AMD Ryzen 5 5600 6-Core Processor 3.50 GHz;
Оперативная память: 16,0 ГБ;
Графический процессор: GeForce RTX 3060 12,0 ГБ, 3584 ядер CUDA.
Заключение
В настоящей работе предложены основные инструменты для обработки матричных операций с использованием технологии CUDA. На примере сложения и умножения матриц показано, что алгоритмы, построенные на распараллеливании вычислений, работают на несколько порядков эффективнее аналогичных последовательных программ на CPU. Определены границы размерностей матриц, которые можно обработать на GPU.
На основе библиотеки devMatrix, построен алгоритм решения СЛАУ методом сопряженных градиентов, который превзошел аналогичный алгоритм на CPU и метод Крамера по скорости решения. Полученный результат демонстрирует актуальность технологии CUDA в контексте научных расчетов, а также предоставляет возможность написания простых, но в то же время эффективных матричных алгоритмов на базе библиотеки devMatrix.
Финансирование
Работа выполнена на основе гранта Министерства образования РФ № FEUZ-2023-0013.
Конфликт интересов
Не указан.
Рецензия
Все статьи проходят рецензирование. Но рецензент или автор статьи предпочли не публиковать рецензию к этой статье в открытом доступе. Рецензия может быть предоставлена компетентным органам по запросу.
Funding
The work was carried out on the basis of a grant from the Ministry of Education of the Russian Federation No. FEUZ-2023-0013.
Conflict of Interest
None declared.
Review
All articles are peer-reviewed. But the reviewer or the author of the article chose not to publish a review of this article in the public domain. The review can be provided to the competent authorities upon request.
Список литературы / References
1. Крылов В. И. Вычислительные методы в 2 т. / В. И. Крылов, В. В. Бобков, П. И. Монастырный. — Москва : Наука, 1976.
2. Некрасов К. А. Параллельные вычисления общего назначения на графических процессорах / К. А. Некрасов, С. И. Поташников, А. С. Боярченков [и др.]. — Изд-во Уральского ун-та, 2016. — 103 с. — ISBN 978-5-7996-1722-6.
3. Тумаков Д. Н. Технология программирования CUDA / Д. Н. Тумаков, Д. Е. Чикрин, А. А. Егорчев [и др.]. — Казань : Казанский федеральный ун-т. — ISBN 978-5-00019-913-8.
4. Локтионов И. К. Численные методы / И. К. Локтионов. — Москва : Инфа-Инженерия, 2022. — 308 с. — ISBN 978-5-9729-0786-1.
5. Огородникова О. М. Вычислительные методы в компьютерном инжиниринге / О. М. Огородникова. — Изд-во Уральского ун-та, 2013. — 130 с. — ISBN 978-5-7996-0816-3.
6. Амосов А. А. Вычислительные методы для инженеров / А. А. Амосов, Ю. А. Дубинский, Н. В. Копченова. — Москва : Высш. шк, 1994. — 543 с. — ISBN 5-06-000625-5.
7. Пирумов У Г. Численные методы / У Г. Пирумов, В. Ю. Гидаспов, И. Э. Иванов [и др.]. — Москва : Юрайт, 2023. — 421 с. — ISBN 978-5-534-03141-6.
8. Снытиков А. В. Математическое моделирование и программная модель CUDA / А. В. Снытников, А. С. Колганов, Н. Н. Попова. — Москва : МАКС Пресс, 2018. — 171 с. — ISBN 978-5-317-05911-8.
9. Тоуманен Б. Программирование GPU при помощи Python и CUDA / Б. Тоуманен. — Москва : ДМК Пресс, 2020. — 235 с. — ISBN 978-5-97060-821-0.
10. Programming Guide // CUDA Toolkit Documentation. — URL: https://docs.nvidia.com/cuda/archive/11A0/cuda-c-programming-guide/index.html (accessed: 24.11.2023).
Список литературы на английском языке / References in English
1. Krylov V. I. Vychislitel'nye metody v 2 t. [Computational methods in 2 volumes] / V. I. Krylov, V. V. Bobkov, P. I. Monastyrny. — Moscow : Nauka, 1976. [in Russian]
2. Nekrasov K. A. Parallel'nye vychislenija obshhego naznachenija na graficheskih processorah [General-purpose parallel computing on graphics processors] / K. A. Nekrasov, S. I. Potashnikov, A. S. Boyarchenkov [et al.]. — Publishing House of the Ural University, 2016. — 103 p. — ISBN 978-5-7996-1722-6. [in Russian]
3. Tumakov D. N. Tehnologija programmirovanija CUDA [CUDA programming technology] / D. N. Tumakov, D. E. Chirkin, A. A. Egorchev [et al.]. — Kazan : Kazan Federal University. — ISBN 978-5-00019-913-8. [in Russian]
4. Loktionov I. K. Chislennye metody [Numerical methods] / I. K. Loktionov. — Moscow : Infa-Engineering, 2022. — 308 p. — ISBN 978-5-9729-0786-1. [in Russian]
5. Ogorodnikova O. M. Vychislitel'nye metody v komp'juternom inzhiniringe [Computational methods in computer engineering] / O. M. Ogorodnikova. — Publishing house of the Ural University, 2013. — 130 p. — ISBN 978-5-7996-0816-3. [in Russian]
6. Amosov A. A. Vychislitel'nye metody dlja inzhenerov [Computational methods for engineers] / A. A. Amosov, Yu. A. Dubinsky, N. V. Kopchenova. — Moscow : Higher School of Economics, 1994. — 543 p. — ISBN 5-06-000625-5. [in Russian]
7. Pirumov U. G. Chislennye metody [Numerical methods] / U. G. Pirumov, V. Y. Gidaspov, I. E. Ivanov [et al.]. — Moscow : Yurait, 2023. — 421 p. — ISBN 978-5-534-03141-6. [in Russian]
8. Snytnikov A. V. Matematicheskoe modelirovanie i programmnaja model' CUDA [Mathematical modeling and the CUDA software model] / A. V. Snytnikov, A. S. Kolganov, N. N. Popova. — Moscow : MAKS Press, 2018. — 171 p. — ISBN 978-5-317-05911-8. [in Russian]
9. Tumanen B. Programmirovanie GPU pri pomoshhi Python i CUDA [GPU programming using Python and CUDA] / B. Tumanen. — Moscow : DMK Press, 2020. — 235 p. — ISBN 978-5-97060-821-0. [in Russian]
10. Programming Guide // CUDA Toolkit Documentation. — URL: https://docs.nvidia.com/cuda/archive/11A0/cuda-c-programming-guide/index.html (accessed: 24.11.2023).
