100 DAYS OF GPU Course preview
Day 36Module 4
in-technical-review

CUDA merge path: merging sorted arrays in parallel

Here is the whole of a sequential merge. Two sorted arrays in, one sorted array out.

while (k < m + n) {
    out[k++] = (j >= n || (i < m && a[i] <= b[j])) ? a[i++] : b[j++];
}

i and j leave one iteration and enter the next. Which array out[500] came from depends on all five hundred comparisons before it, so the loop is a chain of m + n steps and no step can start before the one above it finishes.

Read that as a fact about the answer and you conclude that merging does not parallelize. It is only a fact about the loop. A thread needs two numbers, and it can compute them with one binary search.

Each thread can then own a range of the output and find the matching part of each input. It needs no barrier or shared memory, and it reads nothing that another thread wrote.

The merge matrix, and the diagonal that cuts it

Lay a down the left of a grid and b across the top. A merge is a path from the top-left corner to the bottom-right: on each step you either take the next element of a and move down, or take the next element of b and move right. The path has m + n steps, one per output element, and it is the whole answer.

Take a = 1, 4, 6, 9 and b = 2, 3, 5, 7. The merged sequence is 1, 2, 3, 4, 5, 6, 7, 9, so the path goes down, right, right, down, right, down, right, down.

Now draw a second thing on the grid: the anti-diagonal i + j = k, the squares where exactly k elements have already been written. Every step of the path adds one to i + j, so the path crosses every diagonal exactly once. One crossing, one square, one pair (i, j) with i + j = k.

That square is a statement about the output: the first k elements of the merge are a[0..i) and b[0..j) and nothing else. The i is called the co-rank of k. Because j is k - i, one number names the split.

You do not have to walk the path to find where it crosses diagonal k. You can bisect. The answer lies between max(0, k - n) and min(k, m), and two tests say which side of the crossing a candidate i sits on. Neither reads more than two elements:

  • a[i - 1] > b[j] means the largest element you took from a beats one you left behind in b, so i is too big.
  • b[j - 1] >= a[i] means the largest element you took from b is at least one you left behind in a, so i is too small.

When neither fires you are standing on the path. The cost is a binary search over the shorter input, in global memory, and the answer for output index k owes nothing to the answer for k - 1. That last clause is the whole reason this is a GPU algorithm.

Diagram: the merge matrix and the diagonal that cuts it. A square grid with a = 1, 4, 6, 9 labelling the rows top to bottom and b = 2, 3, 5, 7 labelling the columns left to right, and a staircase from the top-left corner to the bottom-right that steps down when the merge takes from a and right when it takes from b. Band 1, the path alone: the staircase down, right, right, down, right, down, right, down, with the merged sequence written along it. Caption "One path, eight steps, and the path is the answer." Band 2, one diagonal: the anti-diagonal i + j = 3 drawn across the same grid, meeting the staircase at a single square. Caption "Diagonal 3 meets the path once, at i = 1." Band 3, five boundaries: k = 0, 2, 4, 6 and 8, with every crossing marked. The four path stretches between consecutive boundaries are four non-overlapping output slices. Caption "Five crossings delimit four slices." Alt text: "Five diagonal crossings at k equals 0, 2, 4, 6, and 8 delimit 4 non-overlapping merge slices; k equals 3 crosses at i equals 1."

The intuition the loop gives you, and why it is wrong

The loop suggests that you cannot know where output k comes from without performing the first k comparisons. That limit belongs to the loop, not the merge. The loop computes a monotone path, and you can find a diagonal crossing without walking the path from its start.

The second intuition is one this course spent day 11 arguing against, and here it comes back with the sign flipped. Day 11 says never give a thread its own contiguous chunk, because at any instant the 32 lanes of a warp are then a chunk apart and you have built the worst coalescing case on purpose. Merge path requires exactly that: a thread's slice has to be contiguous or its sequential merge has nothing to be sequential about.

Both statements hold. Day 11's rule describes the addresses one warp touches at one instant, and that access pattern still costs time here. You choose the amount of it by changing how many output elements each thread owns.

The program varies only that slice size. Its sweep measures the cost.

Measuring it against one thread

Full program in code/day36-merge/merge_path.cu.

Four checks keep the benchmark valid.

Every row is checked before it is timed. The host merges the same inputs with a plain two-pointer loop, and each kernel's whole output is compared against that before a stopwatch touches it. A row cannot get fast by writing less, which is the discipline day 11 exists to teach.

The comparison is exact, and that is a departure. The house comparator allows a relative tolerance because a GPU fuses a multiply and an add where the host compiler does not. A merge computes nothing, it only permutes, so there is nothing to allow for. On this input a tolerance would be worse than useless: neighbouring values differ by 1, and one part in a hundred thousand of the largest value is over 40, so a split off by twenty positions would pass.

The baseline is the same loop on the same card. mergeOneThread runs the sequential merge in a single thread of a single block, launched <<<1, 1>>>. It is a bad kernel on purpose. Comparing against a host merge instead would fold a PCIe crossing and a different clock into the answer, and day 9 is where that lesson was paid for.

One knob, five values. The slice size arrives as a runtime argument, so every row of the sweep runs the same instructions in the same order and only the number of threads and the length of each thread's run change.

The search is the lesson:

__device__ size_t coRank(size_t k, const float* __restrict__ a, size_t m,
                         const float* __restrict__ b, size_t n) {
    size_t i = (k < m) ? k : m;
    size_t j = k - i;
    size_t iLow = (k > n) ? (k - n) : 0;
    size_t jLow = (k > m) ? (k - m) : 0;

    while (true) {
        if (i > 0 && j < n && a[i - 1] > b[j]) {
            const size_t delta = (i - iLow + 1) >> 1;  // ceiling, not floor
            jLow = j;
            i -= delta;
            j += delta;
        } else if (j > 0 && i < m && b[j - 1] >= a[i]) {
            const size_t delta = (j - jLow + 1) >> 1;
            iLow = i;
            j -= delta;
            i += delta;
        } else {
            return i;
        }
    }
}

That function took two attempts, and the first one is worth showing you. The ceiling, not floor comment is not decoration. The first version read:

const size_t delta = (i - iLow) >> 1;  // wrong: floor, not ceiling

Trace it on four elements. Take a = 7, 8 and b = 1, 2, and ask for the co-rank of k = 2.

The search starts at i = 2, j = 0, iLow = 0, jLow = 0. Test one fires, delta is 1, and the state becomes i = 1, j = 1, jLow = 0.

Test one fires again because a[0] is 7 and b[1] is 2. Now i - iLow is 1, the floor of half of it is 0, and neither i nor j moves. jLow becomes 1, the state stops changing, and the loop does not return.

The ceiling fixes it. When test one fires, the answer's i is strictly below the current i, and iLow is never above the answer. Thus i - iLow is at least 1, and its half rounded up is also at least 1.

The same argument covers test two. Floor division yields 0 at the one width where the invariant still requires a step.

Note. The program never stages a tile of a or b in shared memory. A production merge does, because the reads inside a thread's slice are the chunked pattern. Leaving it out is what keeps the slice-size sweep readable, and it means the sweep below prices the naive kernel rather than the best this card can do.

Results

Re-verified without meaningful drift on a Tesla T4 with driver 580.173.02 and CUDA 13.0 (V13.0.88) on 2026-09-02. Every correctness gate and the four-elements-per-thread winner reproduced; timings stayed within two percent. The original CUDA 12.6 capture below and the CUDA 13 transcript both remain in evidence.

GPU: Tesla T4 (compute capability 7.5)
40 SMs, 1024 resident threads per SM, 40960 thread slots
a: 2097763 elements, b: 2097152, merged: 4194915, 32 MiB moved
co-rank bisects at most 2097152, so at most 22 steps per call

Part 1: the sequential merge, one thread on the device
kernel                          time (ms)         GB/s  vs 1 thread
----------------------------   ----------   ----------   ----------
mergeOneThread, <<<1, 1>>>        421.800         0.08         1.00

Part 2: mergeCoRank, 256 threads per block
 per thread      threads    time (ms)         GB/s  vs 1 thread
 ---------- ------------   ----------   ----------   ----------
          1      4194915        0.717        46.79       588.11
          4      1048729        0.346        97.10      1220.42
         16       262183        0.926        36.24       455.46
         64        65546        5.551         6.05        75.99
        256        16387        4.366         7.69        96.62

Every row above wrote the same 4194915 elements and the host compared
all of them against the reference before the row was timed, so no
row got fast by writing less. The GB/s column counts only the 32
MiB the merge has to move; the co-rank searches read on top of it.

One thread takes 421.800 ms. A thousand times faster is available, and the best configuration is not the most parallel one. Giving each thread a single element runs in 0.717 ms. Giving each thread four elements runs in 0.346 ms, twice as fast, at 1220 times the sequential version.

Past that point, 16 elements per thread costs 0.926 ms and 64 costs 5.551 ms. The co-rank search explains the results at small slice sizes. Each thread can take up to 22 binary-search steps before it merges any elements, so larger slices spread that search cost across more output.

The results at larger slice sizes need more care. At 64 elements per thread, the grid has 65,546 threads for 40,960 thread slots, or 1.6 waves. The second wave therefore leaves 60 percent of the card idle.

That does not explain the whole result. The 256-element row launches only 16,387 threads, fills 40 percent of the slots in one wave, and runs faster than the 64-element row at 4.366 ms. This run does not explain the change, so occupancy alone cannot predict the 64-element result.

The best measured slice size spreads out the search cost while keeping the device busy. On this card, that is 4 elements per thread. Other GPUs may favour a different value.

Run it yourself

Use Colab or a local CUDA GPU with compute capability 7.5 or later. The build command, from the repo's README:

nvcc -std=c++17 -O3 -arch=sm_75 -o merge_path merge_path.cu

What rules out an embed is mergeOneThread. At 421.800 ms per launch, its thirteen launches are about 5.5 seconds before the host reference merge and the five sweep rows run at all, so the 20 second run cap does not survive the whole program. Cut the array lengths by a factor of ten and it fits fine.

Exercise

Change >= to > in the second test of coRank and run the program. Every check still passes. Then make the program able to catch it: merge the same values a second time carrying, for each output position, which input array it came from, and compare that against the host merge.

Time: 30 to 45 minutes. Submit: the changed merge_path.cu, the first output position the new check reports, and one sentence on why the value check never fires.

Check: the harness runs both comparisons. The existing value comparison checks every element exactly. The provenance comparison prints the first output position where your kernel and the host merge disagree about which array supplied it, with both answers.

Failing the second check while passing the first is the intended outcome. Passing both means your input has no ties.

Hint 1

The program's a holds the even numbers and its b holds the odd ones, so no value in one equals any value in the other. Which of the two tests in coRank can only behave differently once that stops being true?

Hint 2

The merge inside a slice takes from a when a[i] <= b[j], so ties go to a. Write down what the split has to satisfy for that rule to be consistent with it, and compare it against the two tests. Then try a = 1, 5, 9 and b = 5, 5, 5 with three elements per thread, and look at output position 3.

Solution

With >= the co-rank of 3 is 2, so the second thread starts at a[2] and b[1]. With > it is 1, so it starts at a[1] and b[2], and a[1] has already been written by the first thread. The output becomes a[0], a[1], b[0], a[1], b[2], a[2]: one element written twice and b[1] never written at all.

The values are 1, 5, 5, 5, 5, 9 either way, and that is not luck. The two splits differ by trading a[i] for b[j - 1] at a point where those two are equal, so the multiset on each side is unchanged and no comparison of values can see it.

Pitfalls

Your co-rank loop never returns and the whole program hangs. A floor division where the search needs a ceiling leaves delta at 0 on the last step, so neither index moves. The symptom is a cudaDeviceSynchronize() that never comes back. Trace the four-element case above rather than adding a step counter: the bug happens at width 1, and a counter tells you it happened without telling you where.

You compute two co-ranks per thread when one will do. The end of a thread's slice is implied: start from a correct split and take exactly as many steps as the slice is long. A second search per thread doubles the part of the kernel that the sweep shows is expensive at small slice sizes.

You use a relative tolerance to compare the output. A merge permutes and computes nothing, so the comparison is exact. On sorted input, where neighbouring values are close together by construction, a tolerance is precisely the wrong tool: it hides index errors, which are the only errors this code can make.

You test with two arrays that share no values. Every tie-breaking bug in coRank is invisible then, including the one in the exercise above. Sorted inputs in the wild have ties, and a merge is usually the inner loop of a sort, where the same key appears in both halves constantly.

You conclude that more threads is better and stop at one element each. Each thread's binary search costs a walk over the shorter input's logarithm, in global memory, and one output element does not amortise it. The knob has two ends and both are bad; the sweep in this program exists to find the middle.

You expect a barrier to be needed somewhere. It is not, and that is the result worth carrying forward. No thread reads anything another thread wrote, so there is no __syncthreads() in either kernel and no shared memory either. A merge that does use shared memory uses it to stage the inputs for bandwidth, never to communicate.

Go deeper

Next

Day 37 keeps the theme and drops the guarantee: a sparse matrix in CSR gives every row a different amount of work, so which thread owns which piece of the input stops having a closed form. The piece this page leaves on the table is the tiled merge, which stages the two spans a block needs before any thread searches, and you already hold both halves of it from day 13 and day 16. The module is parallel patterns, and the GPU kernel engineer path says where it leads.