Итерационный метод решения обобщенной задачи поиска собственных значений и векторов на основе LU-разложения правой матрицы с использованием технологии CUDA

Макошенко Д.М.

Итерационный метод решения обобщенной задачи поиска собственных значений и векторов на основе LU-разложения правой матрицы с использованием технологии CUDA

02.03.02 Фундаментальная информатика и информационные технологии

Информатика

Дипломы

Вуз: Белгородский государственный университет - национальный исследовательский университет (НИУ «БелГУ»)

ID: 5c1a62aa7966e104f6f8571e

UUID: 18149340-e5d0-0136-0166-525400005860

Язык: Русский

Опубликовано: почти 6 лет назад

Просмотры: 441

Макошенко Д.М.

Источник: Федеральное государственное автономное образовательное учреждение высшего образования «Белгородский государственный национальный исследовательский университет»

Комментировать 0

Рецензировать 0

Скачать - 882,1 КБ

Поделиться работой

ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ «БЕЛГОРОДСКИЙ ГОСУДАРСТВЕННЫЙ НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ» ( Н И У «Б е л Г У » ) ИНСТИТУТ ИНЖЕНЕРНЫХ ТЕХНОЛОГИЙ И ЕСТЕСТВЕННЫХ НАУК КАФЕДРА МАТЕМАТИЧЕСКОГО И ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ ИНФОРМАЦИОННЫХ СИСТЕМ Итерационный метод решения обобщенной задачи поиска собственных значений и векторов на основе LU-разложения правой матрицы с использованием технологии CUDA Выпускная квалификационная работа обучающегося по направлению подготовки 02.03.02 Фундаментальная информатика и информационные технологии очной формы обучения, группы 07001401 Макошенко Данила Михайловича Научный руководитель к.т.н., ст.пр. Лихошерстный А.Ю. БЕЛГОРОД 2018

СОДЕРЖАНИЕ ВВЕДЕНИЕ .......................................................................................................... 4 1. ОБЗОР СОВРЕМЕННЫХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧИ ПОИСКА СОБСТВЕННЫХ ЗНАЧЕНИЙ И ВЕКТОРОВ ................................................... 6 1.1 Собственные значения и вектора матриц ..................................................... 6 1.2 LU разложение матриц ................................................................................... 7 1.3 Существующие алгоритмы нахождения собственных значений и векторов ............................................................................................................................... 8 1.3.1 Решения обобщенной задачи поиска собственных значений на основе разложения Холецкого ......................................................................................... 8 1.3.2 Метод итераций для нахождения собственных значений и векторов .... 10 1.3.3 Метод А.М. Данилевского ........................................................................ 11 1.3.4 Метод А.Н. Крылова ................................................................................. 15 1.3.5 Метод Леверрье-Фаддеева ........................................................................ 18 1.4 Обзор технологии NVIDIA CUDA .............................................................. 22 1.4.1 Модель памяти CUDA............................................................................... 23 1.4.2 Математические библиотеки CUDA ........................................................ 24 1.5 Формат Compressed Sparse Row .................................................................. 26 2. РАЗРАБОТКА АЛГОРИТМА РЕШЕНИЯ ОБОБЩЕННОЙ ЗАДАЧИ ПОИСКА СОБСТВЕННЫХ ЗНАЧЕНИЙ И ВЕКТОРОВ ДЛЯ РАЗРЕЖЕННЫХ МАТРИЦ НА ОСНОВЕ LU-РАЗЛОЖЕНИЯ ПРАВОЙ ЧАСТИ УРАВНЕНИЯ ....................................................................................... 27 2.1 Преобразование обобщенной задачи поиска собственных значений и векторов к классической задаче ........................................................................ 27 2.2 Использование неполного LU-разложения матриц .................................... 28 2.3 Степенной метод со сдвигами нахождения собственных значений .......... 31

3. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ РАЗРАБОТАННОГО АЛГОРИТМА С ИСПОЛЬЗОВАНИЕМ ТЕХНОЛОГИИ CUDA ................................................ 33 3.1 Программная реализация параллельной версии алгоритма ....................... 33 3.2 Оценка эффективности разработанного алгоритма на основе проведения вычислительных экспериментов. ...................................................................... 38 ЗАКЛЮЧЕНИЕ .................................................................................................. 42 СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ ............................................ 44

ВВЕДЕНИЕ Задача поиска собственных значений и векторов часто предстает перед инженерами, которые проводят расчеты в таких промышленных сферах как строительство, машиностроение, авиастроение и т.п. Например, частоты колебаний выражаются через собственные значения. При обсчетах конструкций, большие части которых составляют пустоты, получаются разреженные матрицы. В связи с этим появляется проблема хранения и использования большого объёма данных в расчетах. Универсальных алгоритмов решения поставленной задачи не существует. Наиболее высокой скоростью обсчета являются вычисления на графическом процессоре, т.к. подобные решения обладают значительно большим количеством вычислительных ядер по сравнению с центральным процессором. Однако объем памяти графического ускорителя, как правило, составляет всего несколько гигабайт. Существующие методы решения собственных значений и векторов в процессе вычисления получают промежуточные плотные матрицы, что значительно увеличивает расход памяти. Предлагаемый метод сохраняет промежуточные результаты вычислений в разреженном виде, при этом не замедляя процесс вычислений. Целью выпускной квалификационной работы является разработка и программная реализация итерационного алгоритма решения обобщенной задачи поиска собственных значений и векторов для разреженных матриц на основе LU-разложения правой части уравнения с использованием технологии СUDA Поставленная цель достигается решением следующих задач:  Обзор существующих методов решения задачи поиска собственных значений и векторов;

 Разработка алгоритма решения обобщенной задачи поиска собственных значений и векторов для разреженных матриц на основе LU-разложения правой части уравнения;  Программная реализация разработанного алгоритма с использованием технологии CUDA;  Оценка эффективности разработанного алгоритма на основе проведения вычислительных экспериментов. В первой главе осуществлен обзор существующих методов решения задачи поиска собственных значений и векторов, дана базовая информация о собственных значениях и векторах, LU – разложении, формате CSR, а также технологии CUDA, которая использовалась при разработке параллельной версии алгоритма. Во второй главе были описаны все аспекты разработанного алгоритма, а именно преобразование обобщенной задачи поиска собственных значений и векторов к классической задаче, где использовало неполное LU – разложение, которой посвящена отдельная подглава. Так же был описан алгоритм решения классической задачи при помощи степенного метода, который позволяет решить её минимально расходуя память. В третьей главе описаны основные этапы разработки параллельной версии алгоритма с использованием технологии CUDA, а также приведены результаты испытаний программной реализации на матрицах различной размерности. Данная выпускная квалификационная работа состоит из 46 страниц, 3 рисунков, 4 таблиц и 35 использованных литературных источников. 5

1. ОБЗОР СОВРЕМЕННЫХ МЕТОДОВ РЕШЕНИЯ ЗАДАЧИ ПОИСКА СОБСТВЕННЫХ ЗНАЧЕНИЙ И ВЕКТОРОВ 1.1 Собственные значения и вектора матриц Собственный вектор – понятие в линейной алгебре, определяемое для квадратной матрицы как вектор, умножение матрицы на который даёт тот же вектор, умноженный на некоторое скалярное значение, называемое собственным значением матрицы. Таким образом, собственным вектором матрицы А называется такой ненулевой вектор x, что для некоторого λ выполняется равенство (1.1). Ax = λx (1.1) В свою очередь, собственным значением матрицы А называется такое число λ, для которого существует собственный вектор, т.е. равенство (1) имеет ненулевое решение x [1]. Обобщенной задачей поиска собственных значений и векторов называется задача вида (1.2). Ax = λBx, (1.2) где A и B - матрицы, λ – собственное значение, а x – собственный вектор. Для решения этой задачи, её необходимо свести к задаче поиска собственных значений и векторов матрицы общего вида (1.1). Задача нахождения собственных значений и собственных векторов матрицы возникает при решении многих прикладных задач физики, механики, астрономии, в которых требуется определить нетривиальное решение однородной системы линейных алгебраических уравнений вида (1) и тех значение числового параметра λ, при которых такое решение существует [2]. Оценка критических нагрузок при расчете строительных конструкция 6

основана на собственных значения и векторах матриц. Собственные числа и собственные векторы являются важнейшими характеристиками, отражающими существенные стороны линейных моделей. 1.2 LU разложение матриц LU-разложение – представление квадратной матрицы А в виде произведения матриц L и U, где L – нижняя треугольная матрица с единичной диагональю, а U – верхняя треугольная [5]. Такое представление удобно для решения системы алгебраических уравнений, т.к. оно позволяет перейти от решения исходной системы к последовательному решению систем с треугольными матрицами. Будем использовать следующие обозначения для элементов матриц L = (lij ), U = (uij ), i, j = 1, … , n, (1.3) где n – размерность матрицы. Алгоритм нахождения матриц L и U выполняется последовательно, т.к. на каждом следующем шаге для вычислений используются элементы, найденные ранее. Всего алгоритм состоит из m шагов. Для m = 1: u1j = a1j , j = 1,…,n lj1 = aj1 u11 , j = 2,…,n(u11≠0) (1.4) (1.5) Для m = 2, …, n: umj = a1j − ∑m k=1 lmk ukj , j = 1,…,n ljm = 1 umm (ajm − ∑m k=1 ljk ukm ), j = m+1,…,n (1.6) (1.7) 7

В программной реализации в памяти компьютера LU разложение не расходует лишнюю память, т.к. занимает ровно столько места, сколько занимает исходная матрица. Размерность L и U матриц такая же, как у исходной, а т.к. U – это верхняя диагональная матрица, а L –нижняя диагональная с единицами на главной диагонали, которые можно не хранить в памяти, то обе матрицы можно хранить как одну, где ниже главной диагонали будут элементы L матрицы, а выше и на диагонали U матрицы. Кроме того, на каждом шаге результаты вычислений можно сохранять на местах, освобождаемых последовательным затиранием элементов исходной матрицы. 1.3 Существующие алгоритмы нахождения собственных значений и векторов 1.3.1 Решения обобщенной задачи поиска собственных значений на основе разложения Холецкого Данный алгоритм применяется для решения обобщенной симметричной положительно определенной задачей собственных значений вида (1.8). 𝐴𝑥 = 𝜆𝐵𝑥, (1.8) где матрица A - симметричная, а матрица B - симметричная положительно определенная [3]. Очевидно, что эта задача легко сводится к задаче поиска собственных значений несимметричной матрицы общего вида, достаточно умножить обе части системы на B-1. Однако, в результате данного преобразования получится несимметричная задача поиска собственных значений, которая на порядок сложнее симметричной. 8

Поэтому для сведения к классической симметричной задаче воспользуется разложение Холецкого матрицы B (1.9). 𝐵 = 𝐿𝐿𝑇 , (1.9) где L – нижняя треугольная матрицы со строго положительными элементами на диагонали. Тогда получим равенство 1.10. 𝐴𝑥 = 𝜆𝐿𝐿𝑇 𝑥 (1.10) Т.к. при умножении матрицы LT на обратную ей матрицу L-T получится единичная матрицы, то предыдущая система эквивалента следующей: 𝐴𝐿−𝑇 𝐿𝑇 𝑥 = 𝜆𝐿𝐿𝑇 𝑥 (1.11) Умножим обе части уравнения на L-1: 𝐿−1𝐴𝐿−𝑇 𝐿𝑇 𝑥 = 𝜆𝐿𝑇 𝑥 (1.12) Для приведения задачи выполним следующие замены: С = 𝐿−1𝐴𝐿−𝑇 (1.13) 𝑦 = 𝐿𝑇 𝑥 (1.14) Тогда исходное уравнение сведется к следующему: С𝑦 = 𝜆𝑦, (1.15) Где C – симметричная матрица и т.о. получена хорошо известная и эффективно решаемая задача, где для собственные векторы исходной задачи могуть быть получены путем решения системы линейных уравнений с треугольной матрицей LT. 9

1.3.2 Метод итераций для нахождения собственных значений и векторов Для решения частичной проблемы собственных значений и собственных векторов в практических расчетах часто используется метод итераций (степенной метод) [6]. На его основе можно определить приближенно собственные значения матрицы A. Алгоритм метода итераций состоит из следующих шагов: 1. Выбрать произвольное начальное приближение собственного вектора X1(0) (индекс в скобках здесь и далее указывает номер приближения, а индекс без скобок соответствует номеру собственного значения или вектора); 2. Найти первые приближения собственных векторов (1.16) и значения (1.17) и положить k = 0. 𝑋1(1) = 𝐴𝑋1(0) (1) 𝜆1 (1.16) 1(1) = 𝑥𝑖 1(0) 𝑥𝑖 , (1.17) где i – любой номер 1 ≤ i ≤ n; 3. Вычислить (1.18) 𝑋1(𝑘+1) = 𝐴𝑋1(𝑘) (1.18) 4. Найти следующее приближение собственного значения: (𝑘+1) 𝜆1 1(𝑘+1) где 𝑥𝑖 1(𝑘) , 𝑥𝑖 1(𝑘+1) = 𝑥𝑖 1(𝑘) 𝑥𝑖 , (1.19) – соответствующие координаты векторов X1(k+1) и X1(k). При этом может быть использована любая координата с номером i, 1 ≤ i ≤ n; (𝑘+1) 5. Если ∆= |𝜆1 (𝑘) − 𝜆1 | ≤ ε, процесс завершить и положить 10

(𝑘+1) 𝜆1 ≅ 𝜆1 . (1.20) Если ∆> ε, положить k=k+1 и перейти к шагу 3. (𝑘+1) Вместо применяемой в пункте 4 алгоритма формы для 𝜆1 можно взять среднее арифметическое соответствующих отношений для разных координат. Используя 𝜆1, можно определить следующее значение 𝜆2 по формуле 1.21. 1(𝑘+1) 𝜆2 = 𝑥𝑖 1(𝑘) 𝑥𝑖 1(𝑘) −𝜆𝑥𝑖 (1.21) 1(𝑘−1) −𝜆𝑥𝑖 где i = 1, 2, …, n. Эта формула даёт грубые значения для 𝜆2 , так как значение 𝜆1 является приближенным. На основе (1.21) можно вычислять остальные 𝜆𝑗 (j = 3, 4, …, n). 1.3.3 Метод А.М. Данилевского Метод А.М. Данилевского заключается в преобразовании исходной матрицы 𝑎11 𝑎21 𝐴=( … 𝑎𝑛1 𝑎12 𝑎22 … 𝑎𝑛2 … … … … 𝑎1𝑛 𝑎2𝑛 … ) 𝑎𝑛𝑛 (1.22) в подобную ей матрицу Фронбениуса 𝑝1 1 𝑃=( … 0 𝑝2 0 … 0 … … … … 𝑝𝑛−1 0 … 1 𝑝𝑛 0 ) … 1 (1.23) 11

по формуле 𝑃 = 𝐵 −1𝐴𝐵 (1.24) с помощью матрицы подобия B. Переход от матрицы A к подобной ей матрице P осуществляется с помощью n-1 преобразований подобия, последовательно преобразующих строки матрицы A, начиная с последней, в соответствующие строки матрицы P. На первом этапе, предполагая что 𝑎𝑛,𝑛−1 ≠ 0, построим матрицу B1 (1.26), заменив в единичной матрицу порядка n элементы n-1 строки на значения 𝑏𝑛−1,𝑗 = − 𝑎𝑛𝑗 𝑎𝑛,𝑛−1 𝑗 ≠ 𝑛 − 1; (1.25) 𝑏𝑛−1,𝑛−1 = 𝐵1 = 1 0 … 𝑏𝑛−1,1 ( 0 0 1 … 𝑏𝑛−1,2 0 1 𝑎𝑛,𝑛−1 … … … … … 0 0 … 0 0 … 𝑏𝑛−1,𝑛−1 0 (1.26) 𝑏𝑛−1,𝑛 1 ) Умножим справа матрицу A на матрицу B1 𝐴 ∙ 𝐵1 = 𝐶 = 𝑐11 𝑐21 … 𝑐𝑛−1,1 ( 0 𝑐12 𝑐22 … 𝑐𝑛−1,2 0 … … … … … 𝑐1,𝑛−1 𝑐2,𝑛−1 … 𝑐𝑛−1,𝑛−1 1 𝑐1,𝑛 𝑐2,𝑛 … , (1.27) 𝑐𝑛−1,𝑛 0 ) где с𝑖𝑗 = 𝑎𝑖𝑗 + 𝑎𝑖,𝑛−1 ∙ 𝑏𝑛−1,𝑗 1 ≤ 𝑖 ≤ 𝑛, 𝑗 ≠ 𝑛 − 1; (1.28) с𝑖,𝑛−1 = 𝑎𝑖,𝑛−1 ∙ 𝑏𝑛−1,𝑛−1 1 ≤ 𝑖 ≤ 𝑛. 12

Тогда обратная матрица 𝐵1−1 имеет вид 𝐵1−1 1 0 = … 𝑎𝑛1 ( 0 0 1 … 𝑎𝑛2 0 … … … … … 0 0 … 𝑎𝑛,𝑛−1 0 0 0 … 𝑎𝑛𝑛 1 ) (1.29) Пусть 𝐷1 = 𝐵1−1 ∙ 𝐶. Т.к., очевидно, умножение слева матрицы С на матрицу 𝐵1−1 не изменяет последнюю строку C, то матрицы D1 имеет вид 𝐷1 = 𝑑11 𝑑11 … 𝑑𝑛−1,1 ( 0 𝑑12 𝑑22 … 𝑑𝑛−1,2 0 … … … … … 𝑑1,𝑛−1 𝑑2,𝑛−1 … 𝑑𝑛−1,𝑛−1 1 𝑑1,𝑛 𝑑2,𝑛 , … 𝑑𝑛−1,𝑛 0 ) (1.30) где 𝑑𝑖𝑗 = 𝑐𝑖𝑗 1≤ 𝑖 ≤𝑛−2 (1.31) 𝑛 𝑑𝑛−1,𝑗 = ∑ 𝑎𝑛𝑘 𝑐𝑘𝑗 1 ≤ 𝑗 ≤ 𝑛. 𝑘=1 Полученная матрица D1 подобна матрице А и имеет одну преобразованную строку. Этим заканчивается первый этап процесса. На втором этапе, предполагая, что dn-1,n-2 ≠ 0, построим матрицу B2, заменив в единичной матрице порядка n элементы n-2 строки на значения 13

𝑏𝑛−2,𝑗 = − 𝑑𝑛−1,𝑗 𝑑𝑛−1,𝑛−2 𝑗≠ 𝑛−2; (1.32) 𝑏𝑛−2,𝑛−2 = 𝐵2 = 1 0 … 𝑑𝑛−2,1 0 ( 0 0 1 … 𝑑𝑛−2,2 0 0 … … … … … … 1 𝑑𝑛−1,𝑛−2 0 0 … 𝑑𝑛−2,𝑛−1 1 0 0 0 … 𝑑𝑛−2,𝑛 0 1 ) (1.33) Далее, взяв в качестве матрицы A матрицу D1 и проведя вычисления по формулам (1.28) и (1.31), получим матрицу 𝐷2 = 𝐵2−1 ∙ 𝐷1 ∙ 𝐵2 с двумя преобразованными строками. Над матрицей D2 производим те же операции. Продолжая этот процесс, мы получим матрицу Фронбениуса −1 −1 𝑃 = 𝐵𝑛−1 ∙ 𝐵𝑛−2 ∙ … ∙ 𝐵2−1 ∙ 𝐵1−1 ∙ 𝐴 ∙ 𝐵1 ∙ 𝐵2 ∙ … ∙ 𝐵𝑛−2 ∙ 𝐵𝑛−1 , (1.34) Пусть 𝑦⃗ – собственный вектор матрицы P, отвечающий собственному значению 𝜆. Тогда 𝑃𝑦⃗ = 𝜆𝑦⃗ или в координатном виде 𝑝1 𝑦1 + 𝑝2 𝑦2 + ⋯ + 𝑝𝑛−1𝑦𝑛−1 + 𝑝𝑛 𝑦𝑛 = 𝜆𝑦1 𝑦1 = 𝜆𝑦2 𝑦2 = 𝜆𝑦3 … { 𝑦𝑛−1 = 𝜆𝑦𝑛 (1.35) Полагая 𝑦𝑛 = 1, мы получим, используя последовательно эти уравнения с низу вверх, 𝑦𝑛−1 = 𝜆, 𝑦𝑛−2 = 𝜆2 , …, 𝑦1 = 𝜆𝑛−1. 14

𝑦⃗ = (𝜆𝑛−1, 𝜆𝑛−2, … , 𝜆, 1) (1.36) Так как матрицы A и P подобны, то 𝐵 −1𝐴𝐵𝑦⃗ = 𝜆𝑦⃗ или 𝐴𝐵𝑦⃗ = 𝜆𝐵𝑦⃗. Это означает, что вектор 𝜆𝑛−1 𝜆𝑛−2 𝑥⃗ = 𝐵𝑦⃗ = 𝐵 ∙ … 𝜆 ( 1 ) (1.37) является собственным вектором матрица A, отвечающим собственному значению 𝜆. 1.3.4 Метод А.Н. Крылова Согласно тождеству Гамильтона-Кели, матрица A удовлетворяет своему характеристическому уравнению, поэтому 𝐴𝑛 − 𝑝1 𝐴𝑛−1 − 𝑝2 𝐴𝑛−2 − ⋯ − 𝑝𝑛−1𝐴 − 𝑝𝑛 𝐸 = 0 (1.38) 𝑦10 Возьмем произвольный ненулевой вектор 𝑦 ⃗⃗⃗⃗⃗0 = ( ⋮ ) и умножим обе 𝑦𝑛0 части равенства (1.38) справа на ⃗⃗⃗⃗⃗ 𝑦0 𝐴𝑛 ⃗⃗⃗⃗⃗ 𝑦0 − 𝑝1 𝐴𝑛−1⃗⃗⃗⃗⃗ 𝑦0 − 𝑝2 𝐴𝑛−2 ⃗⃗⃗⃗⃗ 𝑦0 − ⋯ − 𝑝𝑛−1𝐴𝑦 ⃗⃗⃗⃗⃗0 − 𝑝𝑛 𝐸𝑦 ⃗⃗⃗⃗⃗0 = 0 (1.39) Положим 15

𝑘 𝑦 ⃗⃗⃗⃗⃗ 𝑦0 𝑘 = 𝐴 ⃗⃗⃗⃗⃗ 𝑘 = 1,2, … , 𝑛, (1.40) тогда равенство (1.39) приобретает вид 𝑦 ⃗⃗⃗⃗⃗ ⃗𝑛−1 − 𝑝2 𝑦⃗𝑛−2 − ⋯ − 𝑝𝑛−1 𝑦⃗1 − 𝑝𝑛 𝑦⃗0 = ⃗0⃗, 𝑛 − 𝑝1 𝑦 (1.41) 𝑦1𝑘 ⋮ ) k = 0, 1, …, n. где 𝑦 ⃗⃗⃗⃗⃗ 𝑘 =( 𝑦𝑛𝑘 Векторы 𝑦 ⃗⃗⃗⃗⃗ 𝑘 удобно находить с помощью рекуррентной формулы 𝑦 ⃗⃗⃗⃗⃗ ⃗𝑘−1 𝑘 = 𝐴𝑦 k = 0, 1, …, n. (1.42) Следовательно, векторное равенство (1.41) эквивалентно системе линейных алгебраических уравнений ⃗⃗ i= 1, 2,…,n 𝑝1 𝑦⃗𝑖,𝑛−1 − 𝑝2 𝑦⃗𝑖,𝑛−2 − ⋯ − 𝑝𝑛−1 𝑦⃗𝑖1 − 𝑝𝑛 𝑦⃗𝑖0 = 0 (1.43) из которой можно найти неизвестные значений p1, p2, …, pn, определяющие коэффициенты характеристического уравнения. Если векторы 𝑦 ⃗⃗⃗⃗⃗ 𝑘 k = 0, 1, …, n-1 линейно независимы, то система (1.43) будет иметь единственное решение, если эти векторы окажутся линейно зависимы, то можно изменить начальный вектор 𝑦 ⃗⃗⃗⃗⃗. 0 Пусть 𝜆1 , 𝜆2 , … , 𝜆𝑛 – корни характеристического уравнения. Разложим вектор 𝑦 ⃗⃗⃗⃗⃗, 0 использовавшийся при вычислении собственных значений, по собственным векторам ⃗⃗⃗⃗, 𝑥1 ⃗⃗⃗⃗⃗, 𝑥2 … , ⃗⃗⃗⃗⃗ 𝑥𝑛 матрицы А 𝑦0 = с1𝑥 ⃗⃗⃗⃗⃗ ⃗⃗⃗⃗1 + с2⃗⃗⃗⃗⃗ 𝑥2 + ⋯ + с𝑛 ⃗⃗⃗⃗⃗ 𝑥𝑛 (1.44) 16

где 𝑐1 ≠ 0 i = 1, 2, …, n – некоторые коэффициенты. Отсюда, учитывая, что 𝐴𝑥⃗𝑖 = 𝜆𝑖 𝑥⃗𝑖 𝐴2𝑥⃗𝑖 = 𝜆2𝑖 𝑥⃗𝑖 … 𝐴𝑛−1 𝑥⃗𝑖 = 𝜆𝑛−1 𝑥⃗𝑖 𝑖 (1.45) i = 1, 2, …, n получим 𝑦⃗1 = 𝑐1 𝜆1𝑥⃗1 + 𝑐2 𝜆2 𝑥⃗2 + ⋯ + 𝑐𝑛 𝜆𝑛 𝑥⃗𝑛 𝑦⃗2 = 𝑐1𝜆21 𝑥⃗1 + 𝑐2𝜆22 𝑥⃗2 + ⋯ + 𝑐𝑛 𝜆2𝑛 𝑥⃗𝑛 (1.46) 𝑦⃗𝑛−1 = 𝑐1𝜆𝑛−1 ⃗1 + 𝑐2 𝜆𝑛−1 ⃗2 + ⋯ + 𝑐𝑛 𝜆𝑛−1 ⃗𝑛 𝑛 𝑥 1 𝑥 2 𝑥 Пусть 𝑄𝑖 (𝜆) = 𝜆𝑛−1 + 𝑞1𝑖 𝜆𝑛−2 + ⋯ + 𝑞𝑛−1,𝑙 i = 1, 2, …, n (1.47) произвольная система многочленов. Составляя линейную комбинацию векторов 𝑦⃗𝑛−1, 𝑦⃗𝑛−2, …,𝑦⃗0 с коэффициентами из (1.47) в силу соотношений (1.44) и (1.46) находим 𝑦⃗𝑛−1 + 𝑞1𝑖 𝑦⃗𝑛−2 + ⋯ + 𝑞𝑛−1,𝑖 𝑦⃗0 = = 𝑐1𝑄𝑖 (𝜆1)𝑥⃗1 + 𝑐2𝑄𝑖 (𝜆2)𝑥⃗2 + ⋯ + 𝑐𝑛 𝑄1 (𝜆𝑛 )𝑥⃗𝑛 (1.48) i = 1, 2, …, n Если положить 𝑄𝑖 (𝜆) = 𝐷(𝜆) 𝜆−𝜆𝑖 i = 1, 2, …, n (1.49) 17

то, очевидно 𝑄𝑖 (𝜆𝑗 ) = { 0, 𝑗 ≠ 𝑖 𝐷`(𝜆𝑖 ) ≠ 0, 𝑗 = 𝑖 (1.50) Формула (1.48) при этом принимает вид 𝑐𝑖 𝑄𝑖 (𝜆𝑖 )𝑥⃗𝑖 = 𝑦⃗𝑛−1 + 𝑞1𝑖 𝑦⃗𝑛−2 + ⋯ + 𝑞𝑛−1,𝑖 𝑦⃗0 i = 1, 2, …, n (1.51) Коэффициенты qji j = 1, 2, …, n-1 могут быть легко определены по схеме Горнера 𝑞0𝑖 = 1 , 𝑞𝑗𝑖 = 𝜆𝑖 𝑞𝑗−1,𝑖 − 𝑝𝑗 𝑗 = 1, 2, … , 𝑛 − 1 (1.52) Таким образом формула (1.51), в которой коэффициенты вычисляются по формуле (1.52), определяет собственный вектор 𝑥⃗𝑖 , отвечающий собственному значению 𝜆𝑖 , с точностью до числового множителя. 1.3.5 Метод Леверрье-Фаддеева Этот метод получения характеристического уравнения матрицы основан на формулах Ньютона для сумм степеней корней алгебраического уравнения. Пусть 𝐷 (𝜆) = 𝜆𝑛 − 𝑝1 𝜆𝑛−1 − 𝑝2 𝜆𝑛−2 − ⋯ − 𝑝𝑛 характеристический многочлен матрицы A и 𝜆1, 𝜆2 , … , 𝜆𝑛 (1.53) – полная совокупность его корней, где каждый корень повторяется столько раз, какова его кратность. 18

Обозначим 𝑆𝑘 = 𝜆𝑘1 + 𝜆𝑘2 + ⋯ + 𝜆𝑘𝑛 , k = 1, 2, …, n. Тогда при k ≤ n справедливы формулы Ньютона: 𝑘𝑝𝑘 = 𝑆𝑘 − 𝑝1 𝑆𝑘−1 − ⋯ − 𝑝𝑘−1𝑆1 , 𝑘 = 1,2, … , 𝑛 (1.54) Если числа Sk известны, то, решая рекуррентную систему (1.54), можно найти нужные коэффициенты pk: 𝑝1 = 𝑆1 1 { 𝑝2 = − (𝑆2 − 𝑝1 𝑆1) 2 … 1 𝑝𝑛 = − (𝑆𝑛 − 𝑝1 𝑆𝑛−1 − ⋯ − 𝑝𝑛−1𝑆1 ) (1.55) 𝑛 Из формул (1.45) следует, что числа 𝜆𝑘𝑖 , I = 1, 2, …, n являются собственными значениями матрицы Ak, а из формул вытекает, что 𝑆𝑘 = 𝑆𝑝𝐴𝑘 , k = 1, 2, …, n, (1.56) где 𝑆𝑝𝐴𝑘 – сумма элементов главной диагонали (след) матрицы Ak. Степени 𝐴𝑘 = 𝐴𝑘−1 ∙ 𝐴 находятся непосредственным перемножением. Таким образом, схема развертывания определителя по методу Леверрье состоит в следующем: сначала вычисляются Аk, k = 1, 2, …, n – степени данной матрицы А, затем находятся соответствующие Sk – суммы элементов главных диагоналей матриц Ak и, наконец, по формулам (1.55) определяются искомые коэффициенты pk, k = 1, 2, …, n. Д.К.Фаддеев предложил видоизменение метода Леверрье, которое кроме упрощений при вычислении коэффициентов характеристического многочлена позволяет определить обратную матрицу и собственные векторы матрицы. 19

Будем вместо последовательности A, A2, …, вычислять An последовательность A1, А2, …, An, построенную следующим образом: 𝐴1 = 𝐴, 𝑆𝑝𝐴1 = 𝑞1 𝐵1 = 𝐴1 − 𝑞1 ∙ 𝐸 𝐴2 = 𝐴𝐵1 𝑆𝑝𝐴2 = 𝑞2 2 𝐵2 = 𝐴2 − 𝑞2 ∙ 𝐸 … … 𝐴𝑛−1 = 𝐴𝐵𝑛−2 𝐴𝑛 = 𝐴𝐵𝑛−1 … 𝑆𝑝𝐴𝑛−1 = 𝑞𝑛−1 𝑛−1 𝑆𝑝𝐴𝑛 = 𝑞𝑛 𝑛 (1.57) 𝐵𝑛−1 = 𝐴𝑛−1 − 𝑞𝑛−1 ∙ 𝐸 𝐵𝑛 = 𝐴𝑛 − 𝑞𝑛 ∙ 𝐸 где E – единичная матрицы того же порядка, что и матрица A. Ниже приведены доказательства следующих утверждений: а) q1=p1, q2=p2, …, qn=pn; б) Bn – нулевая матрица; в) если A – неособенная матрица, то 𝐴−1 = 𝐵𝑛−1 𝑝𝑛 . а) Используем метод математической индукции. Пусть n = 1, тогда 𝑝1 = 𝑆𝑝𝐴 = 𝑞1 . Предположим, что при n = k q1=p1, q2=p2, …, qn=pn, и возьмем n = k+1. Согласно (1.57) 𝐴𝑘+1 = 𝐴𝑘+1 − 𝑞1 𝐴𝑘 − ⋯ − 𝑞𝑘 𝐴 = 𝐴𝑘+1 − 𝑝1𝐴𝑘 − ⋯ − 𝑝𝑘 𝐴. (1.58) Следовательно 𝑆𝑝𝐴𝑘+1 = (𝑘 + 1)𝑞𝑘+1 = 𝑆𝑝𝐴𝑘+1 − 𝑝1 𝑆𝑝𝐴𝑘 − ⋯ − 𝑝𝑘 𝑆𝑝𝐴 = (1.59) = 𝑆𝑘+1 − 𝑝1 𝑆𝑘 − ⋯ − 𝑝𝑘 𝑆1 . 20

Отсюда в силу формул Ньютона (𝑘 + 1)𝑞𝑘+1 = (𝑘 + 1)𝑝𝑘+1 и, следовательно, 𝑞𝑘+1 = 𝑝𝑘+1, что доказывает а). б) В силу тождества Гамильтона-Кели 𝐵𝑛 = 𝐴𝑛 − 𝑝1 𝐴𝑛−1 − ⋯ − 𝑝1 𝐸 = 0 (1.60) в) Из формул (1.59) и (1.60) следует, что 𝐴𝐵𝑛−1 = 𝐴𝑛 = 𝐵𝑛 + 𝑝𝑛 𝐸 = 𝑝𝑛 𝐸 поэтому 𝐴−1 = Таким 𝐵𝑛−1 𝑝𝑛 (1.61) . образом коэффициенты характеристического многочлена матрицы А определяются с помощью формул (1.57). Посмотрим матрицу 𝑄𝑖 = 𝜆𝑛−1 𝐸 + 𝜆𝑛−2 𝐵1 + ⋯ + 𝜆𝑖 𝐵𝑛−2 + 𝐵𝑛−1, 𝑖 𝑖 (1.62) где 𝐵𝑘 – матрицы, вычисленные по формулам (1.57), а 𝜆𝑖 – i-е собственное значение матрицы А. Можно доказать, в предположении, что все 𝜆1 , 𝜆2 , … , 𝜆𝑛 различны, что матрица Qi ненулевая. Покажем, что каждый столбец матрицы Qi, состоит из компонент собственного вектора, отвечающего собственному значению 𝜆𝑖 . Действительно, (𝜆𝑖 𝐸 − 𝐴)𝑄𝑖 = (𝜆𝑖 𝐸 − 𝐴)(𝜆𝑛−1 𝐸 + 𝜆𝑛−2 𝐵1 + 𝜆𝑛−3 𝐵2 + ⋯ + 𝜆𝑖 𝐵𝑛−2 + 𝐵𝑛−1) = 𝑖 𝑖 𝑖 = 𝜆𝑛𝑖 𝐸 + 𝜆𝑛−1 𝐵1 + 𝜆𝑛−2 𝐵2 + ⋯ + 𝜆2𝑖 𝐵𝑛−2 + 𝜆𝑖 𝐵𝑛−1 − 𝜆𝑛−1 𝐴 − 𝜆𝑛−2 𝐴𝐵1 − ⋯ − 𝑖 𝑖 𝑖 𝑖 (𝐵1 − 𝐴) + −𝜆𝑖 𝐴𝐵𝑛−2 − 𝐴𝐵𝑛−1 = 𝜆𝑛𝑖 𝐸 + 𝜆𝑛−1 𝑖 (1.63) (𝐵2 − 𝐴𝐵1 ) + ⋯ + 𝜆𝑖 (𝐵𝑛−1 − 𝐴𝐵𝑛−2) − 𝐴𝐵𝑛−1 = +𝜆𝑛−2 𝑖 = 𝜆𝑛𝑖 𝐸 − 𝑝1 𝜆𝑛−1 𝐸 − 𝑝2 𝜆𝑛−2 𝐸 − ⋯ − 𝑝𝑛−1𝜆𝑖 𝐸 − 𝑝𝑛 𝐸 = 0 𝑖 𝑖 21

Отсюда следует, что для любого столбца 𝑦⃗ построенной матрицы Qi (𝜆𝑖 𝐸 − 𝐴) ∙ 𝑦⃗ = ⃗0⃗, т.е. 𝑦⃗ – собственный вектор матрицы А, отвечающий собственному значению 𝜆𝑖 . Вычисляя собственные векторы описанным образом, нет необходимости находить все столбцы матрицы Qi. Следует ограничиться вычислением одного (любого) столбца, для чего удобно пользоваться рекуррентной формулой: 𝑦⃗0 = 𝑒⃗, 𝑦⃗𝑘 = 𝜆𝑖 𝑦⃗𝑘−1 + 𝑏⃗⃗𝑘 𝑘 = 1,2, … , 𝑛 − 1, (1.64) где 𝑏⃗⃗𝑘 – одноименный с вычисляемым столбцом матрицы Qi столбец матрицы Bk, а 𝑒⃗ – одноименный столбец единичной матрицы. Тогда собственный вектор 𝑥⃗𝑖 матрицы А, отвечающий собственному значению 𝜆𝑖 , есть 𝑥⃗𝑖 = 𝑦⃗𝑛−1. 1.4 Обзор технологии NVIDIA CUDA CUDA – это архитектура параллельных вычислений от NVIDIA, позволяющая существенно увеличить вычислительную производительность благодаря использованию GPU (графических процессоров) [31]. В состав CUDA, кроме самой видеокарты, входят программные компоненты, объединённые в пакет CUDA Toolkit. В ранних версиях так же требовался дополнительный драйвер CUDA, но в настоящее время он включен в стандартные пакет драйверов для видеокарты NVIDIA. CUDA Toolkit содержит компилятор nvcc, транслирующий исходный код программ в промежуточный ассемблерный код, библиотеки, необходимые для работы с платформой, а так же дополнительное программного обеспечение, например Nsight Monitor. Многие разработчики программного обеспечения, ученые и исследователи широко применяют CUDA в различных областях: физика 22

и химия, обработка видео и изображений, в медицине. С помощью нее моделируют динамику жидкостей, восстанавливают изображения, которые были получены при помощи компьютерной томографии, проводят сейсмический анализ, трассировку лучшей и многое другое. Платформа параллельных вычислений CUDA обеспечивает набор расширений для языков С и С++, позволяющих выражать как параллелизм данных, так и параллелизм задач на уровне мелких и крупных структурных единиц. Такой подход позволяет писать код, исполняемый на GPU, так и код, выполняющийся на CPU. Кроме того, программист может выбирать между средствами разработки: использовать высокоуровневые языки или же открытые стандарты, такие как директивы OpenACC. 1.4.1 Модель памяти CUDA В спецификациях к архитектуре CUDA выделяют 6 видов памяти:  Регистровая;  Разделяемая;  Локальная;  Глобальная;  Константная;  Текстурная. Регистровая память представляет собой сверхоперативное запоминающее устройство в мультипроцессоре. Данный вид памяти располагается на мультипроцессоре, не кэшируется и обладает максимальной скоростью среди других видов. Для большей эффективности надо стараться занимать как можно меньше регистров, т.к. на каждый поток выделяется некоторое количество регистров и при их нехватке количество потоков будет ограничено. 23

Константная память физически не отделена от глобальной памяти. Но в отличии от глобальной может кэшироваться и в таком случае скорость работы с данными достаточно высока. Однако, если необходимые данные отсутствуют в кэше, то их чтение будет выполнено с задержками в 400-600 тактов. Локальная память занимает часть DRAM (динамическая память с произвольным доступом) и используется, когда локальные данные процедур занимают слишком большой размер. Этот вид памяти так же не кэшируется и имеет низкую скорость доступа. Глобальная память занимают основную часть DRAM. Обладает высокой пропускной способностью (более 100Гб/с) и возможностью произвольной адресации глобальной памяти. Но также, как и локальная память не кэшируется и обладает низкой скоростью работы. Разделяемая память обладает скоростью, сравнимой с регистровой памятью. Основное назначение — это обеспечение взаимодействия между потоками, в связи с чем эта память доступна на запись и на чтение из всех потоковых процессоров в мультипроцессоре. Текстурная память – блок памяти, доступный только на чтение всеми мультипроцессорами. Выборка из этого типа памяти происходит при помощи текстурных блоков видеочипа. За счёт этого возможно выполнение линейной интерполяции без каких-либо дополнительных затрат. 1.4.2 Математические библиотеки CUDA В состав CUDA Toolkit входит набор таких математических библиотек, как cuBLAS, cuSPARSE и другие [30]. Общей особенностью этих пакетов является то, что все они предназначены для произведения параллельных вычислений на GPU. 24

CuBLAS представляет собой реализацию BLAS (базовых подпрограмм линейной алгебры) поверх среды CUDA. Для использования API cuBLAS программист должен предварительно выделить и заполнить память видеокарты для матриц, векторов, а также, при необходимости, других возможных параметров функций. Ранние версии библиотеки имели собственные функции для выделения памяти и копирования данных, но в последних версиях платформы от данной практики было решено отказаться и теперь для этих целей используются общие функции CUDA. Ещё одной особенностью cuBLAS является индексация массивов. Для максимальной совместимости с существующими средами Fortran библиотека cuBLAS использует индексирование массивов, начинающее с единицы. Но, т.к. в C и C++ используется индексация с 0, то для совместимости с этими языками в некоторые функции библиотеки передаётся дополнительный параметр, который указывает, с какой позиции начинается индексация. В отличии от cuBLAS, библиотека cuSPARSE предназначена для обработки разреженных матриц. Библиотечные функции можно разделить на четыре категории:  Операции между вектором в разреженном формате и вектором в плотном формате;  Операции между матрицей в разреженном формате и вектором в плотном формате;  Операции между матрицей в разреженном формате и набором векторов в плотном формате;  Операции, которые допускают преобразование между различными матричными форматами. 25

1.5 Формат Compressed Sparse Row Формат CSR предназначен для компактного хранения разреженных матриц [16]. Он хранит разреженную матрицу М размером m x n в виде трех одномерных массивов:  Массив A, который имеет длину, равную количеству ненулевых элементов матрицы М и, соответственно, хранящий все не нулевые элементы в порядке слева направо и сверху вниз;  Массив JA, размерности аналогичной массиву А, хранит номера столбцов каждого элемента А;  Массив IA, размером m + 1, в котором i-й элемент указывает, с какой позиции в массивах A и JA начинается i-я строка. Для примера возьмем следующую разреженную матрицу: 1 0 𝐴 =( 0 7 0 3 0 0 0 5 6 0 2 0 ) 8 1 Для нее m = 4 и массивы имеют вид: A = [1, 2, 3, 4, 6, 8, 7, 1] JA = [1, 4, 2, 3, 3, 4, 1, 4] IA = [1, 3, 5, 7, 9] Данный формат имеет следующий преимущества: требует меньше памяти, по сравнению с классическим вариантом хранения матрицы в памяти, обеспечивает быстрый доступ к строке и умножение матрицы на вектор. 26

2. РАЗРАБОТКА АЛГОРИТМА РЕШЕНИЯ ОБОБЩЕННОЙ ЗАДАЧИ ПОИСКА СОБСТВЕННЫХ ЗНАЧЕНИЙ И ВЕКТОРОВ ДЛЯ РАЗРЕЖЕННЫХ МАТРИЦ НА ОСНОВЕ LU-РАЗЛОЖЕНИЯ ПРАВОЙ ЧАСТИ УРАВНЕНИЯ 2.1 Преобразование обобщенной задачи поиска собственных значений и векторов к классической задаче Данный алгоритм применяется для решения обобщенной задачи поиска собственных значений вида 𝐴𝑥 = 𝜆𝐵𝑥 (2.1) где матрицы A и B не вырожденные общего вида (могут быть не положительно определенными, не симметричными), λ – собственное значение, x – собственный вектор. Для решения этой задачи, необходимо свести её к задаче вида 𝐶𝑥 = 𝜆𝑥 (2.2) В качестве первого шага, используем LU разложение матрицы B 𝐵 = 𝐿𝑈 (2.3) 𝐴𝑥 = 𝜆𝐿𝑈𝑥 (2.4) Тогда Произведем замену произведения матрицы U на вектор x: 27

𝑦 = 𝑈𝑥 (2.5) 𝑥 = 𝑈 −1𝑦 (2.6) 𝐴𝑈 −1 𝑦 = 𝜆𝐿𝑦 (2.7) Тогда Умножим обе части уравнения слева на матрицу обратную L 𝐿−1𝐴𝑈 −1𝑦 = 𝜆𝑦 (2.8) Перемножим матрицы L-1, A и U-1 и обозначим результат как матрицу С: 𝐶 = 𝐿−1𝐴𝑈 −1 (2.9) Тогда исходное уравнение сведется к следующему: 𝐶𝑦 = 𝜆𝑦 (2.10) откуда, зная y, можно вычислить собственный вектор x умножив матрицу обратную U на y. 2.2 Использование неполного LU-разложения матриц Одним из широко известных способов разложения матриц на множители является LU-факторизация, позволяющая представить матрицу А в виде 28

𝐴 = 𝐿𝑈, (2.11) где L и U – нижне- и верхнетреугольные матрицы соответственно. Однако, алгоритм факторизации непригоден для разреженных матриц, так как ведет к заполнению портрета, т.е. появлению в матрицах L и U ненулевых элементов в тех позициях, для которых aij = 0, - и как следствие, резкому увеличению объёма памяти, требуемой для хранения матриц. Поэтому, вместо задачи нахождения факторизации (2.11) сформулируем другую задачу. Для заданной матрицы А потребуем представить её в виде 𝐴 = 𝐿𝑈 + 𝑅, (2.12) где матрицы в правой части удовлетворяют следующим свойствам:  Матрицы L и U являются нижнетреугольной и верхнетреугольной соответственно;  𝑃𝐿 ⊂ 𝑃𝐴 и 𝑃𝑈 ⊂ 𝑃𝐴 ;  𝑃𝐴 ∩ 𝑃𝑅 = ∅ . Тогда приближенное предстваление 𝐴 ≈ 𝐿𝑈 называется неполной LUфакторизацией матрицы A или коротко её ILU-разложением. Использование ILU – разложение на CSR матрицей позволит эффективно использовать память. Т.к. портрет матрицы не изменится, то не потребуется выделять дополнительную память для хранения новой матрицы, а результат можно будет записать на место исходной матрицы. Массивы, содержащие номера столбцов элементов и позиции, с которых начинаются новые строки, не потребует изменения. Таким образов будет обновлен только массив, который содержит ненулевые элементы. Алгоритм ILU-разложения представлен на рисунке 2.1. 29

Начало A, n k=0; k < n; k++ j=0; j < k; j++ A[k,j]=0 L[k,j]=0 s=0 i=0; i < j; i++ s+=L[k,i]*U[i,j] L[k,j]=(A[k,j]-s)/U[I,j] L[k,k]=1 j=k-1; j < n; j++ A[k,j]=0 U[k,j]=0 s=0 i=0; i < k-1; i++ s+=L[k,i]*U[i,j] U[k,j]=A[k,j]-s L, U Конец Рисунок 2.1. Алгоритм ILU-разложения 30

2.3 Степенной метод со сдвигами нахождения собственных значений Степенной метод позволяет найти наибольшее по модулю собственное значение и собственный вектор матрицы А. Основным достоинством данного метода является возможность эффективного использования разреженной матрицы, т.к. векторы получаются только с помощью умножения матрицы на вектор, а исходная матрица остаётся разреженной даже после выполнения сдвига. В качестве нулевого приближения берут произвольный вектор x(0) и последовательно строят следующие приближения: 𝑥 (1) = 𝐴 ∙ 𝑥 (0) 𝑥 (2) = 𝐴 ∙ 𝑥 (1) (2.13) … 𝑥 (𝑘) = 𝐴 ∙ 𝑥 (𝑘−1) Для соответствующих координат векторов данной последовательности буду выполняться выражение lim 𝑘→∞ 𝑥 (𝑘) 𝑥 (𝑘−1) = 𝜆(𝑘) (2.14) Для того, чтобы степенной метод можно было использовать для нахождения не только максимального собственного значения применяется следующее преобразование 𝐴′ = 𝐴 − 𝜎𝐸 (2.15) где А – исходная матрица, Е – единичная матрица, 𝜎 – некоторое число. 31

Обозначим 𝜆∗ как следующее собственное значение после 𝜆. Тогда, если 𝜎 = 𝜆, то число 𝜆 − 𝜆∗ станет минимальным по модулю, а максимальным по модулю окажется сдвиг другого собственного значения. 32

3. ПРОГРАММНАЯ РЕАЛИЗАЦИЯ РАЗРАБОТАННОГО АЛГОРИТМА С ИСПОЛЬЗОВАНИЕМ ТЕХНОЛОГИИ CUDA 3.1 Программная реализация параллельной версии алгоритма Для изоляции вычисления собственных значений и векторов от остального кода, все с этим связанные операции удобно вынести в отдельную функцию. Данная функция могла бы иметь только следующие входные значения:  Матрица A в CSR формате;  Матрица B в CSR формате;  Точность вычислений. Однако, ввиду особенностей языка программирования C++, к этим данным необходимо добавить следующие:  Размерность матрицы;  Количество ненулевых элементов в матрице А;  Количество ненулевых элементов в матрице B;  Ссылка на массив, куда будут записаны вычисленные собственные значения;  Ссылка на двумерный массив, куда будут записаны вычисленные собственные вектора. Так как вычисления будут производиться на GPU, то данные, над которыми будут проводиться операции необходимо скопировать в память видеокарты, предварительно выделив для них достаточное количество места. Пример выделения памяти для массива ненулевых элементов матрицы A представлен на листинге 3.1. Листинг 3.1. Выделение памяти для матрицы double *devA; cudaStat = cudaMalloc((void**)&devA, iNNZA * sizeof(double)); if (cudaStat != cudaSuccess) return -1; 33

Пример копирования данных из оперативной памяти в память GPU представлен на рисунке 3.2. Листинг 3.2. Копирование данных из ОЗУ и в память ГПУ cudaStat = cudaMemcpy(devA, a, iNNZA * sizeof(double), cudaMemcpyHostToDevice); if ((cudaStat != cudaSuccess)) { return -2; } Согласно алгоритму, первым этапом вычислений является ILU разложение. В библиотеке cuSPARSE, которая входит в состав Cuda Toolkit, есть специальный метод, который произведет необходимую операцию над данными. Однако перед тем, как передать в него данные, необходимо выделить достаточно памяти буферной памяти. Для этого предусмотрены специальные функции (листинг 3.3), куда необходимо передать информация о размере матрицы, количестве ненулевых элементов матрицы B и ссылку, куда будут записаны прогнозируемые размеры буферов. Листинг 3.3. Выделение буфера для ILU-разложения cusparseDcsrilu02_bufferSize(cusparseHandle, iEqNum, iNNZB, descr_B, devB, devIBcsr, devJB, info_B, &pBufferSize_M); cusparseDcsrsv2_bufferSize(cusparseHandle, trans_L, iEqNum, iNNZB, descr_L, devB, devIBcsr, devJB, info_L, &pBufferSize_L); cusparseDcsrsv2_bufferSize(cusparseHandle, trans_U, iEqNum, iNNZB, descr_U, devB, devIBcsr, devJB, info_U, &pBufferSize_U); pBufferSize = std::max(pBufferSize_M, std::max(pBufferSize_L, pBufferSize_U)); сudaMalloc((void**)&pBuffer, pBufferSize); После того как память для буфера была выделена, вызывается функция, которая произведет неполное lu разложение. Матрица L и матрица U будут записаны на место матрицы B. Таким образом, кроме временно выделенного места для буфера, больше память не расходуется. 34

Следующий этап заключается в вычислении матрицы С, которая получается в результате произведения матрицы обратной L на матрицу A и на матрицу обратную U. Однако, предварительно требуется получить обратные матрицы. После чего необходимо произвести две операции перемножения матриц. В библиотеке cuSPARSE для этих целей предусмотрены два метода: cusparseXcsrgemmNnz и cusparseScsrgemm. Что бы осуществить умножение необходимо использовать их последовательно. Так cusparseXcsrgemmNnz заполнит массив, который содержит позиции начала строк матрицы в CSR формате и количество ненулевых элементов. А метод cusparseScsrgemm основываясь на этих данных, заполнит два оставшихся массива, которые хранят ненулевые элементы и их номера столбцов: Листинг 3.4. Умножение матриц cusparseXcsrgemmNnz(cusparseHandle, trans_L, trans_L, iEqNum, iEqNum, iEqNum, descr_Linv, iNNZB, devIBcsr, devJB, descr_A, iNNZA, devIAcsr, devJA, descr_C, devICcsr, nnzTotalDevHostPtr); cusparseDcsrgemm(cusparseHandle, trans_L, trans_L, iEqNum, iEqNum, iEqNum, descr_Linv, iNNZB, devLinv, devIBcsr, devJB, descr_A, iNNZA, devA, devIAcsr, devJA, descr_C, devC, devICcsr, devJC); Данная операция производится два раза. В первый раз умножая матрицу А слева на матрицу обратную L, а во второй, матрицу- результат предыдущей операции на матрицу обратную U. В результате чего в память будет помещена матрица C, которая является результатом вычислений. Следующим этапом является вычисление собственных значений для матрицы С. Метод для решения данной задачи уже предусмотрен в библиотеке cuSOLVER, но для того, чтобы им воспользоваться, требуется вычислить приближенные собственные значения. Воспользуемся степенным методом со сдвигом. В начале выделяется память для массива собственных векторов во время работы метода и заполняется единицами (листинг 3.5). 35

Листинг 3.5. Выделение памяти для массива собственных векторов и заполнение его единицами double *devX; cudaMalloc((void**)&devX, iEqNum * sizeof(double)); setXone << <(iEqNum - 1 + BLOCK_SIZE) / BLOCK_SIZE, BLOCK_SIZE >> >(devX, iEqNum); cudaDeviceSynchronize(); Теперь в цикле выполняется вычисление собственных значений. Для этого сначала умножается матрица А на приближенный собственный вектор, который на первой итерации заполнен единицами, и записывается результат как собственный вектор. Листинг 3.6. Умножение матрицы на вектор cusparseDcsrmv(cusparseHandle, trans_L, iEqNum, iEqNum, nnzA, &alpha, descr_A, devAA, devIAAcsr, devJAA, devX, &gamma, devX); После этого выполняется нормировка вектора (листинг 3.7), это необходимо для того, чтобы не допустить выхода за диапазон типа данных double. Листинг 3.7. Нормировка вектора cublasDnrm2(cublasHandle, iEqNum, devX, 1.0, &normX); normX = 1.0 / normX; cublasDscal(cublasHandle, iEqNum, &normX, devX, 1.0); Далее повторно выполним аналогичную операцию умножения и запишем результат в другой массив. Проведем скалярное умножение (листинг 3.8) полученного результата с ранее вычисленным приближенным значением собственного вектора. В итоге мы получим приближенное собственное значение: 36

Листинг 3.8. Скалярное умножение векторов cublasDdot(cublasHandle, iEqNum, devC, 1.0, devX, 1.0, &lambda[k]); После того, как работа цикла будет завершена, необходимо произвести сдвиг, чтобы перейти к вычислению следующего собственного значения. Когда будут найдены все приближенные собственные значения, можно воспользоваться функцией cusolverSpDcsreigvsi, которая вычислит собственное значение и собственный вектор с заданной точностью. Данный метод одновременно вычисляет значения только одного собственного значения и собственного вектора, поэтому её необходимо вызывать в цикле, передавая на вход каждый раз следующее приближенное собственное значение, вычисленное ранее: Листинг 3.9. Нахождение собственных значений и векторов с помощью функции cusolverSpDcsreigvsi for (int k = 0; k < iEqNum; k++) { cusolverSpDcsreigvsi(cusolverHandle, iEqNum, NNZA, descr_A, devA, devIAcsr, devJA, lambda[k], devX0, maxite, eps, devMu, devX); cudaMemcpy(egval + k, devMu, sizeof(double), cudaMemcpyDeviceToHost); cudaMemcpy(eigvec[k], devX, iEqNum * sizeof(double), cudaMemcpyDeviceToHost); } По окончанию работы цикла, достаточно освободить память от временных переменных. Собственные значения и вектора будут доступны по ссылкам, которые передавались из управляющего метода. 37

3.2 Оценка эффективности разработанного алгоритма на основе проведения вычислительных экспериментов. Для проведения оценки эффективности проводился ряд испытаний при вычислении собственных значений и векторов матриц различной размерности. Ввиду того, что персональный компьютер, на котором производятся испытания, не обладает достаточной мощностью, чтобы обрабатывать большие объёмы данных (Intel Core i5-4210u, Nvidia GeForce gt840m 2gb), в качестве оптимальных и достаточных для демонстрации эффективности были выбраны следующие размеры матриц: 500, 1000, 2000, 3000 и 4000. Все вычисления производились над различными наборами данных и для каждой размерности проводилось по 3 испытания для того чтобы получить в результате более объективные результаты. Результаты, полученные в результате испытаний, приведены в таблице 3.1. Все значения были округлены до целых секунд. Таблица 3.1 Результаты испытаний параллельной версии алгоритма Порядок Испытание Испытание Испытание Среднее время матрицы №1, с №2, с №3, с 3 испытаний, с 500 10 9 10 10 1000 35 31 29 32 2000 189 186 193 189 3000 420 431 433 428 4000 644 631 633 636 Для наглядной демонстрации эффективности работы программы, наиболее рациональным способом было выбрано сравнение с 38

последовательной версией программы. Были проведены аналогичные испытания, результаты которых приведены в таблице 3.2 Таблица 3.2 Результаты испытаний последовательной версии алгоритма Порядок Испытание Испытание Испытание Среднее время матрицы №1 №2 №3 3 испытаний 500 4 4 4 4 1000 42 36 44 41 2000 382 394 384 386 3000 1328 1300 1318 1316 4000 3520 3543 3537 3533 Для более удобной демонстрации результатов испытаний параллельной и последовательной версии, они представлены в виде графика на рисунке 3.1. Таблица 3.3 Объединенные результаты испытаний программ 500 1000 2000 3000 4000 Параллельная 10с 32с 189с 428с 636c Последовательная 4с 41с 386с 1316с 3533c 39

Рисунок 3.1. Сравнительный график скорости выполнения последовательной и параллельной версии программы Для получения значения ускорения воспользуемся следующей формулой: 𝑆𝑝 (𝑛) = 𝑇1 (𝑛) 𝑇𝑝 (𝑛) , (3.1) где n – размерность матрицы, 𝑆𝑝 – ускорение параллельной версии, 𝑇1 – время выполнения последовательной версии, 𝑇𝑝 – время выполнения параллельной версии. Данные об ускорении, полученные по формуле 3.1 представлены в таблице 3.4 и на рисунке 3.2. Таблица 3.4 Ускорения параллельной версии программы Порядок матрицы 500 Ускорение 0,4 1000 1,28 2000 3000 4000 2,04 3,08 5,56 40

Рисунок 3.2. График, демонстрирующий ускорение параллельной версии Исходя из полученных результатов можно сделать вывод, что при вычислении собственных значений и векторов для матриц размером менее 1000 использование параллельной версии не имеет никаких значимых преимуществ. Однако ситуация меняется при размерности матрицы более 2000. Уже при 3000 ускорение составляет 3 единицы, а при 4000 превышает 5. 41

ЗАКЛЮЧЕНИЕ Целью данной выпускной квалификационной работы была разработка и программная реализация итерационного метод решения обобщенной задачи поиска собственных значений и векторов на основе LU-разложения правой матрицы с использованием технологии CUDA. Все поставленные задачи были достигнуты, а именно:  Произведен обзор существующих методов решения задачи поиска собственных значений и векторов и сделан вывод, что большинство существующих методов предназначены только для нахождения собственных значений симметричный и положительно определенных матриц, а остальные не эффективно расходуют память;  Разработан алгоритм решения обобщенной задачи поиска собственных значений и векторов для разреженных матриц на основе LU-разложения правой части уравнения, который может решать несимметричные и не положительно определенные матрицы;  Выполнена программная реализация разработанного алгоритма с использованием технологии CUDA, в которой используется формат CSR для хранения матриц и минимально расходуется память;  Осуществлена оценка эффективности разработанного алгоритма на основе проведения вычислительных экспериментов, по результатам которых сделан вывод, что параллельный алгоритм позволяет ускорить вычисления в 5,5 раз. По результатам проведенных вычислительных экспериментов можно сделать вывод, что параллельная версия реализованного приложения выполняет задачу поиска собственных значений и векторов быстрее последовательной версии при размерности матриц более 1000. Кроме того, во время выполнения расходуется минимальное количество памяти, за счет чего можно осуществлять вычисления над большими объемами данными, по 42

сравнению с существующими реализациями подобных задач. Исходя из этого, можно заключить, что разработанное приложение можно использоваться для решения различных прикладных задач. 43

СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ 1. NVIDIA Developer Documentation. URL: https://docs.nvidia.com/ 2. Абаффи Й. Математические методы для линейных и нелинейных уравнений: проекционные ABS-алгоритмы / Й. Абаффи, Э. Спедикато. — Москва: Мир, 1996. 3. Баландин М.Ю. Методы решения СЛАУ большой размерности / М.Ю. Баландин, Э.П. Шурина. - Новосибирск: НГТУ, 2000. — 70 с. 4. Бахвалов Н.С. Практикум по алгебре / Н.С. Бахвалов, А.И. Кострикин. – Москва: МГУ, 1983. 5. Березин И.С. Методы вычислений / И.С. Березин, Н.П. Жидков. – Москва: Физматгиз, 1962. – 635 с. 6. Богачев К.Ю. Практикум на ЭВМ. Методы решения линейных систем и нахождения собственных значений / К.Ю. Богачев. – Москва: МГУ, 1998. – 137с. 7. Боглаев Ю.П. Вычислительная математика и программирование. Москва: Высшая школа, 1990. – 544 с. 8. Бордовицина Т.В. Современные численные методы в задачах небесной механики. – Москва: Наука, 1984. – 136 с. 9. Воеводин В.В. Математические модели и методы в параллельных процессах. – Москва: Наука, 1986. – 259 с. 10. Волков Е.А. Численные методы. – Москва: Наука, 1987. – 248 с. 11. Вычисления на графических процессорах (GPU) в задачах математической и теоретической физики / Перепёлкин Е.Е., Садовников Б.И., Иноземцева Н.Г. – Москва: URSS, 2014. – 176 c. 12. Вычислительные методы для инженеров / А.А. Амосов, Ю.А. Дубинский, Н.В. Копченова — Москва: Высшая школа, 1994. — 544 с. 13. Годунов С.К. Современные аспекты линейной алгебры. — Новосибирск: Научн. книга, 1997. 44

14. Демидович Б.П. Основы вычислительной математики / Б.П. Демидович, И.А. Марон. – Москва: Наука, 1966. – 664 с. 15. Долгополов Д.В. Методы нахождения собственных значений и собственных векторов матриц / Д.В. Долгополов. -Санкт-Петербург: СПбГТИ(ТУ), 2005. – 39с 16. Ильин В.П. Методы неполной факторизации для решения линейных систем. — Москва: Физматлит, 1995. 17. Калиткин Н.Н. Численные методы / Н.Н. Калиткин. -Санкт-Петербург: БХВ-Петербург, 2011. – 592 с. 18. Киреев В.И. Численные методы в примерах и задачах: Учеб. пособие / Киреев В.И., Пантелеев А.В. – Москва: Высшая школа, 2006. – 480 с. 19. Курош А.Г. Курс высшей алгебры. – Москва: Наука, 1975. – 431 с. 20. Марчук Г.И. Методы вычислительной математики, 3-е изд. – Москва: Наука, 1989. 21. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. — Москва: Мир, 1991. 22. Параллельные вычисления на GPU. Архитектура и программная модель CUDA: Учебное пособие / Боресков А.В., Марковский Н.Д., Микушин Д.Н. и др. – Москва: МГУ, 2012. – 336 с. 23. Пирумов У. Г. Численные методы: Учеб. пособие. Москва: Дрофа, 2004. – 224 с. 24. Писсанецки С. Технология разреженных матриц. — Москва: Мир,1988. 25. Плис А.И. Лабораторный практикум по высшей математике / Плис А.И. Сливина Н.А. – Москва: Высшая школа, 1994. – 416 с. 26. Плохотников К.Э. Вычислительные методы. Теория и практика в среде MATLAB (курс лекций). – Москва: Телеком, 2009. – 496 с. 27. Самарский А.А. Введение в численные методы. – Москва: Наука, 1987. – 286 с. 28. Современные информационные технологии в задачах навигации и наведения беспилотных маневренных летательных аппаратов / К.К. 45

Веремеенко, С.Ю. Желтов, Н.В. Ким и др. – Москва: ФИЗМАТЛИТ, 2009. – 562с. 29. Тихонов А.Н. Вводные лекции по прикладной математике / Тихонов А.Н., Костомаров Д.П. – Москва: Наука, 1984. – 160 с. 30. Уилкинсон Дж. Алгебраическая проблема собственных значений / Дж. Уилкинсон. -Москва: Наука, 1970. – 565 с. 31. Фаддеев Д.К. Вычислительные методы линейной алгебры / Д.К. Фаддеев, В.Н. Фаддеева. – Москва: Физматгиз, 1963. – 734с. 32. Хемминг Р.В. Численные методы. – Москва: Наука, 1972. – 400 с. 33. Хорн. Р. Матричный анализ / Р. Хорн, Ч. Джонсон. – Москва: Мир, 1989. 34. Численные методы / Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. – Москва: Наука, 1987. 46

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв