Участник:Elijah/Нахождение собственных чисел квадратной матрицы методом QR разложения

Нахождение собственных чисел квадратной матрицы методом QR разложения
Последовательный алгоритм
Последовательная сложность	[math]N * O(n^3)[/math]
Объём входных данных	[math]n^2[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math] N * (12n - 15) [/math]
Ширина ярусно-параллельной формы	[math] O(n^2) [/math]

Основные авторы описания: И.В.Афанасьев, В.А.Шишватов

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

//TODO

1.2 Математическое описание алгоритма

//TODO

1.3 Вычислительное ядро алгоритма

У базового QR алгоритма есть два вычислительных ядра:

1. операция поиска QR разложения с использованием метода Гивенса
   
2. матричное умножения полученных матриц

// TODO about heisenberg

1.4 Макроструктура алгоритма

Как уже описано в описании ядра алгоритма, базовая версия QR-алгоритма на каждой итерации использует следующие алгоритмы:

1.1. Метод Гивенса (вращений) QR-разложения квадратной матрицы

или

1.2. Метод Хаусхолдера(отражений) QR-разложения_квадратной_матрицы

2. Перемножение плотных неособенных матриц

// TODO about heisenberg

1.5 Схема реализации последовательного алгоритма

Последовательная реализация простейшего алгоритма сострит из некоторого числа итераций.

На итерации [math] n [/math]:

1. Для матрицы [math] A_{n} [/math] строится QR разложение любым доступным последовательным алгоритмом на матрицы [math] Q_{n} [/math] и [math] R_{n} [/math]
   
2. Матрицы [math] R_{n} [/math] и [math] Q_{n} [/math] перемножаются, таким образом получается матрица [math] A_{n+1} = R_{n} * Q_{n} [/math] для следующей итерации [math] n + 1 [/math]

По окончании каждой итерации проверятся, приведена ли матрица к диагональной форме.

// TODO about heisenberg form and algorithm with shift

1.6 Последовательная сложность алгоритма

Рассчитаем последовательную сложность базового алгоритма. Пусть [math] A \in \mathbb{R}^{n \times n}[/math], и для приведения матрицы к диагональной форме необходимо произвести N итераций алгоритма.
На каждой итерации алгоритма производится QR разложение (сложностью [math] 2 * n^3 [/math]) и матричное умножение (сложностью [math] n^3 [/math]). Проверка того, имеет ли матрица диагональную форму, может быть проведена за [math] n^2 [/math] операций.
Таким образом итоговая сложность одной итерации составляет: [math] 3*n^3 + n^2 = O(n^3)[/math]
Общая сложность алгоритма при [math] N [/math] итерациях составляет [math] N * O(n^3) [/math]

// TODO about heisenberg complexity

1.7 Информационный граф

//TODO

1.8 Ресурс параллелизма алгоритма

Все итерации базового QR-алгоритма производятся последовательно, поэтому на верхнем уровне алгоритм чисто последователен.

Основной ресурс параллелизма представлен на нижнем уровне, при реализации различных операций, используемых в алгоритме, таких как QR разложение методом вращений и матричное умножение.

Параллельная сложность базового QR-алгоритма
Посчитаем параллельную сложность каждой операции по отдельности, а затем, по полученным данным, и всего алгоритма целиком:

1. Параллельная сложность QR разложения методом вращений составляет [math] 11n - 16 [/math]
   
2. Параллельная сложность матричного умножения составляет [math] n [/math].
   
3. Параллельная сложность проверки, является ли матрица диагональной, равна 1.
   

Таким образом параллельная сложность каждой итерации составляет [math] 12n - 15 [/math]. При N итерациях, которые необходимо производить последовательно, итоговая сложность базового алгоритма составляет [math] N * (12n - 15) [/math]

// TODO for Heisenberg

1.9 Входные и выходные данные алгоритма

• Входные данные: плотная квадратная матрица [math]A[/math] (элементы [math]a_{ij}[/math]).

• Объём входных данных: [math]n^2[/math].

• Выходные данные: [math]n[/math] вещественных собственных чисел [math] | l_{i} | [/math] матрицы [math]A[/math]

• Объём выходных данных: [math]n[/math].

1.10 Свойства алгоритма

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

Вычислительная мощность, равная отношению числа операций [math] N * O(n^3) [/math] к суммарному объему входных и выходных данных [math] n^2 + n [/math], для каждой итерации линейна, а для всего алгоритма в целом равна [math] N*n [/math], что позволяет сделать вывод о том, что перемещение данных для обработки не играет важной роли в данном алгоритме.

Вычислительная погрешность растет линейно, из-за использования метода вращений для QR разложения.

Сбалансированность (??????)

Данный алгоритм недетерминирован, так как число итераций зависит от значений элементов матрицы и выход из основного цикла происходит по достижению некоторой точности. В зависимости от этого, вычислительный поток имеет различную длину. На каждой отдельной итерации в свою очередь алгоритм полностью детерминирован, так как детерминированы его обе базовые макрооперации - QR разложение и матричное умножение.

Степнь исхода вершины информационного графа (??????)

"Длинные" дуги в информационном графе (??????)

2 Программная реализация

2.1 Особенности реализации последовательного алгоритма

Основным типом данных для QR алгоритма является матрица. В данной статье представлена реализация на языке C++, что накладывает определенные ограничения на формат хранени матриц: все матрицы хранятся в одномерном массиве в построчном (ROW_MAJOR) формате.

В простейшем случае QR-алгоритм выглядит следующим образом:

int QR_algorithm(float *A, int size)
{
    float *Q = new double[size * size];
    float *R = new double[size * size];
    while(!check_if_upper_triangular(A, size)))
    {
        sgeqrf (ROW_MAJOR, size, size, A, size, tau);
        set_Q_and_R(A, Q, R, size);
        sorgqr(ROW_MAJOR, size, size, size, Q, size, tau);
        sgemm (RowMajor, NoTrans, NoTrans, size, size, size, 1, R, size, Q, size, 0, A, size);
    }
}

В последовательном случае функции dgeqrf, dorgqr и dgemm представляют собой последовательные функции стандарта BLAS, вычисляющие QR разложение(geqrf), матрицу Q после разложения (dorgqr) и матричное умножение (dgemm). Так же здесь используются простейшие функции set_Q_and_R и check_if_upper_triangular. Функция set_Q_and_R получает после QR разложения из матрицы A данные о матрицах Q и R. Функция check_if_upper_triangular проверяет, приведена ли матрица A к треугольном виду, и если да, происходит выход из цикла.

Для работы программы в памяти необходимо хранить 3 матрицы, таким образом общий объем составляет 3 * n^2. Так же необходим вспомогательный массив tau для QR разложения, однако его размер составляет всего n. Так же важно заметить, что различные реализации функций dgeqrf, dorgqr и dgemm обычно используют внутренние вспомогательные массивы для блочной обработки, однако, их размер не превышает n * C, где C - некоторая архитектурно зависимая константа.

Все вычисления рекомендуется производить в одинарной точности.

Как уже говорилось ранее, внешний цикл представленной программы сугубо последователен (итеративен), и основной ресурс параллелизма расположен в вызываемых функциях (geqrf, gemm). Их реализации бывают крайне различны, поэтому данный алгоритм можно легко преобразовать для любой архитектуры, используя набор стандартных библиотечных функций. К примеру, для реализации алгоритма на GPU можно использовать библиотеку MAGMA, на многоядерных CPU - Intel MKL.

// TODO code with heisinberg

2.2 Локальность данных и вычислений

// не нужно???

2.3 Возможные способы и особенности параллельной реализации алгоритма

Как уже говорилось в разделе 2.1., при использовании различных библиотек можно реализовать данный алгоритм практически на всех архитектурах. В данной статье приводятся две реализации: и использованиям библиотеки MKL для многомерных узлов, и библиотеки Magma для вычислений на GPU.

Для распарарллеливания на многоядерных CPU используются функции LAPACKE_sgeqrf, LAPACKE_sorgqr и cblas_sgemm. Распараллеливание оставшихся функций может быть произведено с использованием openMP, так как внутренние циклы в них независимы по данным.

// Добавить про magma (!!!)