code/day33-compaction/compaction.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 33: stream compaction and partition, by predicate, scan and scatter.
//
// The pattern this file exists to show is three steps and nothing else:
//
// predicate rank[i] = 1 when in[i] survives, 0 when it does not
// scan rank[i] becomes the number of survivors before i, which is
// the slot in the output that in[i] belongs in
// scatter if in[i] survives, out[rank[i]] = in[i]
//
// No thread waits for another thread, nothing is locked, and the number of
// survivors falls out of the top of the scan without anybody counting it.
//
// Two problems, one machine. Removing the zeros is the first. The even/odd
// partition scans the same way and then uses one piece of arithmetic:
// i - rank[i] is the number of odd elements before i, so a single scan
// addresses both halves and the second counter never has to exist.
//
// compactAtomic is the comparison, and it is what most people write first:
// one atomicAdd on a counter the whole grid shares, one write to whatever
// slot that atomic handed back. Two launches instead of seven. It computes a
// different function, because the order in which blocks reach the counter is
// the order the scheduler happens to run them, so this program checks the two
// outputs differently: both must hold the same multiset of survivors, only
// the scanned one must be in the input's order.
//
// The scan is day 32's, restated here so this file stands alone: three
// launches of a tile-local Hillis-Steele scan, then two launches that add
// each tile's carry back in. A tile is 256 elements and three levels reach
// 256^3 = 16,777,216 elements, which is the ceiling kElems sits under and
// which the third static_assert below refuses to let anyone cross quietly.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o compaction compaction.cu
// Run: ./compaction
//
// 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)
// 8 Mi ints minus 611. The 611 is day 5's number, here for the same reason:
// the length is not a multiple of the tile size, so the last tile of every
// scan is a partial one and a tail bug shows up as a wrong count rather than
// as nothing. 32 MiB per buffer is eight times this card's L2, so the reads
// come from DRAM and not from a cache that happens to be holding the array.
constexpr size_t kElems = 8ull * 1024ull * 1024ull - 611ull;
constexpr int kThreadsPerBlock = 256;
// The fixed grid the four grid-stride kernels use. Every one of them walks
// the whole array with a stride loop, so this number sets occupancy and not
// coverage.
constexpr int kBlocks = 1024;
// One scan tile, and the block size the scan kernel must be launched with.
// It is separate from kThreadsPerBlock because the carry kernel divides by it
// once per element and has to agree with the scan kernel regardless of what
// grid the carry kernel itself runs on.
constexpr int kTileElems = 256;
constexpr size_t kTileSpan = static_cast<size_t>(kTileElems);
// Three levels of tiles. Level 0 scans the ranks, level 1 scans level 0's
// tile totals, level 2 scans level 1's, and level 2 has to fit in one tile
// because there is no fourth launch in this file.
constexpr size_t kTiles0 = (kElems + kTileSpan - 1) / kTileSpan;
constexpr size_t kTiles1 = (kTiles0 + kTileSpan - 1) / kTileSpan;
constexpr size_t kTiles2 = (kTiles1 + kTileSpan - 1) / kTileSpan;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// Three inputs, differing only in how many elements survive. The scanned
// compaction does the same work whatever the answer is; the atomic one does
// one atomic per survivor, so its cost is supposed to move down this list.
constexpr int kNumInputs = 3;
constexpr int kMidInput = 1;
constexpr unsigned int kKeepThreshold[kNumInputs] = {26u, 128u, 230u};
// Surviving values are drawn from 1 to 255, so a 256-entry tally is an exact
// multiset test and zero is never a legitimate survivor.
constexpr int kValueRange = 256;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kTileElems & (kTileElems - 1)) == 0,
"addTileCarry divides by kTileElems once per element, and a "
"power of two makes that a shift");
static_assert(kTiles2 == 1,
"this file launches exactly three levels of tile scan, so "
"kElems must be at most kTileElems cubed; raise the level count "
"before raising kElems");
static_assert(kElems <= 2147483647ull,
"ranks, tile totals and the survivor count are all int, so the "
"element count has to fit in one");
// Fills `out` with values in 1 to 255, zeroing every element whose keep draw
// lands at or above keepThreshold. Two draws per element, so the survival
// rate and the surviving values are independent and the multiset test below
// is not measuring the same bits twice.
//
// Fixed constants, no clock and no seed argument, so two people comparing
// numbers are comparing the same bytes.
static void fillInput(unsigned int keepThreshold, int* out, size_t n) {
unsigned int state = 12345u;
for (size_t i = 0; i < n; ++i) {
state = state * 1664525u + 1013904223u;
// The high byte, never the low one. The low bits of a linear
// congruential generator cycle with a very short period, which would
// put a pattern into the survival mask by accident.
const unsigned int keep = state >> 24;
state = state * 1664525u + 1013904223u;
const unsigned int value = 1u + (state >> 24) % 255u;
out[i] = (keep < keepThreshold) ? static_cast<int>(value) : 0;
}
}
// Stable compaction on the host: the nonzero elements, in the order they
// appear. Returns how many survived. Written for obvious correctness, not
// speed, and it never allocates; the caller owns both buffers.
static size_t compactCpu(const int* in, int* out, size_t n) {
size_t m = 0;
for (size_t i = 0; i < n; ++i) {
if (in[i] != 0) {
out[m] = in[i];
++m;
}
}
return m;
}
// Stable two-way partition on the host: every even element in input order,
// then every odd one in input order. Returns how many were even. Two passes
// here, because the host has no reason to be clever and one obvious loop per
// half is easier to check than one loop with two cursors.
static size_t partitionCpu(const int* in, int* out, size_t n) {
size_t evens = 0;
for (size_t i = 0; i < n; ++i) {
if ((in[i] & 1) == 0) {
out[evens] = in[i];
++evens;
}
}
size_t odds = evens;
for (size_t i = 0; i < n; ++i) {
if ((in[i] & 1) != 0) {
out[odds] = in[i];
++odds;
}
}
return evens;
}
// Exact, element by element. These are integers moved from one address to
// another, so there is nothing to round and no tolerance to set. A tolerance
// here would hide the off-by-one that this whole program is about.
static int checkExact(const char* name, const int* got, const int* want,
size_t m) {
for (size_t k = 0; k < m; ++k) {
if (got[k] != want[k]) {
std::fprintf(stderr, "%s wrong at %zu: got %d, want %d\n", name, k,
got[k], want[k]);
return EXIT_FAILURE;
}
}
return EXIT_SUCCESS;
}
// The same values in some order, which is the most that can be asked of a
// kernel whose output order the scheduler decides. One tally, incremented
// from `got` and decremented from `want`, so every entry has to come back to
// zero.
static int checkSameMultiset(const char* name, const int* got, const int* want,
size_t m) {
std::vector<long long> tally(static_cast<size_t>(kValueRange), 0);
for (size_t k = 0; k < m; ++k) {
if (got[k] < 0 || got[k] >= kValueRange) {
std::fprintf(stderr, "%s wrote %d at %zu, outside [0, %d)\n", name,
got[k], k, kValueRange);
return EXIT_FAILURE;
}
// The reference gets the same test. It is in range today by
// construction, but this helper indexes the tally with both sides,
// and an unchecked want[] would turn a generator change into an
// out-of-bounds write that looks like heap corruption.
if (want[k] < 0 || want[k] >= kValueRange) {
std::fprintf(stderr,
"%s reference holds %d at %zu, outside [0, %d)\n",
name, want[k], k, kValueRange);
return EXIT_FAILURE;
}
tally[static_cast<size_t>(got[k])] += 1;
tally[static_cast<size_t>(want[k])] -= 1;
}
for (int v = 0; v < kValueRange; ++v) {
if (tally[static_cast<size_t>(v)] != 0) {
std::fprintf(stderr,
"%s holds the wrong multiset: value %d is off by "
"%lld\n",
name, v, tally[static_cast<size_t>(v)]);
return EXIT_FAILURE;
}
}
return EXIT_SUCCESS;
}
// How many positions of the atomic output hold something other than the
// stable answer. Reported, never gated: a run where the two agreed would be
// luck rather than correctness, and a check that fails on luck is worse than
// no check at all.
static size_t countOutOfOrder(const int* got, const int* want, size_t m) {
size_t differ = 0;
for (size_t k = 0; k < m; ++k) {
if (got[k] != want[k]) {
++differ;
}
}
return differ;
}
// rank[i] = 1 when in[i] survives the predicate, 0 when it does not.
//
// One thread handles many elements. Memory: consecutive lanes take
// consecutive elements on every step of the stride loop, so one warp reads
// 128 contiguous bytes and writes 128 contiguous bytes.
//
// The predicate is the only thing that changes between compaction and
// partition. Everything downstream reads an array of ones and zeros and does
// not know or care what produced it, which is the reason this pattern is
// reusable at all.
//
// Launch assumption: any grid. The loop condition is the bounds check.
__global__ void markNonZero(const int* __restrict__ in, int* __restrict__ rank,
size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
rank[i] = (in[i] != 0) ? 1 : 0;
}
}
// rank[i] = 1 when in[i] is even. The same kernel as markNonZero with one
// character changed, which is the point: the machinery below is the same for
// every predicate anyone will ever write.
//
// Launch assumption: any grid. The loop condition is the bounds check.
__global__ void markEven(const int* __restrict__ in, int* __restrict__ rank,
size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
rank[i] = ((in[i] & 1) == 0) ? 1 : 0;
}
}
// Exclusive scan of one tile of kTileElems values, one value per thread.
//
// Writes each thread's exclusive prefix to out[i] and the whole tile's total
// to sums[blockIdx.x]. Those totals are what the next level scans, which is
// how three launches of this one kernel cover the whole array.
//
// Memory: consecutive lanes take consecutive elements, so the load and the
// store each cover 128 contiguous bytes per warp. Each of the two shared
// reads per step walks consecutive words across the warp, so each one covers
// all 32 banks once and neither conflicts.
//
// in and out may be the same pointer, and every caller here does exactly
// that, which is why nothing in this signature is __restrict__ while every
// other kernel in the file marks all of its pointers. __restrict__ promises
// the memory is reached through that pointer only, and passing one array as
// two restricted parameters is a promise the caller has already broken.
// Scanning in place is safe for a different reason: a thread reads in[i] into
// a register before the first barrier and writes out[i] after the last one,
// and no other thread touches index i.
//
// Hillis-Steele, double buffered. The single-buffer form reads and writes the
// same shared word inside one statement, and there is nowhere to put the
// barrier that would make that safe.
//
// Launch assumption: exactly kTileElems threads per block. The shared array
// is sized from that constant, the step loop counts up to it, and the tile
// total is written by the thread holding the last element of the tile.
__global__ void scanTilesExclusive(const int* in, int* out, int* sums,
size_t n) {
__shared__ int buf[2][kTileElems];
const unsigned int tid = threadIdx.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;
// The guard covers the load, not the barrier. Threads past the end of a
// partial tile contribute a zero and still reach every __syncthreads(),
// because a return above a barrier hangs the block.
const int value = (i < n) ? in[i] : 0;
int cur = 0;
buf[cur][tid] = value;
__syncthreads();
for (unsigned int offset = 1; offset < kTileElems; offset <<= 1) {
const int next = cur ^ 1;
const int left = (tid >= offset) ? buf[cur][tid - offset] : 0;
buf[next][tid] = buf[cur][tid] + left;
__syncthreads();
cur = next;
}
// buf[cur][tid] is now the inclusive prefix. Subtracting this thread's
// own value turns it into the exclusive one, with no second pass and no
// shift of the whole tile.
if (i < n) {
out[i] = buf[cur][tid] - value;
}
if (tid == kTileElems - 1) {
sums[blockIdx.x] = buf[cur][tid];
}
}
// Adds each tile's carry-in to every element of that tile.
//
// carries[t] is the exclusive scan of the tile totals, so adding it turns a
// tile-local prefix into a whole-array one. The tile index is i / kTileElems
// and not anything derived from blockDim, because this kernel runs on a grid
// of its own choosing while the tiles were fixed by the scan.
//
// Memory: one coalesced read-modify-write of out. Every lane of a warp reads
// the same carries entry, which the hardware serves as one broadcast.
//
// Launch assumption: any grid. The loop condition is the bounds check.
__global__ void addTileCarry(int* __restrict__ out,
const int* __restrict__ carries, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] += carries[i / kTileSpan];
}
}
// Writes every surviving element to the slot the scan gave it.
//
// Memory: the reads of in and rank are coalesced. The write is not the random
// scatter the word suggests. Destinations rise with i, so one warp's
// survivors land in a run of at most 32 consecutive slots, and a warp that
// keeps everything writes 128 contiguous bytes exactly like a copy.
//
// The predicate is recomputed here instead of being read back from a flags
// array. That saves a pass over n integers and it cannot disagree with the
// scan, because both answers come from the same in[i].
//
// Launch assumption: any grid, and rank[] already scanned. The loop condition
// is the bounds check.
// snippet: scatter
__global__ void scatterCompact(const int* __restrict__ in,
const int* __restrict__ rank,
int* __restrict__ out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
const int value = in[i];
if (value != 0) {
out[rank[i]] = value;
}
}
}
// end snippet: scatter
// Splits in[] into the even elements then the odd ones, keeping the input
// order inside each half.
//
// rank[i] is the number of even elements before i, so i - rank[i] is the
// number of odd ones before i. One scan therefore addresses both halves and
// the second counter this pattern appears to need never exists.
//
// totalEven is read from device memory rather than passed by value, because
// it is the top of the same scan. Copying it to the host to fill in a launch
// argument would put a synchronising memcpy in the middle of the pipeline.
//
// Memory: the same monotone scatter as scatterCompact, except that two runs
// advance at once, so a warp writes into at most two short runs instead of
// one.
//
// Launch assumption: any grid, rank[] already scanned, and totalEven pointing
// at the total that same scan produced. A total from a different scan makes
// the two halves overlap.
// snippet: partition
__global__ void scatterPartition(const int* __restrict__ in,
const int* __restrict__ rank,
int* __restrict__ out, size_t n,
const int* __restrict__ totalEven) {
const size_t base = static_cast<size_t>(*totalEven);
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
const int value = in[i];
const size_t evensBefore = static_cast<size_t>(rank[i]);
if ((value & 1) == 0) {
out[evensBefore] = value;
} else {
out[base + (i - evensBefore)] = value;
}
}
}
// end snippet: partition
// Compaction with a counter instead of a scan. One pass over the data, and a
// different function of it.
//
// Memory: the read is coalesced. Every surviving lane then issues an atomic
// on the single 4-byte address the whole grid shares and writes one element
// wherever that atomic sent it, so the write pattern is decided at run time
// by the order the atomics resolve.
//
// This does not preserve order and it is not a bug that it does not. Blocks
// reach the counter in whatever order the scheduler runs them, so element
// 8,000,000 can land ahead of element 3. The survivors and their count are
// right; only their arrangement is not, which is why the program checks those
// two things separately.
//
// Launch assumption: any grid, and *count zeroed before the launch*. Running
// it twice against a counter nobody reset walks off the end of out.
// snippet: atomic
__global__ void compactAtomic(const int* __restrict__ in, int* __restrict__ out,
unsigned int* __restrict__ count, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
const int value = in[i];
if (value != 0) {
const unsigned int slot = atomicAdd(count, 1u);
out[slot] = value;
}
}
}
// end snippet: atomic
// Sets the counter back to zero. One thread, and it has to sit inside
// compactAtomic's timed region because the kernel cannot be launched a second
// time without it. The table on the page counts it as a launch and says so.
__global__ void zeroCount(unsigned int* __restrict__ count) {
*count = 0u;
}
// The floor. Reads n ints and writes n ints with no predicate, no scan and no
// scatter, from the same grid over the same buffers in the same process, so
// the rows below have something measured to be a fraction of.
//
// Launch assumption: any grid. The loop condition is the bounds check.
__global__ void copyInts(const int* __restrict__ in, int* __restrict__ out,
size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = in[i];
}
}
// Turns rank[] from an array of ones and zeros into the exclusive prefix sum
// of itself, in place, and leaves the total in sums2[0]. Five launches and no
// host round trip, so the survivor count never leaves the device and nothing
// downstream has to wait for a copy before it can run.
//
// No error check in here. Every caller either wraps this in timeKernel, which
// checks after its own loops, or follows it with cudaGetLastError() and
// cudaDeviceSynchronize().
static void launchScan(int* d_rank, int* d_sums0, int* d_sums1, int* d_sums2) {
scanTilesExclusive<<<static_cast<int>(kTiles0), kTileElems>>>(
d_rank, d_rank, d_sums0, kElems);
scanTilesExclusive<<<static_cast<int>(kTiles1), kTileElems>>>(
d_sums0, d_sums0, d_sums1, kTiles0);
scanTilesExclusive<<<static_cast<int>(kTiles2), kTileElems>>>(
d_sums1, d_sums1, d_sums2, kTiles1);
addTileCarry<<<kBlocks, kThreadsPerBlock>>>(d_sums0, d_sums1, kTiles0);
addTileCarry<<<kBlocks, kThreadsPerBlock>>>(d_rank, d_sums0, kElems);
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
// A host clock around a launch measures the launch and not the kernel,
// because launches are asynchronous. Day 9 takes that apart.
//
// This is the one template and the one lambda allowed in module 1 to 3 code,
// carried forward here unchanged. Copy it 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;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t bytes = kElems * sizeof(int);
const double mib = static_cast<double>(bytes) / (1024.0 * 1024.0);
// The banner prints before anything can fail, like every other file in
// this module, because a failing run is exactly the one where the card,
// grid and sizes need to be on the record.
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("%d SMs, %zu bytes of shared memory per block\n",
prop.multiProcessorCount, prop.sharedMemPerBlock);
std::printf("Input: %zu ints (%.2f MiB per buffer)\n", kElems, mib);
std::printf("Grid: %d blocks x %d threads for the stride kernels\n",
kBlocks, kThreadsPerBlock);
std::printf("Scan: %d elements per tile, %zu then %zu then %zu tiles\n",
kTileElems, kTiles0, kTiles1, kTiles2);
std::printf(
"Mean of %d runs after %d warm-ups. Host copies are "
"outside every timing.\n\n",
kTimedRuns, kWarmupRuns);
std::vector<int> h_in(kElems);
std::vector<int> h_out(kElems);
std::vector<int> h_want(kElems);
std::vector<int> h_atomic(kElems);
int* d_in = nullptr;
int* d_rank = nullptr;
int* d_out = nullptr;
int* d_atomic = nullptr;
int* d_sums0 = nullptr;
int* d_sums1 = nullptr;
int* d_sums2 = nullptr;
unsigned int* d_count = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_rank, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMalloc(&d_atomic, bytes));
CUDA_CHECK(cudaMalloc(&d_sums0, kTiles0 * sizeof(int)));
CUDA_CHECK(cudaMalloc(&d_sums1, kTiles1 * sizeof(int)));
CUDA_CHECK(cudaMalloc(&d_sums2, kTiles2 * sizeof(int)));
CUDA_CHECK(cudaMalloc(&d_count, sizeof(unsigned int)));
int status = EXIT_SUCCESS;
size_t survivors[kNumInputs] = {0, 0, 0};
double keptPercent[kNumInputs] = {0.0, 0.0, 0.0};
size_t outOfOrder[kNumInputs] = {0, 0, 0};
float compactMs[kNumInputs] = {0.0f, 0.0f, 0.0f};
float atomicMs[kNumInputs] = {0.0f, 0.0f, 0.0f};
float copyMs = 0.0f;
float markMs = 0.0f;
float markScanMs = 0.0f;
float scatterMs = 0.0f;
float partitionMs = 0.0f;
size_t evenCount = 0;
for (int which = 0; which < kNumInputs; ++which) {
fillInput(kKeepThreshold[which], h_in.data(), kElems);
const size_t mWant = compactCpu(h_in.data(), h_want.data(), kElems);
survivors[which] = mWant;
keptPercent[which] =
100.0 * static_cast<double>(mWant) / static_cast<double>(kElems);
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
// Predicate, scan, scatter. Seven launches, nothing synchronous
// between them, and the count stays on the card.
markNonZero<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, kElems);
launchScan(d_rank, d_sums0, d_sums1, d_sums2);
scatterCompact<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, d_out,
kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
// Gate 1: the top of the scan is the survivor count. A scan that is
// off by one anywhere shows up here before any element is compared.
int scannedCount = 0;
CUDA_CHECK(cudaMemcpy(&scannedCount, d_sums2, sizeof(int),
cudaMemcpyDeviceToHost));
if (static_cast<size_t>(scannedCount) != mWant) {
std::fprintf(stderr,
"the scan counted %d survivors, the host counted "
"%zu\n",
scannedCount, mWant);
status = EXIT_FAILURE;
break;
}
// Gate 2: the compacted array, exactly, in order. A real branch
// returning EXIT_FAILURE and not an assert(), because CI builds
// Release, Release defines NDEBUG, and NDEBUG deletes assert().
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, mWant * sizeof(int),
cudaMemcpyDeviceToHost));
if (checkExact("compactScanned", h_out.data(), h_want.data(), mWant) !=
EXIT_SUCCESS) {
status = EXIT_FAILURE;
break;
}
zeroCount<<<1, 1>>>(d_count);
compactAtomic<<<kBlocks, kThreadsPerBlock>>>(d_in, d_atomic, d_count,
kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
// Gate 3: the counter has to agree with the host, which is the check
// that catches a lost update.
unsigned int atomicCount = 0u;
CUDA_CHECK(cudaMemcpy(&atomicCount, d_count, sizeof(unsigned int),
cudaMemcpyDeviceToHost));
if (static_cast<size_t>(atomicCount) != mWant) {
std::fprintf(stderr,
"the atomic counter reached %u, the host counted "
"%zu\n",
atomicCount, mWant);
status = EXIT_FAILURE;
break;
}
// Gate 4: the same survivors, in whatever order. This is the most
// that can be asked of compactAtomic and asking for more would be a
// test that fails on scheduling.
CUDA_CHECK(cudaMemcpy(h_atomic.data(), d_atomic, mWant * sizeof(int),
cudaMemcpyDeviceToHost));
if (checkSameMultiset("compactAtomic", h_atomic.data(), h_want.data(),
mWant) != EXIT_SUCCESS) {
status = EXIT_FAILURE;
break;
}
outOfOrder[which] =
countOutOfOrder(h_atomic.data(), h_want.data(), mWant);
compactMs[which] = timeKernel([&] {
markNonZero<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, kElems);
launchScan(d_rank, d_sums0, d_sums1, d_sums2);
scatterCompact<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, d_out,
kElems);
});
atomicMs[which] = timeKernel([&] {
zeroCount<<<1, 1>>>(d_count);
compactAtomic<<<kBlocks, kThreadsPerBlock>>>(d_in, d_atomic,
d_count, kElems);
});
if (which != kMidInput) {
continue;
}
// Everything below runs once, on the middle input, because these
// rows are about where the time goes rather than about how the
// survival rate moves it.
copyMs = timeKernel([&] {
copyInts<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
});
markMs = timeKernel([&] {
markNonZero<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, kElems);
});
markScanMs = timeKernel([&] {
markNonZero<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, kElems);
launchScan(d_rank, d_sums0, d_sums1, d_sums2);
});
// The row above left d_rank holding the scanned ranks, so this one
// times the scatter on its own rather than the scatter plus the work
// that fed it.
scatterMs = timeKernel([&] {
scatterCompact<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, d_out,
kElems);
});
const size_t evenWant =
partitionCpu(h_in.data(), h_want.data(), kElems);
markEven<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, kElems);
launchScan(d_rank, d_sums0, d_sums1, d_sums2);
scatterPartition<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, d_out,
kElems, d_sums2);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
int scannedEven = 0;
CUDA_CHECK(cudaMemcpy(&scannedEven, d_sums2, sizeof(int),
cudaMemcpyDeviceToHost));
if (static_cast<size_t>(scannedEven) != evenWant) {
std::fprintf(stderr,
"the scan counted %d even elements, the host counted "
"%zu\n",
scannedEven, evenWant);
status = EXIT_FAILURE;
break;
}
// Gate 5: the partition is a permutation of the input with both
// halves in order, checked over all kElems and not only the first
// half, because an odd element landing one slot early is exactly the
// bug the i - rank[i] line can have.
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
if (checkExact("partitionEvenOdd", h_out.data(), h_want.data(),
kElems) != EXIT_SUCCESS) {
status = EXIT_FAILURE;
break;
}
evenCount = evenWant;
partitionMs = timeKernel([&] {
markEven<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, kElems);
launchScan(d_rank, d_sums0, d_sums1, d_sums2);
scatterPartition<<<kBlocks, kThreadsPerBlock>>>(d_in, d_rank, d_out,
kElems, d_sums2);
});
}
// Arithmetic, not a claim about the card: a zero or negative elapsed time
// means an event pair never separated, and every ratio derived from it
// would be infinite or nonsense.
if (status == EXIT_SUCCESS) {
const float once[5] = {copyMs, markMs, markScanMs, scatterMs,
partitionMs};
for (int r = 0; r < 5; ++r) {
if (once[r] <= 0.0f) {
std::fprintf(stderr, "timing broken: row %d measured %.6f ms\n",
r, static_cast<double>(once[r]));
status = EXIT_FAILURE;
}
}
for (int which = 0; which < kNumInputs; ++which) {
if (compactMs[which] <= 0.0f || atomicMs[which] <= 0.0f) {
std::fprintf(stderr,
"timing broken on input %d: %.6f and %.6f ms\n",
which, static_cast<double>(compactMs[which]),
static_cast<double>(atomicMs[which]));
status = EXIT_FAILURE;
}
}
}
if (status == EXIT_SUCCESS) {
const double copyGBs = 2.0 * static_cast<double>(bytes) /
(static_cast<double>(copyMs) * 1.0e-3) / 1.0e9;
std::printf("Where the time goes, on the %.1f percent input\n",
keptPercent[kMidInput]);
std::printf("%-36s %9s %11s %8s\n", "row", "launches", "time (ms)",
"x copy");
std::printf("%-36s %9s %11s %8s\n", "-----------------------------",
"--------", "----------", "-------");
std::printf("%-36s %9d %11.3f %8.2f\n", "copy in to out, no compaction",
1, static_cast<double>(copyMs), 1.0);
std::printf("%-36s %9d %11.3f %8.2f\n", "mark the predicate", 1,
static_cast<double>(markMs),
static_cast<double>(markMs / copyMs));
std::printf("%-36s %9d %11.3f %8.2f\n", "mark, then scan", 6,
static_cast<double>(markScanMs),
static_cast<double>(markScanMs / copyMs));
std::printf("%-36s %9d %11.3f %8.2f\n", "scatter alone, ranks ready", 1,
static_cast<double>(scatterMs),
static_cast<double>(scatterMs / copyMs));
std::printf("%-36s %9d %11.3f %8.2f\n", "compaction, all three phases",
7, static_cast<double>(compactMs[kMidInput]),
static_cast<double>(compactMs[kMidInput] / copyMs));
std::printf("%-36s %9d %11.3f %8.2f\n", "partition even then odd", 7,
static_cast<double>(partitionMs),
static_cast<double>(partitionMs / copyMs));
std::printf("%-36s %9d %11.3f %8.2f\n",
"compaction, one atomic counter", 2,
static_cast<double>(atomicMs[kMidInput]),
static_cast<double>(atomicMs[kMidInput] / copyMs));
std::printf("\nThe copy row moved %zu bytes each way at %.1f GB/s.\n",
bytes, copyGBs);
std::printf("The partition put %zu even elements first.\n\n",
evenCount);
std::printf("Survival rate sweep\n");
std::printf("%6s %12s %11s %13s %9s %14s\n", "kept", "survivors",
"scan (ms)", "atomic (ms)", "scan/atom", "out of order");
std::printf("%6s %12s %11s %13s %9s %14s\n", "-----", "-----------",
"----------", "-----------", "--------", "-------------");
for (int which = 0; which < kNumInputs; ++which) {
std::printf("%5.1f%% %12zu %11.3f %13.3f %9.2f %14zu\n",
keptPercent[which], survivors[which],
static_cast<double>(compactMs[which]),
static_cast<double>(atomicMs[which]),
static_cast<double>(compactMs[which] / atomicMs[which]),
outOfOrder[which]);
}
std::printf(
"\nEvery row above was checked before it was timed. The scanned "
"output\nhad to match the host element for element; the atomic "
"output only had\nto hold the same values, and the last column "
"counts how many of its\nslots hold something else.\n");
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_rank));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_atomic));
CUDA_CHECK(cudaFree(d_sums0));
CUDA_CHECK(cudaFree(d_sums1));
CUDA_CHECK(cudaFree(d_sums2));
CUDA_CHECK(cudaFree(d_count));
return status;
}