Алговики - Вклад [ru]

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

2021-12-07T13:06:36Z

Bulatral: /* Выводы */

{{algorithm
| name = QR-разложения квадратной матрицы методом вращений Гивенса
| serial_complexity = <math>2n^3</math>
| pf_height = <math>11n-16</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>n^2</math>
| output_data = <math>n^2</math>
}}
Автор: [[Участник:Bulatral|Б.И.Валиахметов]]

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

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

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

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===

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

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

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

Для удобства будем полагать, что всевозможные размеры и параметры таковы, что все распределения производятся нацело (например, степени 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).
Этим дополнительным разделением мы задействуем дополнительный ресурс параллелизма и, что не менее важно, эффективно используем кэш-память высокого уровня.

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

Реализация алгоритма выполнена на языке C++11 с использованием стандартных библиотек языка, а также пакетов OpenMPI и OpenBLAS.
Визуализация для результатов исследования параллельных свойств выполнена на языке Python3 с использованием библиотек numpy и matplotlib в среде Jupyter Notebook.
Полный код программы находится в репозитории по [https://github.com/bulatral42/mpiQR ссылке], здесь приведем реализацию основной процедуры алгоритма.

[[Файл:Гивенс t(p=1,2).png|thumb|right|500px|Рисунок 1. Разница во времени между 1 и 2 процессами.]]

[[Файл:Гивенс-3Д,Т=1.png|540px|обрамить|справа|Рисунок 2. График зависимости времени от числа процессов и размера матрицы для T=1 нити]]
[[Файл:Гивенс-3Д,Т=2.png|540px|обрамить|справа|Рисунок 3. График зависимости времени от числа процессов и размера матрицы для T=2 нитей]]
[[Файл:Гивенс-3Д,Т=3.png|540px|обрамить|справа|Рисунок 4. График зависимости времени от числа процессов и размера матрицы для T=3 нитей]]

<source lang="c++">

#include <mpi.h>
#include <omp.h>
#include <iostream>
#include <cmath>
#include <random>
#include <chrono>
#include <unistd.h>
#include "BLAS.h"
#include "LAPACK.h"

void block_qrdec(int N, int block_size, int subblock_size, MPI_Comm comm = MPI_COMM_WORLD)
{
double t_start{0}, t_end{0}, mpi_qr_time{0};
int n_proc, proc_id, n_thr;
MPI_Comm_rank(comm, &proc_id);
MPI_Comm_size(comm, &n_proc);
#pragma omp parallel
{
n_thr = omp_get_num_threads();
}
int n_block_cols = N / (n_proc * block_size), n_block_rows = N / block_size; // per proc

if (proc_id == 0) {
std::cout << "Start mpi-omp-QR, n_proc = " << n_proc << ", n_thr = " << n_thr << std::endl;
std::cout << "N = " << N << ", k = " << block_size << ", m = " << subblock_size;
std::cout << ", nbc = " << n_block_cols << ", nbr = " << n_block_rows << std::endl;
}
// blocks are stored rowwise, every block inside is stored columnwise
double *blocks = new double[n_block_rows * n_block_cols * block_size * block_size];

std::srand(proc_id);
for (int i = 0; i < n_block_rows * n_block_cols * block_size * block_size; ++i) {
//blocks[i] = i + proc_id * n_block_rows * n_block_cols * block_size * block_size;
blocks[i] = -1.0 + 2 * static_cast<double>(std::rand()) / RAND_MAX;
}

// Main algorithm start
t_start = MPI_Wtime();

int proc_rots_cnt = (n_block_rows - proc_id) * n_block_cols - n_proc * ((n_block_cols - 1) * n_block_cols) / 2;
proc_rots_cnt *= block_size * block_size;
proc_rots_cnt -= n_block_cols * (block_size * (block_size + 1)) / 2;

double *my_rotations = new double[2 * proc_rots_cnt];
double *out_rotations = new double[2 * block_size * block_size];
double *cur_rots = my_rotations;

// Main cycle
for (int b_col = 0; b_col < n_block_cols * n_proc; ++b_col) {
int sender_id = b_col % n_proc;
int local_b_col = b_col / n_proc;
double *main_block = blocks + (block_size * block_size) * (b_col * n_block_cols + local_b_col);
for (int b_row = b_col; b_row < n_block_rows; ++b_row) {
double *cur_block = blocks + (block_size * block_size) * (b_row * n_block_cols + local_b_col);
int n_rots{0};
bool triang_zero = b_row == b_col; // Zerorize triangular part or full block
if (proc_id == sender_id) { // Sender part
for (int j = 0; j < block_size; ++j) {
for (int i = triang_zero ? j + 1 : 0; i < block_size; ++i) {
double r{0.0};
dlartg(main_block[j * block_size + j], cur_block[j * block_size + i],
cur_rots[2 * n_rots], cur_rots[2 * n_rots + 1], r);
drot(block_size - j - 1, main_block + (j + 1) * block_size + j, block_size,
cur_block + (j + 1) * block_size + i, block_size,
cur_rots[2 * n_rots], cur_rots[2 * n_rots + 1]);
main_block[j * block_size + j] = r;
cur_block[j * block_size + i] = 0.0;
++n_rots;
}
}
} else {
if (triang_zero) {
n_rots = ((block_size - 1) * block_size) / 2;
} else {
n_rots = block_size * block_size;
}
}
// Send/receive rotations
MPI_Request req_bcast;
double *senrecv_rots = proc_id == sender_id ? cur_rots : out_rotations;
MPI_Ibcast(senrecv_rots, 2 * n_rots, MPI_DOUBLE, sender_id, comm, &req_bcast);

// Perform rotations to other blocks in row
double *up_other_blocks = main_block;
double *dn_other_blocks = cur_block;
int n_subblocks = (block_size / subblock_size) * (n_block_cols - b_col);
if (proc_id == sender_id) { // Current block has already been rotated
up_other_blocks += block_size * block_size;
dn_other_blocks += block_size * block_size;
n_subblocks -= block_size / subblock_size;
} else {
MPI_Status status;
MPI_Wait(&req_bcast, &status);
}
#pragma omp parallel for schedule(static)
for (int subbl_id = 0; subbl_id < n_subblocks; ++subbl_id) {
double *up_subblock = up_other_blocks + block_size * subblock_size * subbl_id;
double *dn_subblock = dn_other_blocks + block_size * subblock_size * subbl_id;
int up_off{0};
int dn_off = triang_zero ? up_off + 1 : 0;
for (int i = 0; i < n_rots; ++i) {
drot(subblock_size, up_subblock + up_off, block_size, dn_subblock + dn_off,
block_size, senrecv_rots[2 * i], senrecv_rots[2 * i + 1]);
++dn_off;
if (dn_off % block_size == 0) {
++up_off;
dn_off = triang_zero ? up_off + 1 : 0;
}
}
}
if (proc_id == sender_id) { // Current block has already been rotated
MPI_Status status;
MPI_Wait(&req_bcast, &status);
cur_rots += 2 * n_rots;
}
}
}
MPI_Barrier(comm);
t_end = MPI_Wtime();
mpi_qr_time = t_end - t_start;
if (proc_id == 0) {
std::cout << "MPI QR time = " << mpi_qr_time << " sec" << std::endl;
}

delete[] my_rotations;
delete[] out_rotations;
delete[] blocks;
}

</source>

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

[[Файл:Гивенс-N=4,T=1.png|обрамить|справа|Рисунок 5. График зависимости времени от числа процессов, <math>N=4096, T=1</math>]]
[[Файл:Гивенс-N=5,T=1.png|обрамить|справа|Рисунок 6. График зависимости времени от числа процессов, <math>N=8192, T=1</math>]]
[[Файл:Гивенс-N=8192,T=1,2,3.png|обрамить|справа|Рисунок 7. График зависимости времени от количества процессов при разном числе нитей, <math>N=8192</math>]]

Эксперименты по установлению параллельных свойств алгоритма проводились на кластере ИВМ РАН.
Характеристики узлов кластера на момент написания раздела (полное описание по [https://cluster2.inm.ras.ru/info/nodes/ ссылке]):
* 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).

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

1. На предварительных запусках были установлены оптимальные параметры блочного разбиения: <math>k = 128</math>, <math>m = 16</math>.

2. Размеры матриц выбирались в виде степеней 2: <math>N = 2^s, s \in \overline{8, 16}</math>.

3. Число процессов выбиралось в виде степеней 2:<math>P = 2^l, l \in \overline{0, 6}</math>, при условии <math>P \geq \dfrac{N}{k}</math>.

4. Число нитей выбиралось из набора: <math>T \in \{1, 2, 4\}</math>.

5. Для каждой конфигурации параметров проводилось 3 запуска, измеренное время работы усреднялось.

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

На Рисунках 2, 3, 4 приведены графики зависимости времени от размеров матрицы и количества процессов (для удобства рассмотрения все величины прологарифмированы)

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

На Рисунке 7 приведены графики зависимости времени работы программы от количества процессов при фиксированном размере матрицы <math>N = 8192</math> и различном количестве потоков <math>T \in \{1, 2, 4\}</math>.
Можем заметить, что использование двух потоков дает практически двукратный выигрыш во времени по сравнению с одним, однако 4 потока работают уже дольше.
Вероятно, это связано с накладными расходами на их синхронизацию и менее равномерным распределением задач.

==== Выводы ====

В рамках практической работы мы реализовали параллельную версию алгоритма QR-разложения квадратной вещественной матрицы методом вращений Гивенса и изучили масштабируемость написанной программы.
Из результатов проделанной работы можем сделать следующие выводы:
* Реализованная нами программа выполняет поставленную задачу.
* Параллельная программа по сравнению с последовательной версией ускоряется в большее число раз, чем число процессов - это говорит о хороших параллельный свойствах алгоритма.
* На малом (от 2 до 8) количестве процессов имеет место почти линейное ускорение программы.
* При большом количестве процессов (<math>\geq 16</math>) или потоков (4) на исследованных размерах матриц программа заметно не ускоряется, так как накладные расходы и несбалансированность "перевешивают" ускорение.

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===

=== Существующие реализации алгоритма ===

== Литература ==
<references />

[[En: Givens method]]

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

2021-12-07T13:05:13Z

Bulatral: Основная часть раздела 2.4

{{algorithm
| name = QR-разложения квадратной матрицы методом вращений Гивенса
| serial_complexity = <math>2n^3</math>
| pf_height = <math>11n-16</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>n^2</math>
| output_data = <math>n^2</math>
}}
Автор: [[Участник:Bulatral|Б.И.Валиахметов]]

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

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

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

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===

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

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

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

Для удобства будем полагать, что всевозможные размеры и параметры таковы, что все распределения производятся нацело (например, степени 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).
Этим дополнительным разделением мы задействуем дополнительный ресурс параллелизма и, что не менее важно, эффективно используем кэш-память высокого уровня.

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

Реализация алгоритма выполнена на языке C++11 с использованием стандартных библиотек языка, а также пакетов OpenMPI и OpenBLAS.
Визуализация для результатов исследования параллельных свойств выполнена на языке Python3 с использованием библиотек numpy и matplotlib в среде Jupyter Notebook.
Полный код программы находится в репозитории по [https://github.com/bulatral42/mpiQR ссылке], здесь приведем реализацию основной процедуры алгоритма.

[[Файл:Гивенс t(p=1,2).png|thumb|right|500px|Рисунок 1. Разница во времени между 1 и 2 процессами.]]

[[Файл:Гивенс-3Д,Т=1.png|540px|обрамить|справа|Рисунок 2. График зависимости времени от числа процессов и размера матрицы для T=1 нити]]
[[Файл:Гивенс-3Д,Т=2.png|540px|обрамить|справа|Рисунок 3. График зависимости времени от числа процессов и размера матрицы для T=2 нитей]]
[[Файл:Гивенс-3Д,Т=3.png|540px|обрамить|справа|Рисунок 4. График зависимости времени от числа процессов и размера матрицы для T=3 нитей]]

<source lang="c++">

#include <mpi.h>
#include <omp.h>
#include <iostream>
#include <cmath>
#include <random>
#include <chrono>
#include <unistd.h>
#include "BLAS.h"
#include "LAPACK.h"

void block_qrdec(int N, int block_size, int subblock_size, MPI_Comm comm = MPI_COMM_WORLD)
{
double t_start{0}, t_end{0}, mpi_qr_time{0};
int n_proc, proc_id, n_thr;
MPI_Comm_rank(comm, &proc_id);
MPI_Comm_size(comm, &n_proc);
#pragma omp parallel
{
n_thr = omp_get_num_threads();
}
int n_block_cols = N / (n_proc * block_size), n_block_rows = N / block_size; // per proc

if (proc_id == 0) {
std::cout << "Start mpi-omp-QR, n_proc = " << n_proc << ", n_thr = " << n_thr << std::endl;
std::cout << "N = " << N << ", k = " << block_size << ", m = " << subblock_size;
std::cout << ", nbc = " << n_block_cols << ", nbr = " << n_block_rows << std::endl;
}
// blocks are stored rowwise, every block inside is stored columnwise
double *blocks = new double[n_block_rows * n_block_cols * block_size * block_size];

std::srand(proc_id);
for (int i = 0; i < n_block_rows * n_block_cols * block_size * block_size; ++i) {
//blocks[i] = i + proc_id * n_block_rows * n_block_cols * block_size * block_size;
blocks[i] = -1.0 + 2 * static_cast<double>(std::rand()) / RAND_MAX;
}

// Main algorithm start
t_start = MPI_Wtime();

int proc_rots_cnt = (n_block_rows - proc_id) * n_block_cols - n_proc * ((n_block_cols - 1) * n_block_cols) / 2;
proc_rots_cnt *= block_size * block_size;
proc_rots_cnt -= n_block_cols * (block_size * (block_size + 1)) / 2;

double *my_rotations = new double[2 * proc_rots_cnt];
double *out_rotations = new double[2 * block_size * block_size];
double *cur_rots = my_rotations;

// Main cycle
for (int b_col = 0; b_col < n_block_cols * n_proc; ++b_col) {
int sender_id = b_col % n_proc;
int local_b_col = b_col / n_proc;
double *main_block = blocks + (block_size * block_size) * (b_col * n_block_cols + local_b_col);
for (int b_row = b_col; b_row < n_block_rows; ++b_row) {
double *cur_block = blocks + (block_size * block_size) * (b_row * n_block_cols + local_b_col);
int n_rots{0};
bool triang_zero = b_row == b_col; // Zerorize triangular part or full block
if (proc_id == sender_id) { // Sender part
for (int j = 0; j < block_size; ++j) {
for (int i = triang_zero ? j + 1 : 0; i < block_size; ++i) {
double r{0.0};
dlartg(main_block[j * block_size + j], cur_block[j * block_size + i],
cur_rots[2 * n_rots], cur_rots[2 * n_rots + 1], r);
drot(block_size - j - 1, main_block + (j + 1) * block_size + j, block_size,
cur_block + (j + 1) * block_size + i, block_size,
cur_rots[2 * n_rots], cur_rots[2 * n_rots + 1]);
main_block[j * block_size + j] = r;
cur_block[j * block_size + i] = 0.0;
++n_rots;
}
}
} else {
if (triang_zero) {
n_rots = ((block_size - 1) * block_size) / 2;
} else {
n_rots = block_size * block_size;
}
}
// Send/receive rotations
MPI_Request req_bcast;
double *senrecv_rots = proc_id == sender_id ? cur_rots : out_rotations;
MPI_Ibcast(senrecv_rots, 2 * n_rots, MPI_DOUBLE, sender_id, comm, &req_bcast);

// Perform rotations to other blocks in row
double *up_other_blocks = main_block;
double *dn_other_blocks = cur_block;
int n_subblocks = (block_size / subblock_size) * (n_block_cols - b_col);
if (proc_id == sender_id) { // Current block has already been rotated
up_other_blocks += block_size * block_size;
dn_other_blocks += block_size * block_size;
n_subblocks -= block_size / subblock_size;
} else {
MPI_Status status;
MPI_Wait(&req_bcast, &status);
}
#pragma omp parallel for schedule(static)
for (int subbl_id = 0; subbl_id < n_subblocks; ++subbl_id) {
double *up_subblock = up_other_blocks + block_size * subblock_size * subbl_id;
double *dn_subblock = dn_other_blocks + block_size * subblock_size * subbl_id;
int up_off{0};
int dn_off = triang_zero ? up_off + 1 : 0;
for (int i = 0; i < n_rots; ++i) {
drot(subblock_size, up_subblock + up_off, block_size, dn_subblock + dn_off,
block_size, senrecv_rots[2 * i], senrecv_rots[2 * i + 1]);
++dn_off;
if (dn_off % block_size == 0) {
++up_off;
dn_off = triang_zero ? up_off + 1 : 0;
}
}
}
if (proc_id == sender_id) { // Current block has already been rotated
MPI_Status status;
MPI_Wait(&req_bcast, &status);
cur_rots += 2 * n_rots;
}
}
}
MPI_Barrier(comm);
t_end = MPI_Wtime();
mpi_qr_time = t_end - t_start;
if (proc_id == 0) {
std::cout << "MPI QR time = " << mpi_qr_time << " sec" << std::endl;
}

delete[] my_rotations;
delete[] out_rotations;
delete[] blocks;
}

</source>

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

[[Файл:Гивенс-N=4,T=1.png|обрамить|справа|Рисунок 5. График зависимости времени от числа процессов, <math>N=4096, T=1</math>]]
[[Файл:Гивенс-N=5,T=1.png|обрамить|справа|Рисунок 6. График зависимости времени от числа процессов, <math>N=8192, T=1</math>]]
[[Файл:Гивенс-N=8192,T=1,2,3.png|обрамить|справа|Рисунок 7. График зависимости времени от количества процессов при разном числе нитей, <math>N=8192</math>]]

Эксперименты по установлению параллельных свойств алгоритма проводились на кластере ИВМ РАН.
Характеристики узлов кластера на момент написания раздела (полное описание по [https://cluster2.inm.ras.ru/info/nodes/ ссылке]):
* 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).

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

1. На предварительных запусках были установлены оптимальные параметры блочного разбиения: <math>k = 128</math>, <math>m = 16</math>.

2. Размеры матриц выбирались в виде степеней 2: <math>N = 2^s, s \in \overline{8, 16}</math>.

3. Число процессов выбиралось в виде степеней 2:<math>P = 2^l, l \in \overline{0, 6}</math>, при условии <math>P \geq \dfrac{N}{k}</math>.

4. Число нитей выбиралось из набора: <math>T \in \{1, 2, 4\}</math>.

5. Для каждой конфигурации параметров проводилось 3 запуска, измеренное время работы усреднялось.

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

На Рисунках 2, 3, 4 приведены графики зависимости времени от размеров матрицы и количества процессов (для удобства рассмотрения все величины прологарифмированы)

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

На Рисунке 7 приведены графики зависимости времени работы программы от количества процессов при фиксированном размере матрицы <math>N = 8192</math> и различном количестве потоков <math>T \in \{1, 2, 4\}</math>.
Можем заметить, что использование двух потоков дает практически двукратный выигрыш во времени по сравнению с одним, однако 4 потока работают уже дольше.
Вероятно, это связано с накладными расходами на их синхронизацию и менее равномерным распределением задач.

==== Выводы ====

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

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===

=== Существующие реализации алгоритма ===

== Литература ==
<references />

[[En: Givens method]]

Файл:Гивенс-N=8192,T=1,2,3.png

2021-12-07T12:55:34Z

Bulatral:

График зависимости времени от количества процессов при разном числе нитей

Файл:Гивенс-N=5,T=1.png

2021-12-07T12:45:09Z

Bulatral:

График зависимости времени от размера матрицы

Файл:Гивенс-N=4,T=1.png

2021-12-07T12:43:48Z

Bulatral:

График зависимости времени от размера матрицы

Файл:Гивенс-3Д,T=1,2,3.png

2021-12-07T12:21:31Z

Bulatral:

График зависимости времени от числа процессов и размера матрицы для T=1,2,3 нитей

Файл:Гивенс-3Д,Т=3.png

2021-12-07T12:11:43Z

Bulatral:

График зависимости времени от числа процессов и размера матрицы для T=3 нитей

Файл:Гивенс-3Д,Т=2.png

2021-12-07T12:09:25Z

Bulatral:

График зависимости времени от числа процессов и размера матрицы для 2 нитей

Файл:Гивенс-3Д,Т=1.png

2021-12-07T12:02:38Z

Bulatral:

3Д график зависимости времени от числа процессов и размера матрицы для 1 нити

Файл:Гивенс t(p=1,2).png

2021-12-07T11:23:47Z

Bulatral:

зависимость времени от размера матрицы при 1 и 2 процессах

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

2021-12-07T01:14:40Z

Bulatral: Новая страница: «{{algorithm | name = QR-разложения квадратной матрицы методом вращений Гивенса | serial_complexity...»

{{algorithm
| name = QR-разложения квадратной матрицы методом вращений Гивенса
| serial_complexity = <math>2n^3</math>
| pf_height = <math>11n-16</math>
| pf_width = <math>O(n^2)</math>
| input_data = <math>n^2</math>
| output_data = <math>n^2</math>
}}
Автор: [[Участник:Bulatral|Б.И.Валиахметов]]

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

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

Математическое описание алгоритма и его последовательная версия приведены на портале AlgoWiki: [https://algowiki-project.org/ru/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%93%D0%B8%D0%B2%D0%B5%D0%BD%D1%81%D0%B0_(%D0%B2%D1%80%D0%B0%D1%89%D0%B5%D0%BD%D0%B8%D0%B9)_QR-%D1%80%D0%B0%D0%B7%D0%BB%D0%BE%D0%B6%D0%B5%D0%BD%D0%B8%D1%8F_%D0%BA%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82%D0%BD%D0%BE%D0%B9_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D1%8B_(%D0%B2%D0%B5%D1%89%D0%B5%D1%81%D1%82%D0%B2%D0%B5%D0%BD%D0%BD%D1%8B%D0%B9_%D1%82%D0%BE%D1%87%D0%B5%D1%87%D0%BD%D1%8B%D0%B9_%D0%B2%D0%B0%D1%80%D0%B8%D0%B0%D0%BD%D1%82) Метод Гивенса (вращений) QR-разложения квадратной матрицы (вещественный точечный вариант) ]

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

=== Особенности реализации последовательного алгоритма ===
=== Локальность данных и вычислений ===
=== Возможные способы и особенности параллельной реализации алгоритма ===
=== Масштабируемость алгоритма и его реализации ===

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

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

==== Масштабируемость реализации алгоритма ====
Для удобства будем полагать, что всевозможные размеры и параметры таковы, что все распределения производятся нацело (например, степени 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).
Этим дополнительным разделением мы задействуем дополнительный ресурс параллелизма и, что не менее важно, эффективно используем кэш-память высокого уровня.

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

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

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

Эксперименты по установлению параллельных свойств алгоритма проводились на кластере ИВМ РАН.
Характеристики узлов кластера на момент написания раздела (полное описание по [https://cluster2.inm.ras.ru/info/nodes/ ссылке]):
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}

=== Динамические характеристики и эффективность реализации алгоритма ===
=== Выводы для классов архитектур ===

=== Существующие реализации алгоритма ===

== Литература ==
<references />

[[En: Givens method]]