基于OpenMP的归并排序划分与归并阶段并行化最优实现咨询

我们编写了对归并排序（mergesort）的划分阶段进行并行化的代码，即每个递归调用分配给一个线程，但老师指出还需要对归并函数进行并行化。查阅资料后找到斯坦福大学的相关讲义（算法3），该方案将归并函数的复杂度从O(n)调整为O(nlog(n))，支持并行化。现有两份均可通过g++编译的实现代码（开启pedantic编译选项可能产生警告，但运行结果正确可忽略），特此咨询如何设计代码以获得最快的运行性能：

旧代码（仅划分阶段并行） #

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <sys/time.h>

#include <iostream>
#include <algorithm>

#include <cstdlib>
#include <cstdio>

#include <cmath>
#include <ctime>
#include <cstring>
#include <omp.h>
// Constants.h
#if !defined(MYLIB_CONSTANTS_H)
#define MYLIB_CONSTANTS_H 1

const int CUTOFF =11;

#endif


/**
  * helper routine: check if array is sorted correctly
  */
bool isSorted(int ref[], int data[], const size_t size){
    std::sort(ref, ref + size);
    for (size_t idx = 0; idx < size; ++idx){
        if (ref[idx] != data[idx]) {
            return false;
        }
    }
    return true;
}

/**
  * sequential merge step (straight-forward implementation)
  */
void MsMergeSequential(int *out, int *in, long begin1, long end1, long begin2, long end2, long outBegin) {
    long left = begin1;
    long right = begin2;
    long idx = outBegin;
    
    while (left < end1 && right < end2) {
        if (in[left] <= in[right]) {
            out[idx] = in[left];    
            left++;
        } else {
            out[idx] = in[right];
            right++;
        }
        idx++;
        }
    
    while (left < end1) {
        out[idx] = in[left];
        left++, idx++;
    }
    while (right < end2) {
        out[idx] = in[right];
        right++, idx++;
    }
}
bool myfunc (long i , long j){return (i<j);}
/**
  * sequential MergeSort
  */
void MsSequential(int *array, int *tmp, bool inplace, long begin, long end) {
  if( end <= begin + CUTOFF -1){

    std::sort(array+begin,array + end, myfunc);
  }
  else  if (begin < (end - 1)) {
           long half =(begin+end) / 2;


            #pragma omp taskgroup
         {
           #pragma omp task shared(array) untied if(end-begin >= (1<<15))

             MsSequential(array, tmp, !inplace, begin, half);

             MsSequential(array, tmp, !inplace, half, end);
              }
 if (inplace){
            MsMergeSequential(array, tmp, begin, half, half, end, begin);
 } else {
            MsMergeSequential(tmp, array, begin, half, half, end, begin);
 }
        
    } else if (!inplace) {

        tmp[begin] = array[begin];
    }
}

/**
  * Serial MergeSort
  */
void MsSerial(int *array, int *tmp, const size_t size) {

    MsSequential(array, tmp, true, 0, size);
}

/**
  * @brief program entry point
  */
int main(int argc, char* argv[]) {
    // variables to measure the elapsed time
    struct timeval t1, t2;
    double etime;

    // expect one command line arguments: array size
    if (argc != 2) {
        printf("Usage: MergeSort.exe <array size> \n");
        printf("\n");
        return EXIT_FAILURE;
    }
    else {
        const size_t stSize = strtol(argv[1], NULL, 10);
        int *data = (int*) malloc(stSize * sizeof(int));
        int *tmp = (int*) malloc(stSize * sizeof(int));     
        int *ref = (int*) malloc(stSize * sizeof(int));
        printf("Initialization...\n");

        srand(95);

        #pragma omp parallel for num_threads(100) schedule(static)
        for (size_t idx = 0; idx < stSize; ++idx){
            data[idx] = (int) (stSize * (double(rand()) / RAND_MAX));
        }
        std::copy(data, data + stSize, ref);

        double dSize = (stSize * sizeof(int)) / 1024 / 1024;
        printf("Sorting %zu elements of type int (%f MiB)...\n", stSize, dSize);

        gettimeofday(&t1, NULL);
        #pragma omp parallel num_threads(80) 
        {
        #pragma omp single
        {
        MsSerial(data, tmp, stSize);
        }
        }
        gettimeofday(&t2, NULL);
        etime = (t2.tv_sec - t1.tv_sec) * 1000 + (t2.tv_usec - t1.tv_usec) / 1000;
        etime = etime / 1000;

        printf("done, took %f sec. Verification...", etime);
        if (isSorted(ref, data, stSize)) {
            printf(" successful.\n");
        }
        else {
            printf(" FAILED.\n");
        }

        free(data);
        free(tmp);
        free(ref);
    }

    return EXIT_SUCCESS;
}

新代码（归并阶段支持并行） #

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <sys/time.h>

#include <iostream>
#include <algorithm>

#include <cstdlib>
#include <cstdio>

#include <cmath>
#include <ctime>
#include <cstring>
#include <omp.h>
// Constants.h
#if !defined(MYLIB_CONSTANTS_H)
#define MYLIB_CONSTANTS_H 1

#endif


//Takes a sorted list of size n and a value, puts the value in one of n+1 possible positions,
//if value is same to an element of the list take the position before the first occurence of the same element

int binarysearchfindlowerrank(int *in,int n,int value,int projection){

    int* array= in+projection;
    int L=0;
    int R=n;
    while(R-L>1){
        int middle = (R+L)/2;
        if(array[middle]==value){
            while(array[middle]==value&&middle>0){
                middle=middle-1;
            }
            if(middle==0&&array[middle]>=value){
                return 0;
            }
            else{
            return middle+1;
            }
        }
        if(array[middle]<value){
            L=middle;
        }
        if(array[middle]>value){
            R=middle;
        }
    }
    if(n==1){
        if(array[0]>=value){
            return 0;
        }
        else return 1;
    }
    if(L==0&&array[L]>value){
        return 0;
    }
    if(R==n && array[R-1]< value){
        return n;
    }
    if(R==n&& array[R-1]>=value){
        return R-1;
    }
    if(array[R]<value){
        return R+1;
    }
    if(array[L]<value){
        return R;
    }
    return L;
}


//Takes a sorted list of size n and a value, puts the value in one of n+1 possible positions,
//if value is same to an element of the list take the position after the last occurence of the same element


int binarysearchfinduperrank(int *in,int n,int value, int projection){

    int* array= in+projection;
    int L=0;
    int R=n;
    while(R-L>1){
        int middle = (R+L)/2;
        if(array[middle]==value){
            while(array[middle]==value&&middle<n){
                middle=middle+1;
            }
            return middle;
        }
        if(array[middle]<value){
            L=middle;
        }
        if(array[middle]>value){
            R=middle;
        }
    }
     if(n==1){
         if(array[0]> value){
             return 0;
        }
        else{
            return 1;
        }
    }
    if(L==0&&array[L]>value){
        return 0;
    }
    if(R==n && array[R-1]<= value){
        return n;
    }
    if(R==n&& array[R-1]>value){
        return R-1;
    }
    if(array[R]<=value){
        return R+1;
    }
    if(array[L]<=value){
        return R;
    }
    return L;
}

/**
  * helper routine: check if array is sorted correctly
  */
bool isSorted(int ref[], int data[], const size_t size){
    std::sort(ref, ref + size);
    for (size_t idx = 0; idx < size; ++idx){
        if (ref[idx] != data[idx]) {
            printf("\nFalscher Index:%d\n",idx);
            return false;
        }
    }
    return true;
}

/**
  * sequential merge step (straight-forward implementation)
  */
void MsMergeParallelized(int *out, int *in, long begin1, long end1, long begin2, long end2, long outBegin,int *data,int *tmp) {
    if(begin1==end2){
        out[begin1]=in[begin1];
    }
    if(begin1==begin2||begin2==end2){
        out[begin1+binarysearchfinduperrank(in,1,in[end2],begin1)]=in[end2];
        out[begin1+binarysearchfindlowerrank(in,1,in[begin1],end2)]=in[begin1];
    }
    else{
        for(int i=0;i<(end2-begin2);i++){
            out[begin1+i+binarysearchfinduperrank(in,(end1-begin1),in[begin2+i],begin1)]=in[begin2+i];
        }
        for(int i=0;i<(end1-begin1);i++){
            out[begin1+i+binarysearchfindlowerrank(in,(end2-begin2),in[begin1+i],begin2)]=in[begin1+i];
        }
    }
}
bool myfunc (long i , long j){return (i<j);}
/**
  * sequential MergeSort
  */
void MsParallelized(int *array, int *tmp, bool inplace, long begin, long end) {
  if (begin < (end - 1)) {
        long half =(begin+end) / 2;
        MsParallelized(array, tmp, !inplace, begin, half);
        MsParallelized(array, tmp, !inplace, half, end);
        if (inplace){
            MsMergeParallelized(array, tmp, begin, half, half, end, begin,array,tmp);
        }
        else {
            MsMergeParallelized(tmp, array, begin, half, half, end, begin,array,tmp);
        }
    }
    else if (!inplace) {
        tmp[begin] = array[begin];
    }
}

/**
  * Serial MergeSort
  */
void MsParallel(int *array, int *tmp, const size_t size) {

    MsParallelized(array, tmp, true, 0, size);
}

/**
  * @brief program entry point
  */
int main(int argc, char* argv[]) {
    // variables to measure the elapsed time
    struct timeval t1, t2;
    double etime;

    // expect one command line arguments: array size
    if (argc != 2) {
        printf("Usage: MergeSort.exe <array size> \n");
        printf("\n");
        return EXIT_FAILURE;
    }
    else {
        const size_t stSize = strtol(argv[1], NULL, 10);
        int *data = (int*) malloc(stSize * sizeof(int));
        int *tmp = (int*) malloc(stSize * sizeof(int));
        int *ref = (int*) malloc(stSize * sizeof(int));
        printf("Initialization...\n");

        srand(95);


        for (size_t idx = 0; idx < stSize; ++idx){
            data[idx] = (int) (stSize * (double(rand()) / RAND_MAX));
        }
        std::copy(data, data + stSize, ref);
        double dSize = (stSize * sizeof(int)) / 1024 / 1024;
        printf("Sorting %u elements of type int (%f MiB)...\n", stSize, dSize);
        gettimeofday(&t1, NULL);
        // Mergesort starts
        MsParallel(data, tmp, stSize);

        gettimeofday(&t2, NULL);
        etime = (t2.tv_sec - t1.tv_sec) * 1000 + (t2.tv_usec - t1.tv_usec) / 1000;
        etime = etime / 1000;
        printf("done, took %f sec. Verification...", etime);
        if (isSorted(ref, data, stSize)) {
            printf(" successful.\n");
        }
        else {
            printf(" FAILED.\n");
        }
        free(data);
        free(tmp);
        free(ref);
    }
    return EXIT_SUCCESS;
}

内容的提问来源于stack exchange，提问作者New2Math