Сеть обменной сортировки со слиянием Бэтчера

Здесь и далее будем приводить примеры кода на Си со следующими соглашениями.

1. Массивы типа T обозначаются как array(T).
2. Считаем, что все функции работы с массивами определены в каком-нибудь array.h и делают ровно то, что написано в названии.
3. Кроме функции array_push(): для векторов она добавляет переданное значение в конец массива, при необходимости выделяя память; для обычных массивов — записывает переданное значение после последнего записанного, начиная с нулевого.

void S(array(int) procs_up, array(int) procs_down)
{
    int proc_count = array_size(procs_up) + array_size(procs_down);
    if (proc_count == 1) {
        return;
    } else if (proc_count == 2) {
        array_push(&comparators, ((pair_t){ procs_up[0], procs_down[0] }));
        return;
    }
    
    array(int) procs_up_odd = array_new(array_size(procs_up) / 2 + array_size(procs_up) % 2, int);
    array(int) procs_down_odd = array_new(array_size(procs_down) / 2 + array_size(procs_down) % 2, int);
    array(int) procs_up_even = array_new(array_size(procs_up) / 2, int);
    array(int) procs_down_even = array_new(array_size(procs_down) / 2, int);
    array(int) procs_result = array_new(array_size(procs_up) + array_size(procs_down), int);
    
    for (int i = 0; i < array_size(procs_up); i++) {
        if (i % 2) {
            array_push(&procs_up_even, procs_up[i]);
        } else {
            array_push(&procs_up_odd, procs_up[i]);
        }
    }
    for (int i = 0; i < array_size(procs_down); i++) {
        if (i % 2) {
            array_push(&procs_down_even, procs_down[i]);
        } else {
            array_push(&procs_down_odd, procs_down[i]);
        }
    }
    
    S(procs_up_odd, procs_down_odd);
    S(procs_up_even, procs_down_even);
    
    array_concatenate(&procs_result, procs_up, procs_down);

    for (int i = 1; i + 1 < array_size(procs_result); i += 2) {
        array_push(&comparators, ((pair_t){ procs_result[i], procs_result[i + 1] }));
    }
    
    array_delete(&procs_up_odd);
    array_delete(&procs_down_odd);
    array_delete(&procs_up_even);
    array_delete(&procs_down_even);
    array_delete(&procs_result);
}

void B(array(int) procs)
{
    if (array_size(procs) == 1) {
        return;
    }
    
    array(int) procs_up = array_new(array_size(procs) / 2, int);
    array(int) procs_down = array_new(array_size(procs) / 2 + array_size(procs) % 2, int);
    
    array_copy(procs_up, procs, 0, array_size(procs_up));
    array_copy(procs_down, procs, array_size(procs_up), array_size(procs_down));
    
    B(procs_up);
    B(procs_down);
    S(procs_up, procs_down);
    
    array_delete(&procs_up);
    array_delete(&procs_down);
}

Код функций с нормальной подсветкой синтаксиса на гитхабе.

Для упрощения будем использовать встроенную функцию сортировки qsort(). Размер сортируемого массива на каждом процессоре определяется как:

int proc_count; // всего процессоров
int elems_count; // всего элементов
int elems_count_new = elems_count + (elems_count % proc_count ? proc_count - elems_count % proc_count : 0);
int elems_per_proc_count = elems_count_new / proc_count;

Тогда первый этап состоит всего из одной строки кода:

qsort(elems_result, elems_per_proc_count, array_item_size(elems_result), compare_uint32);

На гитхабе.

Как уже было сказано, данный этап заключается в последовательной обработке каждого компаратора из массива comparators. Процессор из первого входа компаратора отправляет свой массив на второй вход, другой процессор поступает аналогично, в результате чего каждый процессор имеет два массива: свой (elems_result) и другого процессора (elems_current).

Затем на каждом процессоре формируется третий массив (elems_temp), состоящий из наименьших элементов (для первого процессора) или наибольших (для второго процессора).

После чего массив elems_temp записывается на место массива elems_result.

for (int i = 0; i < array_size(comparators); i++) {
    pair_t comparator = comparators[i];
    if (rank == comparator.a) {
        MPI_Send(elems_result, elems_per_proc_count, MPI_UNSIGNED, comparator.b, 0, MPI_COMM_WORLD);
        MPI_Recv(elems_current, elems_per_proc_count, MPI_UNSIGNED, comparator.b, 0, MPI_COMM_WORLD, &status);
        
        for (int res_index = 0, cur_index = 0, tmp_index = 0; tmp_index < elems_per_proc_count; tmp_index++) {
            uint32_t result = elems_result[res_index];
            uint32_t current = elems_current[cur_index];
            if (result < current) {
                elems_temp[tmp_index] = result;
                res_index++;
            } else {
                elems_temp[tmp_index] = current;
                cur_index++;
            }
        }
        
        swap_ptr(&elems_result, &elems_temp);
    } else if (rank == comparator.b) {
        MPI_Recv(elems_current, elems_per_proc_count, MPI_UNSIGNED, comparator.a, 0, MPI_COMM_WORLD, &status);
        MPI_Send(elems_result, elems_per_proc_count, MPI_UNSIGNED, comparator.a, 0, MPI_COMM_WORLD);
        
        int start = elems_per_proc_count - 1;
        for (int res_index = start, cur_index = start, tmp_index = start; tmp_index >= 0; tmp_index--) {
            uint32_t result = elems_result[res_index];
            uint32_t current = elems_current[cur_index];
            if (result > current) {
                elems_temp[tmp_index] = result;
                res_index--;
            } else {
                elems_temp[tmp_index] = current;
                cur_index--;
            }
        }
        
        swap_ptr(&elems_result, &elems_temp);
    }
}

На гитхабе.

Идея заключается в следующем. Пусть у нас есть массив из 15 элементов, по 3 элемента на процессор (переменная elems_per_proc_count) (рисунок 5 ниже). Необходимо пропустить первые 5 элементов (переменная skip).

Рисунок 5. Пример распределения массива по процессорам. Необходимо вывести обведенные синим элементы.

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

Возможны три случая.

1. Процессор не выводит ничего, т. е. смещение — ноль, количество выводимых — ноль.
2. Процессор выводит некоторые свои элементы, т. е. смещение — число пропускаемых элементов, количество выводимых — разность elems_per_proc_count и смещения.
3. Процессор выводит все свои элементы массива, т. е. смещение — ноль, количество выводимых — elems_per_proc_count.

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

Таким образом, после группировки всех условий получаем выражение для print_offset.

По приведенному ниже коду не составит труда разобраться в формуле для вычисления количества элементов.

// Количество элементов из общего числа, которые не нужно выводить.
int skip = elems_count_new - elems_count;

// Номер элемента, с которого начинается вывод для каждого процессора.
int print_offset = (skip / elems_per_proc_count == rank) * (skip % elems_per_proc_count);

// Количество выводимых процессором элементов.
int print_count = (skip / elems_per_proc_count <= rank) * elems_per_proc_count - print_offset;

// Число элементов выводит только один процессор с рангом 0.
MPI_File_write_ordered(output, &elems_count, rank == 0, MPI_UNSIGNED, &status);

// Последовательно выводим элементы на основе рассчитанных смещений и количеств для каждого процессора.
MPI_File_write_ordered(output, (unsigned char *)elems_result + print_offset * array_item_size(elems_result), print_count, MPI_UNSIGNED, &status);

На гитхабе.