Участник: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. Полный код программы находится в репозитории по ссылке, здесь приведем реализацию основной процедуры алгоритма.

Рисунок 1. Разница во времени между 1 и 2 процессами.

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

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

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

#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;
}

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

Рисунок 5. График зависимости времени от числа процессов, [math]N=4096, T=1[/math]

Рисунок 6. График зависимости времени от числа процессов, [math]N=8192, T=1[/math]

Рисунок 7. График зависимости времени от количества процессов при разном числе нитей, [math]N=8192[/math]

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

• 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 потока работают уже дольше. Вероятно, это связано с накладными расходами на их синхронизацию и менее равномерным распределением задач.

2.4.5 Выводы

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

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