Участник:Bulatral/QR-разложение плотной вещественной матрицы методом вращений Гивенса

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

Автор: Б.И.Валиахметов

Необходимо построить QR-разложение квадратной вещественной [math]n \times n[/math] матрицы [math]A[/math] методом вращений Гивенса.

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

Математическое описание алгоритма и его последовательная версия приведены на портале AlgoWiki: Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный точечный вариант)

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

2.4 Масштабируемость алгоритма и его реализации

2.4.1 Масштабируемость алгоритма

Отметим, что по сути, после применения вращений к исходной матрице мы получаем матрицу [math]R[/math] и для получения матрицы [math]Q[/math] необходимо выполнить обратные вращения над единичной матрицей. Однако на практике это не нужно, так как в большинстве случаев требуется умножать на [math]Q[/math] или [math]Q^*[/math], а для этого достаточно применить подсчитанные ранее вращения. Также этот подход дает меньшую численную ошибку и выполняется быстрее.

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

2.4.2 Масштабируемость реализации алгоритма

Для удобства будем полагать, что всевозможные размеры и параметры таковы, что все распределения производятся нацело (например, степени 2).

Зададим размер блока [math]k[/math] и размер подблока [math]m[/math]. Пусть нам доступны [math]p[/math] процессов по [math]t[/math] нити в каждом. Разобьем матрицу [math]A[/math] на блоки размера [math]k[/math] и распределим блочные столбцы по процессам так: столбец с номером [math]j[/math] хранится на процессе с номером [math]j \mod p[/math], чтобы распределение количества вращений было более равномерным и узлы меньше простаивали в ожидании пересылок. Также это способствует активному использованию кэш-памяти. Такое распределение легко произвести, если матрица задается как функция вычисления ее элемента, иначе необходимо провести некоторые пересылки (их мы не учитываем при подсчете времени работы алгоритма). На каждом процессе блоки по блочным строкам (а сами блоки внутри - по стодбцам), так они задействуются совместно на шагах алгоритма.

Далее будем по очереди обнулять поддиагонали диагональных блоков и блоки ниже диагональных. Процесс, на котором хранится блок (ведущий), сначала вычисляет необходимые вращения для данного блока (процедура DLARTG библиотеки LAPACK), а затем рассылает их по остальным процессам (процедурой MPI\_IBcast стандарта MPI). Во время неблокирующей пересылки ведущий процесс производит вращения на оставшиеся блоки в блочной строке и приступает к вычислению вращений для последующих блоков в столбце.

Как проводятся вращения на блоке, который сейчас не обнуляется? Каждый процесс (в том числе и ведущий), получив набор вращений, применяет их в той же блочной строке, что и процесс-отправитель. Каждый блок дополнительно разбивается на блочные подстолбцы ширины [math]m[/math] и данные подблоки статически распределяются по нитям на процессе (директива omp parallel for библиотеки OpenMP), которые и производят вращения (процедура DROT библиотеки BLAS). Этим дополнительным разделением мы задействуем дополнительный ресурс параллелизма и, что не менее важно, эффективно используем кэш-память высокого уровня.

2.4.3 Реализация алгоритма

Реализация алгоритма выполнена на языке C++11 с использованием стандартных библиотек языка, а также пакетов OpenMPI и OpenBLAS. Визуализация для результатов исследования параллельных свойств выполнена на языке Python3 с использованием библиотек numpy и matplotlib в среде Jupyter Notebook. Код программы находится в репозитории по ссылке.

2.4.4 Численные эксперименты

Эксперименты по установлению параллельных свойств алгоритма проводились на кластере ИВМ РАН. Характеристики узлов кластера на момент написания раздела (полное описание по ссылке):

   Compute Node Arbyte Alkazar R2Q50 G5.
   40 ядер (два 20-ядерных процессора Intel Xeon Gold 6230@2.10ГГц).
   Оперативная память: 256/384 Гб.
   Дисковая память: 480 Гб SSD.
   Операционная система: SUSE Linux Enterprise Server 15 SP2 (x86\_64).

Были выбраны следующие параметры запусков.

   На предварительных запусках были установлены оптимальные параметры блочного разбиения: $k = 128$, $m = 16$.
   Размеры матриц выбирались в виде степеней 2: $N = 2^s, s \in \overline{8, 16}$.
   Число процессов выбиралось в виде степеней 2: $P = 2^l, l \in \overline{0, 6}$, при условии $P \geq \dfrac{N}{k}$.
   Число нитей выбиралось из набора: $T \in \{1, 2, 4\}$.
   Для каждой конфигурации параметров проводилось 3 запуска, измеренное время работы усреднялось.

В таблицах \ref{tab:1thr}, \ref{tab:2thr}, \ref{tab:3thr} приведены результаты времени работы программы в зависимости от параметров. Некоторые запуски не проводились ввиду их слишком большой потенциальной длительности или слишком малого необходимого значения размера блока.

\begin{table}[h]

   \centering
   \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
   \hline
       P \textbackslash N & 256 & 512 & 1024 & 2048 & 4096 & 8192 & 16384 & 32768 & 65536  \\ \hline
       1 & 0.0402 & 0.3432 & 2.824 & 29.84 & 498.4 & N/A & N/A & N/A & N/A \\ \hline
       2 & N/A & 0.0459 & 0.2231 & 1.372 & 9.505 & 99.99 & 807.177 & N/A & N/A \\ \hline
       4 & N/A & 0.03902 & 0.1529 & 0.7388 & 4.615 & 34.43 & 250.7 & 1555.2 & N/A \\ \hline
       8 & N/A & N/A & 0.3030 & 0.6821 & 3.046 & 15.33 & 89.11 & 588.2 & N/A \\ \hline
       16 & N/A & N/A & N/A & 0.6065 & 2.502 & 10.68 & 49.22 & 254.7 & N/A \\ \hline
       32 & N/A & N/A & N/A & N/A & 2.445 & 9.951 & 41.64 & 181.14 & N/A \\ \hline
       64 & N/A & N/A & N/A & N/A & N/A & 10.24 & 39.56 & 161.6 & 675.9 \\ \hline
   \end{tabular}
   \caption{Время работы в сек., $T = 1$ нить}
   \label{tab:1thr}

\end{table} \begin{table}[h]

   \centering
   \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
   \hline
       P \textbackslash N & 256 & 512 & 1024 & 2048 & 4096 & 8192 & 16384 & 32768 & 65536  \\ \hline
       2 & N/A & 0.08189 & 0.2260 & 1.274 & 8.117 & 56.26 & 416.2 & N/A & N/A \\ \hline
       4 & N/A & 0.06885 & 0.1688 & 0.7482 & 3.926 & 22.35 & 146.8 & 990.4 & N/A \\ \hline
       8 & N/A & N/A & 0.1630 & 0.6467 & 2.744 & 13.14 & 67.82 & 374.315 & N/A \\ \hline
       16 & N/A & N/A & N/A & 0.9567 & 3.663 & 14.84 & 64.07 & 294.9 & N/A \\ \hline
       32 & N/A & N/A & N/A & N/A & 3.785 & 14.70 & 57.52 & 239.5 & N/A \\ \hline
       64 & N/A & N/A & N/A & N/A & N/A & 14.54 & 55.85 & 223.2 & 911.6 \\ \hline
   \end{tabular}
   \caption{Время работы в сек., $T = 2$ нити}
   \label{tab:2thr}

\end{table} \begin{table}[h]

   \centering
   \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
   \hline
       P \textbackslash N & 256 & 512 & 1024 & 2048 & 4096 & 8192 & 16384 & 32768 & 65536  \\ \hline
       2 & N/A & 0.1687 & 0.6239 & 2.654 & 13.17 & 76.14 & 487.3 & N/A & N/A \\ \hline
       4 & N/A & 0.1299 & 0.4685 & 1.905 & 8.141 & 40.02 & 214.4 & 1154.6 & N/A \\ \hline
       8 & N/A & N/A & 0.2759 & 0.9752 & 3.870 & 16.59 & 75.56 & 374.2 & N/A \\ \hline
       16 & N/A & N/A & N/A & 1.058 & 3.783 & 15.13 & 62.32 & 269.9 & N/A \\ \hline
       32 & N/A & N/A & N/A & N/A & 3.915 & 14.96 & 59.65 & 242.9 & N/A \\ \hline
       64 & N/A & N/A & N/A & N/A & N/A & 15.28 & 59.28 & 236.4 & 955.5 \\ \hline
   \end{tabular}
   \caption{Время работы в сек., $T = 3$ нитей}
   \label{tab:3thr}

\end{table}

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

Ниже на Рисунках \ref{graph:3d_1}, \ref{graph:3d_2}, \ref{graph:3d_4} приведены зависимости из таблиц выше.

\begin{figure}[h] \centering \includegraphics[scale=0.45]{pics/pic_T1.png} \caption{$T = 1$} \label{graph:3d_1} \end{figure} \begin{figure}[h] \centering \includegraphics[scale=0.45]{pics/pic_T2.png} \caption{$T = 1$} \label{graph:3d_2} \end{figure} \begin{figure}[h] \centering \includegraphics[scale=0.45]{pics/pic_T3.png} \caption{$T = 1$} \label{graph:3d_4} \end{figure}

На Рисунке \ref{graph:2d_P} приведены графики зависимости времени работы программы от количества процессов при фиксированном размере матрицы $N \in \{4096, 81921\}$ и количество потоков $T = 1$. Видим, что на малом числе процессов ускорение почти линейное, затем достигается некоторая граница и время начинает несколько расти из-за увеличения накладных расходов. Также это может быть связано с тем, что процессы располагаются на большем числе отдельных узлов и скорость передачи сообщений между ними меньше, чем внутри одного узла.

\begin{figure}[h] \centering \begin{subfigure}[b]{0.45\textwidth} \includegraphics[scale=0.55]{pics/graph_N4.png} \caption{$N = 4096$} \end{subfigure} \hfill \begin{subfigure}[b]{0.45\textwidth} \includegraphics[scale=0.55]{pics/graph_N5.png} \caption{$N = 8192$} \end{subfigure} \caption{$time(P)$} \label{graph:2d_P} \end{figure}

На Рисунке \ref{graph:2d_T} приведены графики зависимости времени работы программы от количества процессов при фиксированном размере матрицы $N = 81921$ и различном количестве потоков $T \in \{1, 2, 4\}$. Можем заметить, что использование двух потоков дает практически двукратный выигрыш во времени по сравнению с одним, однако 4 потока работают уже дольше. Вероятно, это связано с накладными расходами на их синхронизацию и менее равномерным распределением задач.

\begin{figure}[h] \centering \includegraphics[scale=0.75]{pics/graph_T123.png} \caption{$N = 8192$} \label{graph:2d_T} \end{figure}

\section{Выводы}

В рамках практической работы мы реализовали параллельную версию алгоритма QR-разложения квадратной вещественной матрицы методом вращений Гивенса и изучили масштабируемость написанной программы. Из результатов проделанной работы можем сделать следующие выводы: \begin{enumerate}

   \item Реализованная нами программа выполняет поставленную задачу.
   \item Параллельная программа по сравнению с последовательной версией ускоряется в большее число раз, чем число процессов - это говорит о хороших параллельный свойствах алгоритма.
   \item На малом (от 2 до 8) количестве процессов имеет место почти линейное ускорение программы.
   \item При большом количестве процессов ($\geq 16$) или потоков (4) на исследованных размерах матриц программа заметно не ускоряется, так как накладные расходы и несбалансированность "перевешивают" ускорение.

\end{enumerate}