code/day36-merge/merge_path.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
//
// Day 36: merging two sorted arrays in parallel, with co-rank.
//
// Three ways to produce the same merged sequence:
//
// mergeCpu the host reference. Two pointers and one loop.
// mergeOneThread the same loop, on the device, in a single thread. This
// is the baseline the parallel version has to beat, and it
// runs on the same memory system, so the comparison does
// not smuggle in a PCIe crossing or a different clock.
// mergeCoRank one thread per slice of the output. Each thread runs a
// binary search on the merge diagonal to find where its
// slice starts in a and in b, then merges its slice
// sequentially. No barrier, no shared memory, nothing read
// that another thread wrote.
//
// mergeCoRank is measured at five slice sizes, because the slice size is the
// only knob and it moves two costs in opposite directions. A small slice
// means many threads and one binary search per output element. A large slice
// means few threads, and each thread then walks a long contiguous run, which
// is the per-thread chunk pattern day 11 measured at a seventh of this
// card's coalesced copy bandwidth.
//
// Every row produces the same output, and the host compares all of it
// against the reference before the row is timed, so a row cannot get fast by
// writing less. That is the discipline day 11 exists to teach.
//
// What this program does not do: it 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 that out is what makes the slice-size
// sweep readable, and the lesson says so rather than quietly implying that
// this is the fast version.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o merge_path merge_path.cu
// Run: ./merge_path
//
// Verified 2026-08-30 on a Tesla T4 (compute capability 7.5), driver
// 595.84, CUDA 12.6 (V12.6.85). Transcript: evidence/run-2026-08-30.txt
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, the way a Compiler Explorer
// embed is, so it carries its own verbatim copy. `err_` carries a trailing
// underscore so it cannot collide with a variable at the call site, and the
// do/while makes the macro one statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// On names. research/CUDA-CODE-STYLE.md reserves `i` for the global element
// index and asks a second index to be named for what it indexes. A merge has
// three arrays and needs three indices, and every published account of this
// algorithm calls them i into a, j into b and k into the output. Renaming
// them here would cost the reader every other source they will read on the
// subject, so the file keeps i, j and k and says so once, here.
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// a holds the even values 0, 2, 4, ... and b the odd values 1, 3, 5, ... so
// every value in the merge is distinct and a split bug cannot hide behind a
// repeated value. a is 611 elements longer than b, which does two jobs: the
// two co-rank bounds stop being mirror images of each other, and the last
// 611 outputs come from a alone, so the "one side is exhausted" path runs on
// every launch instead of only at the very end of the array.
constexpr size_t kElemsA = 2ull * 1024ull * 1024ull + 611ull;
constexpr size_t kElemsB = 2ull * 1024ull * 1024ull;
constexpr size_t kTotal = kElemsA + kElemsB;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kElemsA != kElemsB,
"the inputs must differ in length, or the two co-rank bounds "
"max(0, k - n) and max(0, k - m) are mirror images and a bug "
"in one of them hides behind the other");
static_assert(kTotal % 2 == 1,
"an odd total means every slice size in the sweep leaves a "
"short last slice, so the tail path runs on every launch");
static_assert(2 * (kElemsA - 1) < (1ull << 24),
"every value has to be an integer float holds exactly, or two "
"different elements could compare equal and hide a bug");
// Given an output index k, return the i with i + j = k such that a[0..i) and
// b[0..j) are exactly the k smallest elements of the merged sequence. One
// number names the split, because j is k - i.
//
// The search bisects i. Two invariants hold at the top of every iteration:
// iLow is never above the answer's i, and jLow is never above the answer's
// j. They are what make the loop finish.
//
// Test one fires when a[i - 1] > b[j], which says the answer's i is
// strictly below i. With iLow at or below the answer, that gives
// i - iLow >= 1. Test two fires when b[j - 1] >= a[i], which says the
// answer's j is strictly below j, giving j - jLow >= 1 the same way.
//
// delta is therefore at least 1 on every step, and only because the + 1
// before the shift makes it a ceiling. A floor division returns 0 when the
// interval is 1 wide, neither i nor j moves, and the loop never ends.
// README.md traces the four-element input where that happens.
//
// The two tests are not mirror images and the asymmetry is deliberate: `>`
// on one side and `>=` on the other put a tie in a ahead of the same value
// in b, which is what makes this merge stable. With keys alone the two
// choices give the same output values, so nothing here can see the
// difference; the exercise is where it becomes visible.
//
// Cost: the interval starts no wider than min(m, n) and at least halves each
// step, so the loop is logarithmic in the shorter input and every step is a
// dependent global load.
// snippet: corank
__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;
}
}
}
// end snippet
// The whole merge in one thread, so the parallel version has something on
// the same card to be a multiple of.
//
// One thread: walks a and b with two pointers and writes every one of the
// m + n outputs.
//
// One warp: there is no warp. Lane 0 of block 0 does all of it and the other
// 31 lanes are masked off from the first instruction, so every load fetches
// a 32-byte sector to deliver the 4 bytes lane 0 asked for. The two input
// streams are sequential, so most of those sectors are already in cache;
// this row is slow because it is one thread, not because it is uncoalesced.
//
// Launch assumption: <<<1, 1>>>. The guard makes any other launch write
// nothing rather than write the whole array once per thread.
// snippet: one-thread
__global__ void mergeOneThread(const float* __restrict__ a, size_t m,
const float* __restrict__ b, size_t n,
float* __restrict__ out) {
if (blockIdx.x == 0 && threadIdx.x == 0) {
size_t i = 0;
size_t j = 0;
for (size_t k = 0; k < m + n; ++k) {
if (j >= n || (i < m && a[i] <= b[j])) {
out[k] = a[i];
++i;
} else {
out[k] = b[j];
++j;
}
}
}
}
// end snippet
// One thread owns out[kStart, kEnd), a slice of perThread outputs.
//
// One thread: one co-rank search to find where its slice begins in a and in
// b, then the sequential merge above for exactly kEnd - kStart steps. It
// needs no second search: starting from a correct split and taking that many
// steps lands on the same elements the whole-array merge would have written
// there.
//
// One warp: the write is coalesced when perThread is 1, because lane L
// writes out[t]. Past that the 32 lanes are perThread elements apart at
// every step of the loop, which is the chunk pattern from day 11. The reads
// inside coRank are the opposite case and better than they look: adjacent
// output indices have almost the same co-rank, so the 32 lanes of a warp
// bisect nearly the same interval and land on nearly the same addresses.
//
// Launch assumption: gridDim.x * blockDim.x * perThread >= m + n. Nothing
// here reads what another thread wrote, so there is no barrier and no shared
// memory, which is the point of the whole page.
// snippet: merge-corank
__global__ void mergeCoRank(const float* __restrict__ a, size_t m,
const float* __restrict__ b, size_t n,
float* __restrict__ out, int perThread) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t total = m + n;
const size_t kStart = t * static_cast<size_t>(perThread);
if (kStart < total) {
const size_t kSpan = static_cast<size_t>(perThread);
const size_t kEnd = (kStart + kSpan < total) ? (kStart + kSpan) : total;
size_t i = coRank(kStart, a, m, b, n);
size_t j = kStart - i;
for (size_t k = kStart; k < kEnd; ++k) {
if (j >= n || (i < m && a[i] <= b[j])) {
out[k] = a[i];
++i;
} else {
out[k] = b[j];
++j;
}
}
}
}
// end snippet
// CPU reference. Written for obvious correctness, not speed: plain loop, no
// OpenMP, no intrinsics. It never allocates; the caller owns every buffer.
// `a[i] <= b[j]` rather than `<` is what makes it stable, and it has to
// match the rule both kernels use or the comparison is between two different
// answers.
static void mergeCpu(const float* a, size_t m, const float* b, size_t n,
float* out) {
size_t i = 0;
size_t j = 0;
for (size_t k = 0; k < m + n; ++k) {
if (j >= n || (i < m && a[i] <= b[j])) {
out[k] = a[i];
++i;
} else {
out[k] = b[j];
++j;
}
}
}
// Returns the first index where got and want differ at all, or n if they
// agree everywhere. Exact, not within a tolerance, and that is a deliberate
// departure from the house comparator: a merge permutes its input and
// computes nothing, so there is no rounding to allow for. A 1e-5 relative
// tolerance would be actively wrong on this input, because neighbouring
// values differ by 1 while 1e-5 of the largest value is over 40, so a split
// off by twenty positions would pass.
static size_t firstDifference(const float* got, const float* want, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
return i;
}
}
return n;
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
// A host clock around a launch measures the launch, because launches are
// asynchronous. Day 9 takes that apart. Copy this helper verbatim.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
// Every kernel here reads kTotal floats and writes kTotal floats, so the
// figure comes from one constant rather than from each kernel, where the two
// could drift apart. It counts the traffic the merge cannot avoid and
// nothing else: the co-rank searches read on top of it, which is part of
// what the sweep is measuring.
static double mergeGBs(float ms) {
const double bytes = 2.0 * static_cast<double>(kTotal) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Copies the output back and compares all of it against the reference.
// Returns 1 on a mismatch and 0 on agreement, so the caller adds it to a
// failure count and carries on to the next row instead of returning early
// with device memory still allocated.
static int checkOutput(const char* label, const float* d_out, float* h_out,
const float* want) {
CUDA_CHECK(cudaMemcpy(h_out, d_out, kTotal * sizeof(float),
cudaMemcpyDeviceToHost));
const size_t bad = firstDifference(h_out, want, kTotal);
if (bad != kTotal) {
std::fprintf(stderr, "%s wrong at %zu: got %.9g, want %.9g\n", label,
bad, h_out[bad], want[bad]);
return 1;
}
return 0;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("%d SMs, %d resident threads per SM, %d thread slots\n",
prop.multiProcessorCount, prop.maxThreadsPerMultiProcessor,
prop.multiProcessorCount * prop.maxThreadsPerMultiProcessor);
const double movedMiB =
2.0 * static_cast<double>(kTotal) * sizeof(float) / (1024.0 * 1024.0);
std::printf("a: %zu elements, b: %zu, merged: %zu, %.0f MiB moved\n",
kElemsA, kElemsB, kTotal, movedMiB);
// Arithmetic, not a measurement. The co-rank interval starts no wider
// than min(m, n) and at least halves on every step, so this bounds the
// loop without anyone timing it.
const size_t coRankSpan = (kElemsA < kElemsB) ? kElemsA : kElemsB;
size_t span = coRankSpan;
int maxSteps = 0;
while (span > 0) {
span >>= 1;
++maxSteps;
}
std::printf("co-rank bisects at most %zu, so at most %d steps per call\n",
coRankSpan, maxSteps);
// Every failure below records itself and falls through to the one
// cleanup block at the bottom, so no path can return with device memory
// still allocated.
int failures = 0;
float* d_a = nullptr;
float* d_b = nullptr;
float* d_out = nullptr;
const size_t bytesA = kElemsA * sizeof(float);
const size_t bytesB = kElemsB * sizeof(float);
const size_t bytesOut = kTotal * sizeof(float);
CUDA_CHECK(cudaMalloc(&d_a, bytesA));
CUDA_CHECK(cudaMalloc(&d_b, bytesB));
CUDA_CHECK(cudaMalloc(&d_out, bytesOut));
std::vector<float> h_a(kElemsA);
std::vector<float> h_b(kElemsB);
std::vector<float> h_out(kTotal);
std::vector<float> h_want(kTotal);
for (size_t i = 0; i < kElemsA; ++i) {
h_a[i] = static_cast<float>(2 * i);
}
for (size_t i = 0; i < kElemsB; ++i) {
h_b[i] = static_cast<float>(2 * i + 1);
}
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytesA, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytesB, cudaMemcpyHostToDevice));
mergeCpu(h_a.data(), kElemsA, h_b.data(), kElemsB, h_want.data());
// The big input keeps every merged value distinct, so it can never
// exercise the a[i - 1] == b[j] case the co-rank comments spend five
// lines on. This small case exists to run that branch: duplicated
// values across both inputs, one element per thread, checked against
// the host merge that uses the matching <= rule. Without it the
// stability the `>` versus `>=` asymmetry buys is a claim no run of
// this program tests.
{
const float tinyA[] = {1.0f, 1.0f, 2.0f, 3.0f};
const float tinyB[] = {1.0f, 2.0f, 2.0f, 4.0f};
const size_t tinyM = sizeof(tinyA) / sizeof(tinyA[0]);
const size_t tinyN = sizeof(tinyB) / sizeof(tinyB[0]);
float tinyWant[tinyM + tinyN];
float tinyGot[tinyM + tinyN];
mergeCpu(tinyA, tinyM, tinyB, tinyN, tinyWant);
CUDA_CHECK(
cudaMemcpy(d_a, tinyA, sizeof(tinyA), cudaMemcpyHostToDevice));
CUDA_CHECK(
cudaMemcpy(d_b, tinyB, sizeof(tinyB), cudaMemcpyHostToDevice));
mergeCoRank<<<1, kThreadsPerBlock>>>(d_a, tinyM, d_b, tinyN, d_out, 1);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(tinyGot, d_out, sizeof(tinyWant),
cudaMemcpyDeviceToHost));
const size_t bad = firstDifference(tinyGot, tinyWant, tinyM + tinyN);
if (bad != tinyM + tinyN) {
std::fprintf(stderr,
"equal-element merge wrong at %zu: got %.9g, want "
"%.9g\n",
bad, tinyGot[bad], tinyWant[bad]);
++failures;
}
// Put the big input back; everything below assumes it is resident.
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytesA, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytesB, cudaMemcpyHostToDevice));
}
std::printf("\nPart 1: the sequential merge, one thread on the device\n");
std::printf("%-28s %12s %12s %12s\n", "kernel", "time (ms)", "GB/s",
"vs 1 thread");
std::printf("%-28s %12s %12s %12s\n", "----------------------------",
"----------", "----------", "----------");
// One untimed launch first, checked, so a wrong answer is reported
// before a number that came from it reaches the table.
mergeOneThread<<<1, 1>>>(d_a, kElemsA, d_b, kElemsB, d_out);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures +=
checkOutput("mergeOneThread", d_out, h_out.data(), h_want.data());
const float oneThreadMs = timeKernel(
[&] { mergeOneThread<<<1, 1>>>(d_a, kElemsA, d_b, kElemsB, d_out); });
int rows = 1;
std::printf("%-28s %12.3f %12.2f %12.2f\n", "mergeOneThread, <<<1, 1>>>",
oneThreadMs, mergeGBs(oneThreadMs), 1.0);
std::printf("\nPart 2: mergeCoRank, %d threads per block\n",
kThreadsPerBlock);
std::printf("%11s %12s %12s %12s %12s\n", "per thread", "threads",
"time (ms)", "GB/s", "vs 1 thread");
std::printf("%11s %12s %12s %12s %12s\n", "----------", "------------",
"----------", "----------", "----------");
const int slices[] = {1, 4, 16, 64, 256};
const int kNumSlices = sizeof(slices) / sizeof(slices[0]);
for (int s = 0; s < kNumSlices; ++s) {
const int perThread = slices[s];
const size_t threads = (kTotal + static_cast<size_t>(perThread) - 1) /
static_cast<size_t>(perThread);
const int blocks = static_cast<int>((threads + kThreadsPerBlock - 1) /
kThreadsPerBlock);
mergeCoRank<<<blocks, kThreadsPerBlock>>>(d_a, kElemsA, d_b, kElemsB,
d_out, perThread);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
char label[64];
std::snprintf(label, sizeof(label), "mergeCoRank, %d per thread",
perThread);
failures += checkOutput(label, d_out, h_out.data(), h_want.data());
const float ms = timeKernel([&] {
mergeCoRank<<<blocks, kThreadsPerBlock>>>(
d_a, kElemsA, d_b, kElemsB, d_out, perThread);
});
++rows;
std::printf("%11d %12zu %12.3f %12.2f %12.2f\n", perThread, threads, ms,
mergeGBs(ms), oneThreadMs / ms);
}
std::printf(
"\nEvery row above wrote the same %zu elements and the host compared\n"
"all of them against the reference before the row was timed, so no\n"
"row got fast by writing less. The GB/s column counts only the %.0f\n"
"MiB the merge has to move; the co-rank searches read on top of it.\n",
kTotal, movedMiB);
// The row count is checked so the lesson's tables cannot drift from what
// the program prints. A real branch, not an assert: CI builds Release,
// Release defines NDEBUG, and NDEBUG deletes assert(), so the check would
// be missing from exactly the build that matters.
const int kExpectedRows = kNumSlices + 1;
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's tables and "
"this program disagree\n",
rows, kExpectedRows);
++failures;
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_out));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}