MPI实现psrs

PSRS(Parallel Sorting by Regular Sampling,并行正则采样排序)的基本原理如下：
    begin
        (1)均匀划分：将n个元素A[1..n]均匀划分成p段，每个pi处理A[(i-1)n/p+1..in/p]
        (2)局部排序：pi调用串行排序算法对A[(i-1)n/p+1..in/p]排序
        (3)选取样本：pi从其有序子序列A[(i-1)n/p+1..in/p]中选取p个样本元素
        (4)样本排序：用一台处理器对p2个样本元素进行串行排序
        (5)选择主元：用一台处理器从排好序的样本序列中选取p-1个主元，并播送给其他pi
        (6)主元划分： pi按主元将有序段A[(i-1)n/p+1..in/p]划分成p段
        (7)全局交换：各处理器将其有序段按段号交换到对应的处理器中

        (8)归并排序：各处理器对接收到的元素进行归并排序

基于MPI的并行设计

1、局部排序，并选取样本

  qsort(array + startIndex, subArraySize, sizeof(array[0]), cmp);

  // 正则采样
  for (i = 0; i < p; i++) {
    //pivots[i] = array[startIndex + ((i+1) * (N / (p * (p+1))))];  
    pivots[i] = array[startIndex + (i * (N / (p * p)))];	
  }

2、样本排序，选择主元，并将其发送给其他进程

  double *collectedPivots = (double *) malloc(p * p * sizeof(pivots[0]));
  double *phase2Pivots = (double *) malloc((p - 1) * sizeof(pivots[0]));          //主元
  int index = 0;
  MPI_Barrier(MPI_COMM_WORLD);
  //收集消息，根进程在它的接受缓冲区中包含所有进程的发送缓冲区的连接。
  MPI_Gather(pivots, p, MPI_DOUBLE, collectedPivots, p, MPI_DOUBLE, 0, MPI_COMM_WORLD);       
  if (myId == 0) {

    qsort(collectedPivots, p * p, sizeof(pivots[0]), cmp);          //对正则采样的样本进行排序

    // 采样排序后进行主元的选择
    for (i = 0; i < (p -1); i++) {
      phase2Pivots[i] = collectedPivots[(((i+1) * p) + (p / 2)) - 1];
    }
  }
  //发送广播
  MPI_Bcast(phase2Pivots, p - 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

3、根据主元进行划分

  // 进行主元划分，并计算划分部分的大小
  for ( i = 0; i < subArraySize; i++) {
    if (array[startIndex + i] > phase2Pivots[index]) {
      //如果当前位置的数字大小超过主元位置，则进行下一个划分
      index += 1;
    }
    if (index == p-1) {
      //最后一次划分，子数组总长减掉当前位置即可得到最后一个子数组划分的大小
      partitionSizes[p - 1] = subArraySize - i;
      break;
    }
    partitionSizes[index]++ ;   //划分大小自增
  }

4、全局交换

int *sendDisp = (int *) malloc(p * sizeof(int));
  int *recvDisp = (int *) malloc(p * sizeof(int));
  MPI_Barrier(MPI_COMM_WORLD);
  // 全局到全局的发送，每个进程可以向每个接收者发送数目不同的数据.
  MPI_Alltoall(partitionSizes, 1, MPI_INT, newPartitionSizes, 1, MPI_INT, MPI_COMM_WORLD);

  // 计算划分的总大小，并给新划分分配空间
  for ( i = 0; i < p; i++) {
    totalSize += newPartitionSizes[i];
  }
  *newPartitions = (double *) malloc(totalSize * sizeof(double));

  // 在发送划分之前计算相对于sendbuf的位移，此位移处存放着输出到进程的数据
  sendDisp[0] = 0;
  recvDisp[0] = 0;      //计算相对于recvbuf的位移，此位移处存放着从进程接受到的数据
  for ( i = 1; i < p; i++) {
    sendDisp[i] = partitionSizes[i - 1] + sendDisp[i - 1];
    recvDisp[i] = newPartitionSizes[i - 1] + recvDisp[i - 1];
  }
  MPI_Barrier(MPI_COMM_WORLD);
  //发送数据，实现n次点对点通信
  MPI_Alltoallv(&(array[startIndex]), partitionSizes, sendDisp, MPI_DOUBLE, 
	*newPartitions, newPartitionSizes, recvDisp, MPI_DOUBLE, MPI_COMM_WORLD);

5、归并排序

  // 归并排序
  for ( i = 0; i < totalListSize; i++) {
    double lowest = 134217728.0;
    int ind = -1;
    for (j = 0; j < p; j++) {
      if ((indexes[j] < partitionEnds[j]) && (partitions[indexes[j]] < lowest)) {
    lowest = partitions[indexes[j]];
    ind = j;
      }
    }
    sortedSubList[i] = lowest;
    indexes[ind] += 1;
  }

6、将排序好的子列表发送给根进程

MPI_Barrier(MPI_COMM_WORLD);
  // 发送各子列表的大小回根进程中
  MPI_Gather(&totalListSize, 1, MPI_INT, subListSizes, 1, MPI_INT, 0, MPI_COMM_WORLD);
 
  // 计算根进程上的相对于recvbuf的偏移量
  if (myId == 0) {
    recvDisp[0] = 0;
    for ( i = 1; i < p; i++) {
      recvDisp[i] = subListSizes[i - 1] + recvDisp[i - 1];
    }
  }
  
  MPI_Barrier(MPI_COMM_WORLD);
  //发送各排好序的子列表回根进程中
  MPI_Gatherv(sortedSubList, totalListSize, MPI_DOUBLE, array, subListSizes, 
	recvDisp, MPI_DOUBLE, 0, MPI_COMM_WORLD);

7、验证排序结果

if (myId == 0) {
		finish = MPI_Wtime();
		total = finish - start;
		int ver = 1;
		for( i = 0; i < N-1; i++){
			if(array[i] > array[i+1]){
				ver = 0;
				break;
			}
		}
		 if(ver){
			printf("排序成功\n"); 
		 } else {
			printf("排序失败\n"); 
		 }
		printf("总共用时:%lfs\n",total);
	}

运行结果

进程数	1	2	4	8	16
加速比	1.000	1.513	3.437	6.192	7.918