MPI并行编程系列三：并行正则采样排序PSRS

    快速排序算法的效率相对较高，并行算法在理想的情况下时间复杂度可达到o(n)，但并行快速排序算法有一个严重的问题：会造成严重的负载不平衡，最差情况下算法的复杂度可达o(n^2)。本篇我们介绍一种基于均匀划分的负载平衡的并行排序算法------并行正则采样排序（Parallel Sorting by Regular Sampling）。

一、算法的基本思想

      假设待排序的元素n个，处理器p个。

      首先将这n个元素均匀的分成p部分，每部分包含n/p个元素。每个处理器负责其中的一部分，并对其进行局部排序。为确定局部有序序列在整个序列中的位置，每个处理器从各自的局部有序序列中选取几个代表元素，将这些代表元素进行排序后选出p-1个主元。每个处理器根据这p-1个主元将自己的局部有序序列分成p段。然后通过全局交换的方式，将p段有序序列分发给对应的处理器，使第i个处理器都拥有各个处理器的第i段，共p段有序序列。每个处理器对着p段有序序列进行排序。最后，将各个处理器的有序段依次汇合起来，就是全局有序序列了。

二、算法描述

    根据算法的基本思想，我们对算法的描述如下：

输入：n个待排序的序列

输出：分布在各个处理器上，得到全局有序的数据序列

1）无序序列的划分及局部排序

    根据数据快的划分方法（请看系列一），将无序序列划分成p部分，每个处理器对其中的一部分进行串行快速排序，这样每个处理器就会拥有一个局部有序序列。

2）选取代表元素

    每个处理器从局部有序序列中选取第w,2w,...,(p-1)w共p-1个代表元素。其中w = n/p^2。

3）确定主元

    每个处理器都将自己选取好的代表元素发送给处理器p0。p0对这p段有序序列做多路归并排序，再从这排序后的序列中选取第p-1,2(p-1), ...,(p-1)(p-1)共p-1个元素作为主元。

4）分发主元

    p0将这p-1个主元分发给各个处理器。

5）局部有序序列划分

    每个处理器在接收到主元后，根据主元将自己的局部有序序列划分成p段。

6）p段有序序列的分发

    每个处理器将自己的第i段发送给第i个处理器，是处理器i都拥有所有处理器的第i段。

7）多路排序

     每个处理器将上一步得到的p段有序序列做多路归并。

经过这7步后，一次将每个处理器的数据取出，这些数据是有序的。

三、算法分析

   1）负载均衡分析：

     因为这个算法是一个负载平衡的算法，者从第1）步中就可以看出来，但却不是完美的，因为在第6）步的划分很可能会引起负载的不平衡。

    2）时间复杂度分析

     PSRS算法适合处理大批量的数据（呵呵，数据量不大，何必并行乎）。当n>p^3时，算法的时间复杂度可达n/p*logn。具体每一步的时间复杂度的分析在这里就不一一描述了，因为每一步的排序都是普通的串行排序算法。

四、算法实现

    因为算法比较复杂，代码较长，本文仅仅列出主代码，代码如下：

void psrs_mpi(int *argc, char ***argv){

    int process_id;
    int process_size;

    int *init_array;             //初始数组
    int init_array_length;         //初始数组长度

    int *local_sample;             //每个进程选取的代表元素数组
    int local_sample_length;    //代表元素数组长度

    int *sample;                 //代表元素集合（0号进程使用）
    int *sorted_sample;         //排序后的代表元素的集合
    int sample_length;             //代表预算的长度

    int *primary_sample;         //主元

    int *resp_array;             //偏移数组，主要用户指定个进程数组的各分段的长度

    int *section_resp_array;    //偏移数组,用于指定进程从其他进程获得的数组的长度

    int *section_array;             //从各个进程中获得分段数组的集合
    int *sorted_section_array;
    int section_array_length;     //总长

    int section_index;

    int i, j ;                     //循环变量

    MPI_Request handle;
    MPI_Status  status;

    mpi_start(argc, argv, &process_size, &process_id, MPI_COMM_WORLD);
    resp_array = (int *)my_mpi_malloc(process_id, sizeof(int) * process_size);

    //为每个进程构建一个数组
    //并对改进型的数组进行串行快速排序
    init_array_length = ARRAY_LENGTH;
    init_array = (int *)my_mpi_malloc(process_id, sizeof(int) * init_array_length);
    array_builder_seed(init_array, init_array_length, process_id);

    quick_sort(init_array, 0, init_array_length -1);

    //每个处理器从排序号的序列中选取process_size-1个元素
    //并发送到0号进程中
    local_sample_length = process_size - 1;
    local_sample = array_sample(init_array, local_sample_length, init_array_length/process_size, process_id);

    if(process_id)
         MPI_Send(local_sample, local_sample_length, MPI_INT, 0, SAMPLE_DATA, MPI_COMM_WORLD);


    //0号进程接收各处理器发送过来的代表元素，并将这些元素做多路归并排序
    if(!process_id){
        sample = (int *)my_mpi_malloc(0, sizeof(int) * process_size * local_sample_length);
        sorted_sample = (int *)my_mpi_malloc(0, sizeof(int) * process_size * local_sample_length);
        array_copy(sample, local_sample, local_sample_length);

        for(i = 1; i < process_size; i++)
            MPI_Irecv(sample + local_sample_length * i, local_sample_length, MPI_INT, i, SAMPLE_DATA,
                    MPI_COMM_WORLD, &handle);

        MPI_Wait(&handle, &status);

        for(i = 0; i < process_size; i++)
            resp_array[i] = local_sample_length;

        mul_merger(sample, sorted_sample, resp_array, process_size);

        //从排序好的代表元素中选取process_size-1个主元,并将这些主元广播道其他的处理器中
        primary_sample = array_sample(sorted_sample, process_size -1, process_size -1, process_id);
    }
    if(process_id)
        primary_sample = (int *)my_mpi_malloc(process_id, sizeof(int) * process_size -1);

    MPI_Bcast(primary_sample, process_size-1, MPI_INT, 0, MPI_COMM_WORLD);

    //将处理器上的数据根据主元分成process_size 端
    get_array_sepator_resp(init_array, primary_sample, resp_array, init_array_length, process_size);
    if(process_id == ID){
        printf("process %d resp array is:" ,process_id);
        array_int_print(process_size, resp_array);
    }

    //每个处理器将自己的第i段发送给第i个处理器
    section_resp_array = (int *)my_mpi_malloc(process_id, sizeof(int) * process_size);
    section_resp_array[process_id] = resp_array[process_id];

    //每个进程将要发送的数据的个数发送给哥哥处理器
    for(i = 0; i < process_size; i++){
        if(i == process_id){
            for(j = 0; j < process_size; j++)
               if(i != j)
                      MPI_Send(&(resp_array[j]), 1, MPI_INT, j, SECTION_INDEX ,
                           MPI_COMM_WORLD);
        }
        else
            MPI_Recv(&(section_resp_array[i]), 1, MPI_INT, i, SECTION_INDEX,
                    MPI_COMM_WORLD, &status);
    }

    MPI_Barrier(MPI_COMM_WORLD);

    section_array_length = get_array_element_total(section_resp_array, 0, process_size - 1);
    section_array = (int *)my_mpi_malloc(process_id, sizeof(int) * section_array_length);
    sorted_section_array = (int *)my_mpi_malloc(process_id, sizeof(int) * section_array_length);
    section_index = 0;

    for(i = 0; i < process_size; i++){
        if(i == process_id){
            for(j = 0; j < process_size; j++){
                if(j)
                    section_index = get_array_element_total(resp_array, 0 , j-1);
                if(i == j)
                    array_int_copy(section_array, init_array, section_index, section_index+resp_array[j]);
                if(i != j){
                    if(j)
                        section_index = get_array_element_total(resp_array, 0 , j-1);
                    MPI_Send(&(init_array[section_index]), resp_array[j], MPI_INT,
                            j, SECTION_DATA, MPI_COMM_WORLD);
                }
            }
        }
        else{
            if(i)
                section_index = get_array_element_total(section_resp_array, 0, i-1);
            MPI_Recv(&(section_array[section_index]), section_resp_array[i], MPI_INT,
                    i, SECTION_DATA, MPI_COMM_WORLD, &status);
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);

    //进行多路归并排序
    mul_merger(section_array, sorted_section_array, section_resp_array, process_size);

    array_int_print(section_array_length, sorted_section_array);

    //释放内存
    free(resp_array);
    free(init_array);
    free(local_sample);
    free(primary_sample);
    free(section_array);
    free(sorted_section_array);
    free(section_resp_array);

    if(!process_id){
        free(sample);
        free(sorted_sample);
    }

    MPI_Finalize();
}


五、MPI函数分析

     在上述算法中，用到了MPI的非阻塞通信函数：MPI_IRecv，其对应的是MPI_Isend。这连个函数用于进程间的非阻塞通信，使通信和运算能够同时进行。有这两
个非阻塞通信函数，就不能不提MPI_Wait函数，该函数的作用是阻塞进程执行，直到想对应的所有进程操作都执行的这个地方为止。一般是这个三函数一起使用。

    下篇，我们将介绍kmp字符串匹配算法及其并行化。