Parallezing merge step and divide step in mergesort algorithm with OpenMP

Question

We wrote a code that parallelizes the divide phase of the mergesort algorithm. I.e. every recursive call is assigned to a thread but our teacher was disappointed because he said we should also parallelize the merge function. I have researched how one can do this and I found these (algorithm 3) lecture notes which changes the complexity of the merge function from O(n) to O(n(log(n)) but which one can now parallelize. I wanted to ask for suggestions how to design the code for the fastest result possible:

(both should compile with g++ /there might be some warnings if one uses pedantic flag but since both yield the correct result one could ignore them)

Old code:

#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);
        //delete[] data;
        free(tmp);
        //delete[] tmp;
        free(ref);
        //delete[] ref;
    }

    return EXIT_SUCCESS;
}

New code - merge is parallelizable:

#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);
        //delete[] data;
        free(tmp);
        //delete[] tmp;
        free(ref);
        //delete[] ref;
    }
    return EXIT_SUCCESS;
}

Answer 1

/**
  * helper routine: check if array is sorted correctly
  */
bool isSorted(int ref[], int data[], const size_t size){

Review the standard library to be familiar with what all is in it. Specifically, look at sorting operations where you'll find the standard algorithm is_sorted. Don't re-implement standard library contents that's not directly part of the exercise!

As for your implementation...
The int ref[] syntax in a parameter actually means int*. You should avoid the faux array syntax. But more seriously, there are two "arrays" passed and they are not marked const?

The first line is a call to std::sort ????
So, it appears it wants you to make a copy of the data, pass both copies, then this function sorts one of the copies and checks to see if the two copies are now the same (manually re-implementing std::equals).

It is generally true that you can check your "fancy" version of a function by comparing the results against the "normal" version. But in this case, checking to see if the array is sorted is much more straightforward, and what the name of the function implies.

int *data = (int*) malloc(stSize * sizeof(int));
int *tmp = (int*) malloc(stSize * sizeof(int));     
int *ref = (int*) malloc(stSize * sizeof(int));
printf("Initialization...\n");

Why are you using malloc (and printf)?

data[idx] = (int) (stSize * (double(rand()) / RAND_MAX));

One might be more forgiving of rand() since we don't need particularly high-quality randomness here; but you are doing all the conditioning work and explicit casting to get a random int. Use the C++ random number library, and simply ask it for random integers in the desired range.

Note that C++ has a time library, also. You seem to be using C, not C++, except for calling std::sort and std::copy.

The other call to sort, I see, is:

bool myfunc (long i , long j){return (i<j);}
  ...
std::sort(array+begin,array + end, myfunc);

Passing in myfunc is slowing it down. Simple ascending via < is the default, so you don't need to pass this at all. But had you actually needed to supply a comparison function, a pointer to a plain function is de-optimizing because the optimizer generally doesn't follow pointers to functions even when they are known at compile time. That is, the call to myfunc won't be inlined.

Your function takes long when you are sorting int. The function this appears in is declaring begin and end as long, rather than size_t. It's like you spliced in parts from different programs, or forgot what types you were using.

Answer 2

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>

#include <cstdlib>
#include <cstdio>

Prefer the C++ header files, which declare their identifiers in the std namespace. There's no need to include the old C-compatibility headers when writing new C++ code (they have some limited use, when writing headers for extern "C" functions, but leave them alone for normal code).

---

bool isSorted(int ref[], int data[], const size_t size){

Oops - that's a reserved name (is followed by a letter) in the global namespace.

Things you could do to fix that:

• change the name
• move it to another namespace (perhaps an anonymous namespace)
• change it to static linkage (and check that neither <ctype.h> nor <wctype.h> is included)

---

        printf("\nFalscher Index:%d\n",idx);

That looks like an error message, so should go to stderr. It's probably better to use C++ iostreams (std::cerr) than C-style std::printf() here, though.

---

    const size_t stSize = strtol(argv[1], NULL, 10);

We should be checking that conversion was successful, and also that it resulted in a positive number. Why are we using a signed conversion instead of std::strtoul() for assignment to (unsigned) size_t? I think that std::stul() is a better choice, as it gives improved error reporting compared to the C-style conversion.

---

Why are you using C-style std::malloc() and std::free()? C++ has better memory allocation facilities. Prefer (in priority order):

• a standard container such as std::vector,
• a smart pointer, from std::make_unique() or std::make_shared(),
• (last resort) new/delete or new[]/delete[].

Answer 3

changes the complexity of the merge function from O(n) to O(n(log(n)) but which one can now parallelize

For a merge of two sorted runs of size n, worst case for each instance of a binary search is ⌈log2(n)⌉ reads, and for larger n (>=256 or so), this will result in the parallelized binary search merge being slower than a single threaded conventional merge.

A faster approach for t threads would be to split up the array into t parts, and use t threads to independently sort each of the t parts. Then merge even and odd pairs of sorted runs, the first merge using t/2 threads, the next merge using t/4 threads, and finally a single thread to merge two sorted runs. An example of this using Windows native threads is shown in this question:

Multithreaded bottom up merge sort