Merge sort

Example C-like code using indices for top-down merge sort algorithm that recursively splits the list (called runs in this example) into sublists until sublist size is 1, then merges those sublists to produce a sorted list. The copy back step is avoided with alternating the direction of the merge with each level of recursion (except for an initial one-time copy, that can be avoided too). To help understand this, consider an array with two elements. The elements are copied to B[], then merged back to A[]. If there are four elements, when the bottom of the recursion level is reached, single element runs from A[] are merged to B[], and then at the next higher level of recursion, those two-element runs are merged to A[]. This pattern continues with each level of recursion.

// Array A[] has the items to sort; array B[] is a work array.
void TopDownMergeSort(A[], B[], n)
{
    CopyArray(A, 0, n, B);           // one time copy of A[] to B[]
    TopDownSplitMerge(A, 0, n, B);   // sort data from B[] into A[]
}

// Split A[] into 2 runs, sort both runs into B[], merge both runs from B[] to A[]
// iBegin is inclusive; iEnd is exclusive (A[iEnd] is not in the set).
void TopDownSplitMerge(B[], iBegin, iEnd, A[])
{
    if (iEnd - iBegin <= 1)                     // if run size == 1
        return;                                 //   consider it sorted
    // split the run longer than 1 item into halves
    iMiddle = (iEnd + iBegin) / 2;              // iMiddle = mid point
    // recursively sort both runs from array A[] into B[]
    TopDownSplitMerge(A, iBegin,  iMiddle, B);  // sort the left  run
    TopDownSplitMerge(A, iMiddle,    iEnd, B);  // sort the right run
    // merge the resulting runs from array B[] into A[]
    TopDownMerge(B, iBegin, iMiddle, iEnd, A);
}

//  Left source half is A[ iBegin:iMiddle-1].
// Right source half is A[iMiddle:iEnd-1   ].
// Result is            B[ iBegin:iEnd-1   ].
void TopDownMerge(B[], iBegin, iMiddle, iEnd, A[])
{
    i = iBegin, j = iMiddle;
 
    // While there are elements in the left or right runs...
    for (k = iBegin; k < iEnd; k++) {
        // If left run head exists and is <= existing right run head.
        if (i < iMiddle && (j >= iEnd || A[i] <= A[j])) {
            B[k] = A[i];
            i = i + 1;
        } else {
            B[k] = A[j];
            j = j + 1;
        }
    }
}

void CopyArray(A[], iBegin, iEnd, B[])
{
    for (k = iBegin; k < iEnd; k++)
        B[k] = A[k];
}

Sorting the entire array is accomplished by TopDownMergeSort(A, B, length(A)).

Example C-like code using indices for bottom-up merge sort algorithm which treats the list as an array of n sublists (called runs in this example) of size 1, and iteratively merges sub-lists back and forth between two buffers:

// array A[] has the items to sort; array B[] is a work array
void BottomUpMergeSort(A[], B[], n)
{
    // Each 1-element run in A is already "sorted".
    // Make successively longer sorted runs of length 2, 4, 8, 16... until the whole array is sorted.
    for (width = 1; width < n; width = 2 * width)
    {
        // Array A is full of runs of length width.
        for (i = 0; i < n; i = i + 2 * width)
        {
            // Merge two runs: A[i:i+width-1] and A[i+width:i+2*width-1] to B[]
            // or copy A[i:n-1] to B[] ( if (i+width >= n) )
            BottomUpMerge(A, i, min(i+width, n), min(i+2*width, n), B);
        }
        // Now work array B is full of runs of length 2*width.
        // Copy array B to array A for the next iteration.
        // A more efficient implementation would swap the roles of A and B.
        CopyArray(B, A, n);
        // Now array A is full of runs of length 2*width.
    }
}

//  Left run is A[iLeft :iRight-1].
// Right run is A[iRight:iEnd-1  ].
void BottomUpMerge(A[], iLeft, iRight, iEnd, B[])
{
    i = iLeft, j = iRight;
    // While there are elements in the left or right runs...
    for (k = iLeft; k < iEnd; k++) {
        // If left run head exists and is <= existing right run head.
        if (i < iRight && (j >= iEnd || A[i] <= A[j])) {
            B[k] = A[i];
            i = i + 1;
        } else {
            B[k] = A[j];
            j = j + 1;    
        }
    } 
}

void CopyArray(B[], A[], n)
{
    for (i = 0; i < n; i++)
        A[i] = B[i];
}

Pseudocode for top-down merge sort algorithm which recursively divides the input list into smaller sublists until the sublists are trivially sorted, and then merges the sublists while returning up the call chain.

function merge_sort(list m) is
    // Base case. A list of zero or one elements is sorted, by definition.
    if length of m ≤ 1 then
        return m

    // Recursive case. First, divide the list into equal-sized sublists
    // consisting of the first half and second half of the list.
    // This assumes lists start at index 0.
    var left := empty list
    var right := empty list
    for each x with index i in m do
        if i < (length of m)/2 then
            add x to left
        else
            add x to right

    // Recursively sort both sublists.
    left := merge_sort(left)
    right := merge_sort(right)

    // Then merge the now-sorted sublists.
    return merge(left, right)

In this example, the merge function merges the left and right sublists.

function merge(left, right) is
    var result := empty list

    while left is not empty and right is not empty do
        if first(left) ≤ first(right) then
            append first(left) to result
            left := rest(left)
        else
            append first(right) to result
            right := rest(right)

    // Either left or right may have elements left; consume them.
    // (Only one of the following loops will actually be entered.)
    while left is not empty do
        append first(left) to result
        left := rest(left)
    while right is not empty do
        append first(right) to result
        right := rest(right)
    return result

Pseudocode for bottom-up merge sort algorithm which uses a small fixed size array of references to nodes, where array[i] is either a reference to a list of size 2ⁱ or nil. node is a reference or pointer to a node. The merge() function would be similar to the one shown in the top-down merge lists example, it merges two already sorted lists, and handles empty lists. In this case, merge() would use node for its input parameters and return value.

function merge_sort(node head) is
    // return if empty list
    if head = nil then
        return nil
    var node array[32]; initially all nil
    var node result
    var node next
    var int  i
    result := head
    // merge nodes into array
    while result ≠ nil do
        next := result.next;
        result.next := nil
        for (i = 0; (i < 32) && (array[i] ≠ nil); i += 1) do
            result := merge(array[i], result)
            array[i] := nil
        // do not go past end of array
        if i = 32 then
            i -= 1
        array[i] := result
        result := next
    // merge array into single list
    result := nil
    for (i = 0; i < 32; i += 1) do
        result := merge(array[i], result)
    return result

The sequential merge sort procedure can be described in two phases, the divide phase and the merge phase. The first consists of many recursive calls that repeatedly perform the same division process until the subsequences are trivially sorted (containing one or no element). An intuitive approach is the parallelization of those recursive calls.[18] Following pseudocode describes the merge sort with parallel recursion using the fork and join keywords:

// Sort elements lo through hi (exclusive) of array A.
algorithm mergesort(A, lo, hi) is
    if lo+1 < hi then  // Two or more elements.
        mid := ⌊(lo + hi) / 2⌋
        fork mergesort(A, lo, mid)
        mergesort(A, mid, hi)
        join
        merge(A, lo, mid, hi)

This algorithm is the trivial modification of the sequential version and does not parallelize well. Therefore, its speedup is not very impressive. It has a span of ${\displaystyle \Theta (n)}$, which is only an improvement of ${\displaystyle \Theta (\log n)}$ compared to the sequential version (see Introduction to Algorithms). This is mainly due to the sequential merge method, as it is the bottleneck of the parallel executions.

Better parallelism can be achieved by using a parallel merge algorithm. Cormen et al. present a binary variant that merges two sorted sub-sequences into one sorted output sequence.[18]

In one of the sequences (the longer one if unequal length), the element of the middle index is selected. Its position in the other sequence is determined in such a way that this sequence would remain sorted if this element were inserted at this position. Thus, one knows how many other elements from both sequences are smaller and the position of the selected element in the output sequence can be calculated. For the partial sequences of the smaller and larger elements created in this way, the merge algorithm is again executed in parallel until the base case of the recursion is reached.

The following pseudocode shows the modified parallel merge sort method using the parallel merge algorithm (adopted from Cormen et al.).

/**
 * A: Input array
 * B: Output array
 * lo: lower bound
 * hi: upper bound
 * off: offset
 */
algorithm parallelMergesort(A, lo, hi, B, off) is
    len := hi - lo + 1
    if len == 1 then
        B[off] := A[lo]
    else let T[1..len] be a new array
        mid := ⌊(lo + hi) / 2⌋ 
        mid' := mid - lo + 1
        fork parallelMergesort(A, lo, mid, T, 1)
        parallelMergesort(A, mid + 1, hi, T, mid' + 1) 
        join 
        parallelMerge(T, 1, mid', mid' + 1, len, B, off)

In order to analyze a recurrence relation for the worst case span, the recursive calls of parallelMergesort have to be incorporated only once due to their parallel execution, obtaining

$${\displaystyle T_{\infty }^{\text{sort}}(n)=T_{\infty }^{\text{sort}}\left({\frac {n}{2}}\right)+T_{\infty }^{\text{merge}}(n)=T_{\infty }^{\text{sort}}\left({\frac {n}{2}}\right)+\Theta \left(\log(n)^{2}\right).}$$

For detailed information about the complexity of the parallel merge procedure, see Merge algorithm.

The solution of this recurrence is given by

$${\displaystyle T_{\infty }^{\text{sort}}=\Theta \left(\log(n)^{3}\right).}$$

This parallel merge algorithm reaches a parallelism of ${\textstyle \Theta \left({\frac {n}{(\log n)^{2}}}\right)}$, which is much higher than the parallelism of the previous algorithm. Such a sort can perform well in practice when combined with a fast stable sequential sort, such as insertion sort, and a fast sequential merge as a base case for merging small arrays.[19]

It seems arbitrary to restrict the merge sort algorithms to a binary merge method, since there are usually p > 2 processors available. A better approach may be to use a K-way merge method, a generalization of binary merge, in which ${\displaystyle k}$ sorted sequences are merged. This merge variant is well suited to describe a sorting algorithm on a PRAM.[20][21]

The parallel multiway mergesort process on four processors ${\displaystyle t_{0}}$ to ${\displaystyle t_{3}}$.

Given an unsorted sequence of ${\displaystyle n}$ elements, the goal is to sort the sequence with ${\displaystyle p}$ available processors. These elements are distributed equally among all processors and sorted locally using a sequential Sorting algorithm. Hence, the sequence consists of sorted sequences ${\displaystyle S_{1},...,S_{p}}$ of length ${\textstyle \lceil {\frac {n}{p}}\rceil }$. For simplification let ${\displaystyle n}$ be a multiple of ${\displaystyle p}$, so that ${\textstyle \left\vert S_{i}\right\vert ={\frac {n}{p}}}$ for ${\displaystyle i=1,...,p}$.

These sequences will be used to perform a multisequence selection/splitter selection. For ${\displaystyle j=1,...,p}$, the algorithm determines splitter elements ${\displaystyle v_{j}}$ with global rank ${\textstyle k=j{\frac {n}{p}}}$. Then the corresponding positions of ${\displaystyle v_{1},...,v_{p}}$ in each sequence ${\displaystyle S_{i}}$ are determined with binary search and thus the ${\displaystyle S_{i}}$ are further partitioned into ${\displaystyle p}$ subsequences ${\displaystyle S_{i,1},...,S_{i,p}}$ with ${\textstyle S_{i,j}:=\{x\in S_{i}|rank(v_{j-1})<rank(x)\leq rank(v_{j})\}}$.

Furthermore, the elements of ${\displaystyle S_{1,i},...,S_{p,i}}$ are assigned to processor ${\displaystyle i}$, means all elements between rank ${\textstyle (i-1){\frac {n}{p}}}$ and rank ${\textstyle i{\frac {n}{p}}}$, which are distributed over all ${\displaystyle S_{i}}$. Thus, each processor receives a sequence of sorted sequences. The fact that the rank ${\displaystyle k}$ of the splitter elements ${\displaystyle v_{i}}$ was chosen globally, provides two important properties: On the one hand, ${\displaystyle k}$ was chosen so that each processor can still operate on ${\textstyle n/p}$ elements after assignment. The algorithm is perfectly load-balanced. On the other hand, all elements on processor ${\displaystyle i}$ are less than or equal to all elements on processor ${\displaystyle i+1}$. Hence, each processor performs the p-way merge locally and thus obtains a sorted sequence from its sub-sequences. Because of the second property, no further p-way-merge has to be performed, the results only have to be put together in the order of the processor number.

In its simplest form, given ${\displaystyle p}$ sorted sequences ${\displaystyle S_{1},...,S_{p}}$ distributed evenly on ${\displaystyle p}$ processors and a rank ${\displaystyle k}$, the task is to find an element ${\displaystyle x}$ with a global rank ${\displaystyle k}$ in the union of the sequences. Hence, this can be used to divide each ${\displaystyle S_{i}}$ in two parts at a splitter index ${\displaystyle l_{i}}$, where the lower part contains only elements which are smaller than ${\displaystyle x}$, while the elements bigger than ${\displaystyle x}$ are located in the upper part.

The presented sequential algorithm returns the indices of the splits in each sequence, e.g. the indices ${\displaystyle l_{i}}$ in sequences ${\displaystyle S_{i}}$ such that ${\displaystyle S_{i}[l_{i}]}$ has a global rank less than ${\displaystyle k}$ and ${\displaystyle \mathrm {rank} \left(S_{i}[l_{i}+1]\right)\geq k}$.[22]

algorithm msSelect(S : Array of sorted Sequences [S_1,..,S_p], k : int) is
    for i = 1 to p do 
	(l_i, r_i) = (0, |S_i|-1)
	
    while there exists i: l_i < r_i do
	// pick Pivot Element in S_j[l_j], .., S_j[r_j], chose random j uniformly
	v := pickPivot(S, l, r)
	for i = 1 to p do 
	    m_i = binarySearch(v, S_i[l_i, r_i]) // sequentially
	if m_1 + ... + m_p >= k then // m_1+ ... + m_p is the global rank of v
	    r := m  // vector assignment
	else
	    l := m 
	
    return l

For the complexity analysis the PRAM model is chosen. If the data is evenly distributed over all ${\displaystyle p}$, the p-fold execution of the binarySearch method has a running time of ${\displaystyle {\mathcal {O}}\left(p\log \left(n/p\right)\right)}$. The expected recursion depth is ${\displaystyle {\mathcal {O}}\left(\log \left(\textstyle \sum _{i}|S_{i}|\right)\right)={\mathcal {O}}(\log(n))}$ as in the ordinary Quickselect. Thus the overall expected running time is ${\displaystyle {\mathcal {O}}\left(p\log(n/p)\log(n)\right)}$.

Applied on the parallel multiway merge sort, this algorithm has to be invoked in parallel such that all splitter elements of rank ${\textstyle i{\frac {n}{p}}}$ for ${\displaystyle i=1,..,p}$ are found simultaneously. These splitter elements can then be used to partition each sequence in ${\displaystyle p}$ parts, with the same total running time of ${\displaystyle {\mathcal {O}}\left(p\,\log(n/p)\log(n)\right)}$.

Below, the complete pseudocode of the parallel multiway merge sort algorithm is given. We assume that there is a barrier synchronization before and after the multisequence selection such that every processor can determine the splitting elements and the sequence partition properly.

/**
 * d: Unsorted Array of Elements
 * n: Number of Elements
 * p: Number of Processors
 * return Sorted Array
 */
algorithm parallelMultiwayMergesort(d : Array, n : int, p : int) is
    o := new Array[0, n]                         // the output array
    for i = 1 to p do in parallel                // each processor in parallel
        S_i := d[(i-1) * n/p, i * n/p] 	         // Sequence of length n/p
	sort(S_i)                                // sort locally
        synch
	v_i := msSelect([S_1,...,S_p], i * n/p)          // element with global rank i * n/p
        synch
	(S_i,1, ..., S_i,p) := sequence_partitioning(si, v_1, ..., v_p) // split s_i into subsequences
	    
	o[(i-1) * n/p, i * n/p] := kWayMerge(s_1,i, ..., s_p,i)  // merge and assign to output array
	
    return o

Firstly, each processor sorts the assigned ${\displaystyle n/p}$ elements locally using a sorting algorithm with complexity ${\displaystyle {\mathcal {O}}\left(n/p\;\log(n/p)\right)}$. After that, the splitter elements have to be calculated in time ${\displaystyle {\mathcal {O}}\left(p\,\log(n/p)\log(n)\right)}$. Finally, each group of ${\displaystyle p}$ splits have to be merged in parallel by each processor with a running time of ${\displaystyle {\mathcal {O}}(\log(p)\;n/p)}$ using a sequential p-way merge algorithm. Thus, the overall running time is given by

${\displaystyle {\mathcal {O}}\left({\frac {n}{p}}\log \left({\frac {n}{p}}\right)+p\log \left({\frac {n}{p}}\right)\log(n)+{\frac {n}{p}}\log(p)\right)}$.

The multiway merge sort algorithm is very scalable through its high parallelization capability, which allows the use of many processors. This makes the algorithm a viable candidate for sorting large amounts of data, such as those processed in computer clusters. Also, since in such systems memory is usually not a limiting resource, the disadvantage of space complexity of merge sort is negligible. However, other factors become important in such systems, which are not taken into account when modelling on a PRAM. Here, the following aspects need to be considered: Memory hierarchy, when the data does not fit into the processors cache, or the communication overhead of exchanging data between processors, which could become a bottleneck when the data can no longer be accessed via the shared memory.

Sanders et al. have presented in their paper a bulk synchronous parallel algorithm for multilevel multiway mergesort, which divides ${\displaystyle p}$ processors into ${\displaystyle r}$ groups of size ${\displaystyle p'}$. All processors sort locally first. Unlike single level multiway mergesort, these sequences are then partitioned into ${\displaystyle r}$ parts and assigned to the appropriate processor groups. These steps are repeated recursively in those groups. This reduces communication and especially avoids problems with many small messages. The hierarchical structure of the underlying real network can be used to define the processor groups (e.g. racks, clusters,...).[21]

Merge sort was one of the first sorting algorithms where optimal speed up was achieved, with Richard Cole using a clever subsampling algorithm to ensure O(1) merge.[23] Other sophisticated parallel sorting algorithms can achieve the same or better time bounds with a lower constant. For example, in 1991 David Powers described a parallelized quicksort (and a related radix sort) that can operate in O(log n) time on a CRCW parallel random-access machine (PRAM) with n processors by performing partitioning implicitly.[24] Powers further shows that a pipelined version of Batcher's Bitonic Mergesort at O((log n)²) time on a butterfly sorting network is in practice actually faster than his O(log n) sorts on a PRAM, and he provides detailed discussion of the hidden overheads in comparison, radix and parallel sorting.[25]