MPI并行编程系列一：枚举排序

前言

    貌似有一年多的时间没有在博客园上发过随笔了。原因有三：1、懒，这是最主要的原因。哎，人老了，心和手也都老了。， 2、时间紧，也许这个理由这叫做借口。拿没时间作为任何事情的借口都是站不住脚的，我一直是这么认为的。但当时的确是特殊情况：原本计划一年完成的项目因为特殊原因必须在两个半月内完成。曾记得那时候从早上到晚上没有任何节假日的工作时的痛苦和激情，也忘不了项目成功上线后的喜悦和成就感。孙中山先生解释什么为革命的时候，说过：要享文明之幸福，不得不经文明之痛苦，这痛苦就是革命。我深有感触，同样要享成功之喜悦，不得不经成功之痛苦，这痛苦就是奋斗。3、转行。说转行也不太合适。从事了近两年的Web系统的开发后，导师让我研究并行计算，学好MPI编程，固最近一直在MPI中游荡。说实话，我并不确定这篇文章发在博客园中是否合适，因为博客园是.net社区。但我在别的技术社区中实在是找不到博客园的感觉，所以还是发在博客园中吧。希望博客园中的并行程序设计高手给指点一二。

    以后我会从易到难写MPI程序设计的一系列的例子，开篇算法为枚举排序并行算法。

    编程环境如下：操作系统：Ubuntu， 编程语言 c， 并行库 MPI， 编译器 gcc

一、枚举排序算法说明：

    枚举排序（Enumeration Sort）是一种最为简单的排序算法，通常也被叫做秩排序（Rank Sort）。

    该算法的基本思想是：对每一个要排序的元素统计小于它的所有元素的个数，从而得到该元素在整个序列中的位置。其时间复杂度为o(n^2)。其伪代码为：

输入为：a[1], a[2] , ... , a[n]

输出为：b[1], b[2] , ..., b[n]

for i =1 to n do

    1)k =1

    2)for j = 1 to n do

          if a[i] > a[j] then

             k= k + 1

         endif

      endfor

      3)b[k] = a[i]

endfor

     算法思想很简单，现将其主要代码总结如下：

    1、数组自动生成，随机生成长度为n的数组：

   1:  void array_builder(int *a, int n)

   2:  {

   3:      int i;

   4:   

   5:      srand((int)time(0));

   6:      

   7:      for(i = 0; i < n; i++)

   8:          a[i] = (int)rand() % 100;

   9:      

  10:      return;

  11:  }

  2、取得每个元素在数组中的秩，即统计每个元素按小于它的其他所有元素的个数：

int *get_array_elem_position(int *init_array, int array_length, int start, int size){

    int i,j,k;
    int *position;

    position = (int *)my_malloc(sizeof(int) * size);
    for(i = start; i < start+size; i++){
        k = 0;
        for(j = 0; j < array_length; j++){
            if(init_array[i] < init_array[j])
                k++;
            if((init_array[i] == init_array[j]) && i >j)
                k++;
        }

        position[i-start] = k;
    }

    return position;
}
      其中my_malloc函数的作用为动态分配大小为size的空间，如分配失败，则终止程序：

   1:  void *my_malloc(int malloc_size){

   2:      void *buffer;

   3:   

   4:      if((buffer = (void *)malloc((size_t)malloc_size)) == NULL){

   5:          printf("malloc failure");

   6:          exit(EXIT_FAILURE);

   7:      }

   8:   

   9:      return buffer;

  10:  }

  3、 算法主函数：

void serial(){

    int i;
    int array_length = ARRAY_LENGTH;

    int *init_array;
    int *sort_array;
    int *position;

//    array_length = get_array_length(4);

    sort_array = (int *)my_malloc(sizeof(int) * array_length);
    init_array = (int *)my_malloc(sizeof(int) * array_length);

    array_builder(init_array, array_length);

    position = get_array_elem_position(init_array, array_length, 0, array_length);

    for(i = 0; i < array_length; i++)
        sort_array[position[i]] = init_array[i];

    printf("串行实行结果:\n");
    init_sort_array_print(init_array, sort_array, array_length);
}
       其中的数组打印函数就不列出来了。

二、枚举排序并行化分析
    该算法的并行化算法比较的简单，一般的并行化无非就为数据并行和功能并行。显然该算法适合数据并行算法。即：如果我们有n个处理器，则我们将数据平均
分为n份，让每个处理器来处理其中一份数据。该算法的描述为：

输入：a[1], ... , a[n]

输入：b[1], ... , b[n]

1)进程p0将数组广播到所有进程中

2)每个进程对其相应的数组分块计算其中每个元素的秩

3）p0收集各进程的中相应元素的秩，并给出最终的排序结果。

    算法分析：假设我们在排序中用了m个处理器，则算法的复杂度为o（n^2/m),通信复杂度为o（mn）。

三、并行算法实现

    并行程序设计中一般分为四步：划分，通信，聚集和映射。

        1、划分是将一个大任务划分为一系列小任务的过程。在本例中，计算每个元素的秩是互不相关的，因此我们可以将该任务划分为n个小任务，每个任务的作用就是

           计算其中一个元素的秩。

        2、通信就是要确定个进程的通信模式。在本例中，进程间的通信主要是进程p0将原始数组广播到各个进程中，同时收集各个进程的计算结果。

        3、聚集。如果有m个进程，我们可以将这n个小任务聚集为m个原始任务，每个进程执行其中的一个原始任务。

        4、映射。映射是将原始任务映射到各处理器的过程。这里由于没有使用集群，固映射为系统自动。

     该算法的一个重要问题就是如何对数据分块，即上面第三步中，如何经n个小任务聚集成m个原始任务。我们的数据分块策略为：按块数据分解，这就意味着我们将

数组分为m个连续的快，每个块的大小基本相等。在对数组进行分解的过程中，我们必须要确定两个问题：1、给出进程，能够算出该进程所负责的数据。2、给出数据，

能够计算出该数据被那个一进程所负责。我们将这两个问题定义在头文件中：

//n为数组长度，p为进程数，id为进城id
#define BLOCK_LOW(id, p, n) ((id)*(n)/(p))
#define BLOCK_HIGH(id, p, n) (BLOCK_LOW((id)+1, p, n)-1)
#define BLOCK_SIZE(id, p, n) (BLOCK_HIGH(id, p, n) - BLOCK_LOW(id, p, n) + 1)
#define BLOCK_OWNER(j, p, n) (((P) * ((j) +1)-1)/(n))
 并行程序为：

   1:  #include <stdlib.h>

   2:  #include <stdio.h>

   3:   

   4:  #include "mpi.h"

   5:  #include "my_mpi.h"

   6:  #include "common.h"

   7:  #include "sort/enum_sort.h"

   8:   

   9:  void  mpi_sort(

  10:          int *argc,

  11:          char ***argv){

  12:   

  13:      int process_id;

  14:      int process_size;

  15:   

  16:      int *count;

  17:      int *resp;

  18:   

  19:      int array_low;

  20:      int array_size;

  21:   

  22:      int *local_position;

  23:      int *position;

  24:   

  25:      int *sort_array;

  26:      int *init_array;

  27:      int array_length = ARRAY_LENGTH;

  28:   

  29:      int i;

  30:   

  31:      my_mpi_struct mpi_struct;

  32:  

  33:      mpi_start(argc, argv, &process_size, &process_id, MPI_COMM_WORLD);

  34:   

  35:      if(!process_id){

  36:          position = (int *)my_mpi_malloc(0, sizeof(int) * array_length);

  37:          sort_array = (int *)my_mpi_malloc(0, sizeof(int) * array_length);

  38:      }

  39:   

  40:      init_array = (int *) my_mpi_malloc(process_id, sizeof(int) * array_length);

  41:      if(!process_id)

  42:          array_builder(init_array, array_length);

  43:   

  44:      MPI_Bcast(init_array, array_length, MPI_INT, 0, MPI_COMM_WORLD);

  45:   

  46:      array_low = BLOCK_LOW(process_id, process_size, array_length);

  47:      array_size = BLOCK_SIZE(process_id, process_size, array_length);

  48:   

  49:      get_my_mpi_struct(&mpi_struct, process_id, process_size, array_length);

  50:      get_resp_count_array(&resp, &count, mpi_struct);

  51:   

  52:      local_position = get_array_elem_position(init_array, array_length, 

  53:              array_low, array_size);

  54:   

  55:      MPI_Gatherv(local_position, array_size, MPI_INT, position, count, resp,

  56:             MPI_INT, 0, MPI_COMM_WORLD);    

  57:   

  58:      if(!process_id){

  59:          for(i = 0; i< array_length; i++)

  60:              sort_array[position[i]] = init_array[i];

  61:   

  62:          printf("the result of mpi:\n");

  63:          init_sort_array_print(init_array, sort_array, array_length);

  64:   

  65:          free(sort_array);

  66:          free(position);

  67:      }

  68:   

  69:      free(count);

  70:      free(resp);

  71:      free(local_position);

  72:      free(init_array);

  73:   

  74:      MPI_Finalize();

  75:  }

  76: 

  其中mpi_start函数就是将编写MPI程序时三个必须函数MPI_Init,MPI_Comm_rank和MPI_Comm_size给写在一个函数中了，方便调用。

  函数get_resp_count_array的作用是计算MPI收集函数所需要的偏移和组大小数组，其代码如下：

void get_resp_count_array(
        int **resp,                 //the resp array
        int **count,                 //the count array
        my_mpi_struct mpi_struct){ //my mpi struct

    int i;

    int process_id = mpi_struct -> process_id;
    int process_size = mpi_struct -> process_size;
    int number_count = mpi_struct -> number_count;

    *resp = (int *)my_mpi_malloc(process_id, sizeof(int) * process_size);
    *count = (int *)my_mpi_malloc(process_id, sizeof(int) * process_size);

    (*resp)[0] = 0;
    (*count)[0] = BLOCK_SIZE(0, process_size, number_count);

    for(i = 1; i < process_size; i++){
        (*resp)[i] = (*resp)[i-1] + (*count)[i-1];
        (*count)[i] = BLOCK_SIZE(i, process_size, number_count);
    }
}

四、MPI主要函数说明

    在本例中一共用到了两个主要的MPI函数：MPI_Bcast和MPI_Gatherv。即一个广播函数和一个收集函数。

    MPI_Bcast是一个组通信操作，用来完成一个进程向通信域中所有的进程广播消息

    MPI_Gatherv是完成收集的组通信函数。跟进程从进程i收集count[i]个元素，并将手机的i个元素放在接收缓冲区的resp[i]个位置开始的地方。

五、后记

    这个MPI程序中一个比较简单的例子，作为我下MPI并行程序的一个开端。还请各位高手积极拍砖。

    下篇将对快速排序的并行算法进行简单的介绍和实现。