code/day32-scan-2/scan_2.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 32: the work-efficient scan, and the three kernels that take it past
// one block.
//
// One copy and three tile scans over the same 10,000,000 floats, timed in
// one process on one card:
//
// copy reads n floats and writes n, the ceiling
// hillis-steele day 31's scan at tile scale, one add per element per
// level, no tree
// blelloch, plain up-sweep and down-sweep, 2 * (kTileElems - 1) adds
// blelloch, padded the same tree with a spare word every 32 and one
// more every 1024
//
// Every scan row runs the same three-kernel pipeline and every row is
// checked against the CPU reference before it is timed, so the table
// compares algorithms rather than amounts of work finished. Day 11 is the
// lesson about the other kind of benchmark.
//
// The tile is 4096 elements because a two-level scan reaches kTileElems
// squared elements and no further. At 4096 that is 16,777,216, which holds
// ten million with room over; at 1024 it is 1,048,576, which does not, and
// main() refuses to run rather than print a wrong answer.
//
// What this program does not check: that any variant is faster than any
// other. That is the claim the run exists to test, so encoding it as a gate
// would make the test unfalsifiable.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o scan_2 scan_2.cu
// Run: ./scan_2
//
// 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)
// 32 on every GPU this course targets, and 32 shared memory banks to match.
// main() reads the device's own warpSize back and fails if it disagrees,
// because every count below assumes this value.
constexpr int kWarpSize = 32;
constexpr int kBanks = 32;
// 1024 threads, not the course default of 256. The tile has to be big
// enough that one launch's block sums fit back inside a single tile, which
// is what keeps this to three kernels instead of a recursion, and 4096
// elements under 256 threads would put 16 of them in every thread.
// CUDA-CODE-STYLE.md allows the change where the tile shape forces it.
constexpr int kThreadsPerBlock = 1024;
constexpr int kTileElems = 4096;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kVariants = 3;
// The copy row keeps the course default. It is the ceiling every other row
// is measured against, so it should not be handicapped by a block size the
// scan chose for its tile.
constexpr int kCopyThreads = 256;
constexpr double kMiB = 1024.0 * 1024.0;
// Ten million, from the curriculum row. Not a multiple of the tile: 2,441
// full tiles leave 1,664 elements for the last block, so the bounds check in
// every kernel runs on every launch instead of never.
constexpr size_t kElems = 10000000ull;
// Every fourth element is 2.0f and the rest are 1.0f, so the total is
// 12,500,000. Two properties, both load bearing. The total is under 2^24, so
// every prefix lands on a float exactly and the check below can compare with
// no tolerance at all. And the pattern has period 4, so a kernel that reads
// the right number of elements from the wrong addresses shifts a prefix and
// is caught, which an all-ones input would not do.
constexpr size_t kTwos = (kElems + 3ull) / 4ull;
constexpr size_t kInputTotal = kElems + kTwos;
// log2(kTileElems), the number of levels in one sweep.
constexpr int treeLevels(int n) {
int levels = 0;
while ((1 << levels) < n) {
++levels;
}
return levels;
}
constexpr int kLevels = treeLevels(kTileElems);
// The three padding schemes the model below compares. kPadOne is how most
// code on the web modernises GPU Gems 3's macro: one spare word every 32.
// kPadTwo is what that macro actually intends, a second spare word every
// 1024, and it is the one the padded kernel uses.
constexpr int kPadNone = 0;
constexpr int kPadOne = 1;
constexpr int kPadTwo = 2;
constexpr int padWord(int i, int mode) {
if (mode == kPadOne) {
return i + i / kBanks;
}
if (mode == kPadTwo) {
return i + i / kBanks + i / (kBanks * kBanks);
}
return i;
}
// How many words the padded tile needs. padWord(kTileElems - 1, kPadTwo) is
// the largest index any thread computes.
constexpr int kTileWords =
kTileElems + kTileElems / kBanks + kTileElems / (kBanks * kBanks);
constexpr size_t kHillisBytes = 2ull * kTileElems * sizeof(float);
constexpr size_t kPaddedBytes = kTileWords * sizeof(float);
constexpr size_t kMaxSharedBytes =
(kHillisBytes > kPaddedBytes) ? kHillisBytes : kPaddedBytes;
// The shared word active lane L touches at tree level `stride`.
//
// Lane L writes word (2L + 2) * stride - 1 and reads word
// (2L + 1) * stride - 1, so consecutive lanes sit 2 * stride words apart.
constexpr int treeWord(int lane, int stride) {
return (2 * lane + 2) * stride - 1;
}
// How many lanes the level has work for.
constexpr int activeLanes(int stride) {
return kTileElems / (2 * stride);
}
// The most lanes any single bank must answer, for one warp of one level.
// Every word here is distinct, so there is no broadcast to discount: this is
// the number of separate requests the access splits into, which is what day
// 15 measured on this card.
//
// Counted per warp rather than once for the level because the i / 32 term is
// not linear, so two warps of the same padded level can split differently.
constexpr int warpDegree(int stride, int warp, int lanes, int mode) {
int perBank[kBanks] = {};
const int first = warp * kWarpSize;
for (int lane = first; lane < first + kWarpSize && lane < lanes; ++lane) {
++perBank[padWord(treeWord(lane, stride), mode) % kBanks];
}
int worst = 0;
for (int bank = 0; bank < kBanks; ++bank) {
if (perBank[bank] > worst) {
worst = perBank[bank];
}
}
return worst;
}
// The worst degree over every warp of a level.
constexpr int worstDegree(int stride, int mode) {
const int lanes = activeLanes(stride);
const int warps = (lanes + kWarpSize - 1) / kWarpSize;
int worst = 0;
for (int warp = 0; warp < warps; ++warp) {
const int degree = warpDegree(stride, warp, lanes, mode);
if (degree > worst) {
worst = degree;
}
}
return worst;
}
// Warp-wide shared requests one block issues walking the tree, after
// conflicts split them, summed over both sweeps. A conflict-free warp costs
// one request; a 32-way conflicted warp costs 32. This is the number to
// compare between the padding schemes, and it is the number that says
// whether a work-efficient tree is cheap in shared memory as well as in
// adds.
constexpr int treeRequests(int mode) {
int total = 0;
for (int stride = 1; stride < kTileElems; stride *= 2) {
const int lanes = activeLanes(stride);
const int warps = (lanes + kWarpSize - 1) / kWarpSize;
for (int warp = 0; warp < warps; ++warp) {
total += warpDegree(stride, warp, lanes, mode);
}
}
return 2 * total;
}
// The same count for the Hillis-Steele scan, which has no tree: every level
// touches every word, lanes stay adjacent, and nothing conflicts.
constexpr int hillisRequests() {
return kLevels * (kTileElems / kWarpSize);
}
// Adds per tile. The Blelloch figure is GPU Gems 3's 2 * (n - 1); the
// Hillis-Steele figure counts only the lanes that add rather than copy.
constexpr int blellochAdds() {
return 2 * (kTileElems - 1);
}
constexpr int hillisAdds() {
int total = 0;
for (int stride = 1; stride < kTileElems; stride *= 2) {
total += kTileElems - stride;
}
return total;
}
static_assert(kThreadsPerBlock % kWarpSize == 0,
"the tile loops step by whole warps");
static_assert(kThreadsPerBlock <= 1024,
"1024 threads is the block ceiling on every compute capability "
"this course targets");
static_assert((kTileElems & (kTileElems - 1)) == 0,
"the tree halves the active count at every level, so the tile "
"must be a power of two");
static_assert(kTileElems >= 2 * kWarpSize,
"the conflict model reads at least one whole warp per level");
static_assert(kInputTotal < (1ull << 24),
"every prefix has to land on a float exactly, so the total must "
"stay under 2^24");
static_assert(kElems <= static_cast<size_t>(kTileElems) * kTileElems,
"two levels of tiles reach kTileElems squared elements and no "
"further");
static_assert(padWord(kTileElems - 1, kPadTwo) < kTileWords,
"the padded tile must hold the largest padded index");
// The four claims this file's tables rest on, checked by the compiler rather
// than asserted in prose. The exact values are guarded on the tile size
// because they describe a 4096-element tree; the two orderings below are not,
// because padding can never cost a tree more requests than no padding.
static_assert(kTileElems != 4096 || worstDegree(16, kPadNone) == 32,
"unpadded, the level whose lanes are 32 words apart puts all 32 "
"of them in one bank");
static_assert(kTileElems != 4096 || worstDegree(16, kPadTwo) == 1,
"padding makes that level conflict free");
static_assert(kTileElems != 4096 || worstDegree(64, kPadOne) == 4,
"one spare word every 32 still leaves a 4-way conflict at the "
"deep levels, which is why the macro has a second term");
static_assert(kTileElems != 4096 || worstDegree(64, kPadTwo) == 1,
"the second term clears it");
static_assert(treeRequests(kPadOne) <= treeRequests(kPadNone),
"padding cannot make the tree cost more requests");
static_assert(treeRequests(kPadTwo) <= treeRequests(kPadOne),
"the two-term offset cannot cost more than the one-term one");
// Where logical word i of the tile sits once the tile is padded.
//
// This repeats padWord(i, kPadTwo) instead of calling it. A constexpr
// function is not callable from device code unless nvcc is passed
// --expt-relaxed-constexpr, and the build line this lesson publishes does
// not pass it. The two are kept in step by the static_assert above, which
// fails the build if the padded tile stops holding the largest index.
__device__ int paddedWord(int i) {
return i + i / kBanks + i / (kBanks * kBanks);
}
// Exclusive scan of one tile of kTileElems floats, Hillis-Steele, with the
// tile's total left in blockSums[blockIdx.x].
//
// One thread: loads kTileElems / blockDim.x elements, then runs kLevels
// rounds over the whole tile, then stores the same elements back.
//
// Memory: consecutive threads take consecutive elements on every global load
// and store, so a warp's 32 addresses cover 128 contiguous bytes. In shared
// memory the lanes stay adjacent at every level too, so no level conflicts.
// It pays for that in adds, one per element per level.
//
// Launch assumption: kTileElems elements per block. Every thread reaches
// every barrier; the guards cover the loads and the stores, never the
// barrier.
// snippet: hillis-steele
__global__ void scanTileHillisSteele(const float* __restrict__ in,
float* __restrict__ out,
float* __restrict__ blockSums, size_t n) {
__shared__ float buf[2][kTileElems];
const unsigned int tid = threadIdx.x;
const int first = static_cast<int>(tid);
const int span = static_cast<int>(blockDim.x);
const size_t base = blockIdx.x * static_cast<size_t>(kTileElems);
for (int k = first; k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
buf[0][k] = (i < n) ? in[i] : 0.0f;
}
__syncthreads();
int cur = 0;
for (int stride = 1; stride < kTileElems; stride <<= 1) {
const int nxt = 1 - cur;
for (int k = first; k < kTileElems; k += span) {
buf[nxt][k] = (k >= stride) ? buf[cur][k] + buf[cur][k - stride]
: buf[cur][k];
}
cur = nxt;
__syncthreads();
}
// buf[cur] holds the inclusive scan, so the exclusive answer is the
// element to the left and the tile total is the last inclusive value.
for (int k = first; k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
if (i < n) {
out[i] = (k == 0) ? 0.0f : buf[cur][k - 1];
}
}
if (tid == 0) {
blockSums[blockIdx.x] = buf[cur][kTileElems - 1];
}
}
// end snippet
// The same output, built with Blelloch's balanced tree instead.
//
// One thread: loads its elements, walks kLevels levels up the tree and
// kLevels back down, then stores. The tree does 2 * (kTileElems - 1) adds
// against Hillis-Steele's kLevels * kTileElems, which is the whole reason
// the algorithm exists.
//
// Memory: the global loads and stores are the same contiguous runs as
// above. The shared accesses are not. At level `stride` consecutive active
// lanes sit 2 * stride words apart, and day 15 measured that a warp whose
// lanes are s words apart splits into gcd(s, 32) requests, so from stride 16
// up every lane in the warp lands in one bank.
//
// Launch assumption: kTileElems elements per block, and no early return
// anywhere above a barrier.
__global__ void scanTileBlelloch(const float* __restrict__ in,
float* __restrict__ out,
float* __restrict__ blockSums, size_t n) {
__shared__ float tile[kTileElems];
const unsigned int tid = threadIdx.x;
const int first = static_cast<int>(tid);
const int span = static_cast<int>(blockDim.x);
const size_t base = blockIdx.x * static_cast<size_t>(kTileElems);
for (int k = first; k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
tile[k] = (i < n) ? in[i] : 0.0f;
}
__syncthreads();
// snippet: blelloch-sweeps
// Up-sweep. Level `stride` has kTileElems / (2 * stride) adds to make,
// and a block of kThreadsPerBlock threads takes them in as many passes
// as that needs.
for (int stride = 1; stride < kTileElems; stride <<= 1) {
const int adds = kTileElems / (2 * stride);
for (int k = first; k < adds; k += span) {
const int a = (2 * k + 1) * stride - 1;
const int b = (2 * k + 2) * stride - 1;
tile[b] += tile[a];
}
__syncthreads();
}
// The tile total leaves before the root is cleared. Clearing it is what
// makes the down-sweep produce an exclusive scan, and the total is what
// the second kernel scans.
if (tid == 0) {
blockSums[blockIdx.x] = tile[kTileElems - 1];
tile[kTileElems - 1] = 0.0f;
}
__syncthreads();
// Down-sweep. The same levels in reverse: every node hands its own value
// to its left child and the sum of both children to its right child.
for (int stride = kTileElems / 2; stride > 0; stride >>= 1) {
const int adds = kTileElems / (2 * stride);
for (int k = first; k < adds; k += span) {
const int a = (2 * k + 1) * stride - 1;
const int b = (2 * k + 2) * stride - 1;
const float left = tile[a];
tile[a] = tile[b];
tile[b] += left;
}
__syncthreads();
}
// end snippet
for (int k = first; k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
if (i < n) {
out[i] = tile[k];
}
}
}
// scanTileBlelloch with every shared index sent through paddedWord and the
// array grown to match. Nothing else changes: the same levels, the same
// adds, the same answer.
//
// The padding is not free. It costs kTileWords - kTileElems words of shared
// memory per block and one shift and one add per shared access, because
// neither 33 nor 1057 is a power of two the compiler can fold into the
// index. If it wins, it won against that.
__global__ void scanTileBlellochPadded(const float* __restrict__ in,
float* __restrict__ out,
float* __restrict__ blockSums,
size_t n) {
// snippet: padded-tile
__shared__ float tile[kTileWords];
// end snippet
const unsigned int tid = threadIdx.x;
const int first = static_cast<int>(tid);
const int span = static_cast<int>(blockDim.x);
const size_t base = blockIdx.x * static_cast<size_t>(kTileElems);
for (int k = first; k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
tile[paddedWord(k)] = (i < n) ? in[i] : 0.0f;
}
__syncthreads();
for (int stride = 1; stride < kTileElems; stride <<= 1) {
const int adds = kTileElems / (2 * stride);
for (int k = first; k < adds; k += span) {
const int a = paddedWord((2 * k + 1) * stride - 1);
const int b = paddedWord((2 * k + 2) * stride - 1);
tile[b] += tile[a];
}
__syncthreads();
}
if (tid == 0) {
blockSums[blockIdx.x] = tile[paddedWord(kTileElems - 1)];
tile[paddedWord(kTileElems - 1)] = 0.0f;
}
__syncthreads();
for (int stride = kTileElems / 2; stride > 0; stride >>= 1) {
const int adds = kTileElems / (2 * stride);
for (int k = first; k < adds; k += span) {
const int a = paddedWord((2 * k + 1) * stride - 1);
const int b = paddedWord((2 * k + 2) * stride - 1);
const float left = tile[a];
tile[a] = tile[b];
tile[b] += left;
}
__syncthreads();
}
for (int k = first; k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
if (i < n) {
out[i] = tile[paddedWord(k)];
}
}
}
// The third kernel: out[i] += offsets[blockIdx.x], where offsets is the
// scanned array of block sums.
//
// One thread: adds one float to kTileElems / blockDim.x elements. The block
// covers exactly the tile the scan kernels covered, so the offset is one
// value read once per block and there is no division by the tile size
// anywhere in the kernel.
//
// Memory: consecutive threads take consecutive elements, so the read and the
// write are both 128 contiguous bytes per warp.
// snippet: add-offsets
__global__ void addBlockOffsets(const float* __restrict__ offsets,
float* __restrict__ out, size_t n) {
const float offset = offsets[blockIdx.x];
const size_t base = blockIdx.x * static_cast<size_t>(kTileElems);
const int span = static_cast<int>(blockDim.x);
for (int k = static_cast<int>(threadIdx.x); k < kTileElems; k += span) {
const size_t i = base + static_cast<size_t>(k);
if (i < n) {
out[i] += offset;
}
}
}
// end snippet
// The ceiling. No tree, no shared memory, one element per thread.
//
// Memory: one warp's 32 addresses cover 128 contiguous bytes on both the
// read and the write, which is the pattern day 11 measured at the top of
// this card's bandwidth.
__global__ void copyFloats(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = in[i];
}
}
// CPU reference. Written for obvious correctness, not speed: one loop, no
// OpenMP, no blocking. It never allocates; the caller owns both buffers.
//
// The accumulator is a double even though the kernels accumulate in float,
// which is the house rule. On this input it changes nothing, because every
// prefix is a whole number under 2^24 and float holds those exactly.
static void scanExclusiveCpu(const float* in, float* out, size_t n) {
double running = 0.0;
for (size_t i = 0; i < n; ++i) {
out[i] = static_cast<float>(running);
running += static_cast<double>(in[i]);
}
}
// Returns the first index where the device and the reference differ, or n if
// they agree everywhere.
//
// The comparison is exact and there is no tolerance, which is deliberate.
// Every prefix is a whole number under 2^24, so no ordering of the adds can
// change a single bit, and a kernel that drops or double-counts one element
// is off by exactly one here instead of hiding inside a relative tolerance.
static size_t firstMismatch(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;
}
static void launchTileScan(int variant, int blocks, const float* d_src,
float* d_dst, float* d_sums, size_t n) {
switch (variant) {
case 0:
scanTileHillisSteele<<<blocks, kThreadsPerBlock>>>(d_src, d_dst,
d_sums, n);
break;
case 1:
scanTileBlelloch<<<blocks, kThreadsPerBlock>>>(d_src, d_dst, d_sums,
n);
break;
default:
scanTileBlellochPadded<<<blocks, kThreadsPerBlock>>>(d_src, d_dst,
d_sums, n);
break;
}
}
// One complete scan of n elements, in three launches: scan every tile and
// leave its total behind, scan the totals, add each tile's offset back.
//
// Launches only. There is no synchronise and no error check inside, because
// this runs inside the timed region and a synchronise there would measure
// something other than the kernels. main() checks the launches on the
// untimed call it makes first.
//
// The second launch is the same kernel over the block sums, with one block,
// which is why the whole thing needs blocks <= kTileElems.
// snippet: pipeline
static void scanAll(int variant, const float* d_in, float* d_out, float* d_sums,
float* d_offsets, float* d_total, size_t n, int blocks) {
launchTileScan(variant, blocks, d_in, d_out, d_sums, n);
launchTileScan(variant, 1, d_sums, d_offsets, d_total,
static_cast<size_t>(blocks));
addBlockOffsets<<<blocks, kThreadsPerBlock>>>(d_offsets, d_out, n);
}
// end snippet
// Bytes the pipeline moves for one scan of kElems elements. Counted rather
// than quoted: the two small launches are a rounding error, but they are not
// zero and the page should not have to promise that.
static double pipelineBytes(int blocks) {
const double n = static_cast<double>(kElems);
const double b = static_cast<double>(blocks);
// tile scan reads n and writes n plus b sums; the block-sum scan reads b
// and writes b plus one total; the offset pass reads n and b and writes
// n.
return (4.0 * n + 4.0 * b + 1.0) * sizeof(float);
}
// 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;
}
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), %d SMs, warp size %d\n",
prop.name, prop.major, prop.minor, prop.multiProcessorCount,
prop.warpSize);
std::printf("Max threads per block: %d\n", prop.maxThreadsPerBlock);
std::printf("Shared memory per block: %zu bytes\n", prop.sharedMemPerBlock);
// The gates that can fail before anything is allocated are checked here,
// so this early return has nothing to free. Every failure after the
// allocations records itself and falls through to the one cleanup block
// at the bottom.
int failures = 0;
if (prop.warpSize != kWarpSize) {
std::fprintf(stderr,
"this card reports warp size %d; every count in this "
"program assumes %d\n",
prop.warpSize, kWarpSize);
++failures;
}
if (prop.maxThreadsPerBlock < kThreadsPerBlock) {
std::fprintf(stderr,
"this card takes %d threads per block; the tile needs "
"%d\n",
prop.maxThreadsPerBlock, kThreadsPerBlock);
++failures;
}
if (prop.sharedMemPerBlock < kMaxSharedBytes) {
std::fprintf(stderr,
"the widest kernel wants %zu bytes of shared memory and "
"this card offers %zu per block\n",
kMaxSharedBytes, prop.sharedMemPerBlock);
++failures;
}
const int blocks = static_cast<int>((kElems + kTileElems - 1) / kTileElems);
if (blocks > kTileElems) {
std::fprintf(stderr,
"%zu elements need %d tiles and one tile scans %d, so "
"the block sums do not fit in a single tile; a two-level "
"scan reaches %zu elements and this needs a third "
"level\n",
kElems, blocks, kTileElems,
static_cast<size_t>(kTileElems) * kTileElems);
++failures;
}
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed before any allocation\n",
failures);
return EXIT_FAILURE;
}
std::printf("\nTile model, computed at compile time\n");
std::printf(" %d elements per tile, %d threads per block, %d levels\n",
kTileElems, kThreadsPerBlock, kLevels);
std::printf(" %-20s %11s %18s\n", "kernel", "adds/tile", "requests/tile");
std::printf(" %-20s %11s %18s\n", "--------------------", "-----------",
"------------------");
std::printf(" %-20s %11d %18d\n", "hillis-steele", hillisAdds(),
hillisRequests());
std::printf(" %-20s %11d %18d\n", "blelloch, plain", blellochAdds(),
treeRequests(kPadNone));
std::printf(" %-20s %11d %18d\n", "blelloch, padded", blellochAdds(),
treeRequests(kPadTwo));
std::printf("\nUp-sweep levels. The down-sweep repeats them in reverse.\n");
std::printf(" %8s %8s %7s %8s %8s %8s\n", "stride", "lanes", "warps",
"plain", "i/32", "two-term");
std::printf(" %8s %8s %7s %8s %8s %8s\n", "--------", "--------",
"-------", "--------", "--------", "--------");
for (int stride = 1; stride < kTileElems; stride <<= 1) {
const int lanes = activeLanes(stride);
const int warps = (lanes + kWarpSize - 1) / kWarpSize;
std::printf(" %8d %8d %7d %8d %8d %8d\n", stride, lanes, warps,
worstDegree(stride, kPadNone), worstDegree(stride, kPadOne),
worstDegree(stride, kPadTwo));
}
const size_t inputBytes = kElems * sizeof(float);
const size_t sumBytes = static_cast<size_t>(blocks) * sizeof(float);
float* d_in = nullptr;
float* d_out = nullptr;
float* d_sums = nullptr;
float* d_offsets = nullptr;
float* d_total = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, inputBytes));
CUDA_CHECK(cudaMalloc(&d_out, inputBytes));
CUDA_CHECK(cudaMalloc(&d_sums, sumBytes));
CUDA_CHECK(cudaMalloc(&d_offsets, sumBytes));
CUDA_CHECK(cudaMalloc(&d_total, sizeof(float)));
std::vector<float> h_in(kElems);
std::vector<float> h_out(kElems);
std::vector<float> h_want(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = (i % 4ull == 0ull) ? 2.0f : 1.0f;
}
scanExclusiveCpu(h_in.data(), h_want.data(), kElems);
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), inputBytes, cudaMemcpyHostToDevice));
const double inputMiB = static_cast<double>(inputBytes) / kMiB;
std::printf("\n%zu elements, %d tiles of %d, %.0f MiB in\n", kElems, blocks,
kTileElems, inputMiB);
std::printf("every fourth element is 2.0f, the rest 1.0f, total %zu\n",
kInputTotal);
const int copyBlocks =
static_cast<int>((kElems + kCopyThreads - 1) / kCopyThreads);
const float copyMs = timeKernel(
[&] { copyFloats<<<copyBlocks, kCopyThreads>>>(d_in, d_out, kElems); });
const double copyBytes = 2.0 * static_cast<double>(inputBytes);
const double copyGbps =
copyBytes / (static_cast<double>(copyMs) * 1.0e-3) / 1.0e9;
std::printf("\ncopy, no scan: %.3f ms, %.1f GB/s over %.0f MiB\n", copyMs,
copyGbps, copyBytes / kMiB);
const char* names[kVariants] = {"hillis-steele", "blelloch, plain",
"blelloch, padded"};
std::printf("\n%-20s %12s %14s %10s %11s\n", "kernel", "tile (ms)",
"pipeline (ms)", "GB/s", "% of copy");
std::printf("%-20s %12s %14s %10s %11s\n", "--------------------",
"------------", "--------------", "----------", "-----------");
float pipelineMs[kVariants] = {0.0f, 0.0f, 0.0f};
int rows = 0;
for (int v = 0; v < kVariants; ++v) {
// One untimed pipeline first, checked against the reference, so a
// wrong answer is reported before a number that came from it reaches
// the table.
scanAll(v, d_in, d_out, d_sums, d_offsets, d_total, kElems, blocks);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out.data(), d_out, inputBytes,
cudaMemcpyDeviceToHost));
const size_t bad = firstMismatch(h_out.data(), h_want.data(), kElems);
if (bad != kElems) {
std::fprintf(stderr,
"%s wrong at %zu: got %.9g, want %.9g (block %zu, "
"element %zu of its tile)\n",
names[v], bad, h_out[bad], h_want[bad],
bad / kTileElems, bad % kTileElems);
++failures;
}
float total = 0.0f;
CUDA_CHECK(
cudaMemcpy(&total, d_total, sizeof(float), cudaMemcpyDeviceToHost));
if (total != static_cast<float>(kInputTotal)) {
std::fprintf(stderr,
"%s totalled %.9g, want %zu; that is a difference of "
"%.9g elements\n",
names[v], total, kInputTotal,
static_cast<double>(kInputTotal) - total);
++failures;
}
const float tileMs = timeKernel(
[&] { launchTileScan(v, blocks, d_in, d_out, d_sums, kElems); });
pipelineMs[v] = timeKernel([&] {
scanAll(v, d_in, d_out, d_sums, d_offsets, d_total, kElems, blocks);
});
++rows;
const double gbps = pipelineBytes(blocks) /
(static_cast<double>(pipelineMs[v]) * 1.0e-3) /
1.0e9;
std::printf("%-20s %12.3f %14.3f %10.1f %11.1f\n", names[v], tileMs,
pipelineMs[v], gbps, 100.0 * gbps / copyGbps);
}
std::printf("\npadded is %.2fx the plain tree's pipeline time\n",
static_cast<double>(pipelineMs[2]) /
static_cast<double>(pipelineMs[1]));
std::printf(
"the copy row moves %.0f MiB and a scan pipeline moves %.0f, "
"so the last column compares rates and not times\n",
copyBytes / kMiB, pipelineBytes(blocks) / kMiB);
// The row count is checked so the lesson's table and this program cannot
// quietly disagree. 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.
if (rows != kVariants) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kVariants);
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_sums));
CUDA_CHECK(cudaFree(d_offsets));
CUDA_CHECK(cudaFree(d_total));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}