code/day25-reduction-2/reduction_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 25: parallel reduction, versions 5 to 8.
//
// Day 24 ended at v4, a shared-memory tree with sequential addressing and two
// elements per thread. That kernel is rebuilt here as row two so this page's
// rows are comparable to each other, then given a grid-stride loop, then
// carried up four more rungs and put beside the library call.
//
// v4 reduceSequentialTwoLoads day 24's endpoint, two elements a thread
// v4c reduceCascaded the same tree over a grid-stride loop
// v5 reduceShuffleTail the last six levels in registers
// v6 reduceUnrolled v5 with the block size as a template arg
// v7 reduceVectorized v6 with float4 loads
// v8 two launches v7, then one block over the partials
//
// A copy kernel runs first as the ceiling, measured on the same card in the
// same process. A reduction is bandwidth bound, so the only useful question
// about one is how close to the card's streaming rate it gets, and a ceiling
// quoted from another page or another card does not answer it.
//
// Every reduction row reads the same kElems floats and every GB/s figure
// comes from one constant, so no row can look fast by reading less. That
// discipline is day 11's, and so is the bug it exists to catch.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o reduction_2 reduction_2.cu
// Run: ./reduction_2
//
// Verified: run on a Tesla T4 (compute capability 7.5), driver 595.84,
// CUDA 12.6 (V12.6.85), on 2026-08-30. The only output that may be
// published as this program's output is the transcript in
// evidence/run-2026-08-30.txt. See research/REVIEW-PROCESS.md.
#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. `warpSize` is a built-in run-time
// value, so it cannot size an array or fill a template argument. This can.
constexpr int kWarpSize = 32;
constexpr unsigned int kFullMask = 0xffffffffu; // all 32 lanes, day 23's name
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// A fixed grid, not one derived from kElems. Every thread then walks the same
// number of elements, which is what day 24 called cascading and what makes
// the exactness argument below hold.
constexpr int kBlocks = 1024;
// 2^26 floats, 256 MiB. Large enough that the load loop, not the tree, is
// what the clock sees: each thread reads 256 elements and then folds eight
// levels once.
constexpr size_t kElems = 64ull * 1024ull * 1024ull;
constexpr size_t kTotalThreads =
static_cast<size_t>(kBlocks) * kThreadsPerBlock;
constexpr size_t kElemsPerThread = kElems / kTotalThreads;
// Day 24's v4 gives every thread exactly two elements, so its grid is fixed
// by the problem size rather than chosen. It writes far more partials than
// the cascading rows do, and one buffer has to hold the larger of the two.
constexpr int kTwoLoadBlocks =
static_cast<int>(kElems / (2 * static_cast<size_t>(kThreadsPerBlock)));
constexpr size_t kPartialElems = static_cast<size_t>(kTwoLoadBlocks);
// Counts the halving levels a block of `blockSize` threads folds through.
// constexpr so the claims below are checked by the compiler instead of
// asserted in a comment.
constexpr int treeLevels(int blockSize) {
int levels = 0;
for (int half = blockSize / 2; half > 0; half >>= 1) {
++levels;
}
return levels;
}
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"the tree halves the stride, so the block size must be a power "
"of two");
static_assert(kThreadsPerBlock >= 2 * kWarpSize,
"the shuffle tail folds 64 values into 32 before its first "
"shuffle, so a block of one warp would read past the live part "
"of the tile");
static_assert(treeLevels(kThreadsPerBlock) == 8,
"256 threads fold in eight levels");
static_assert(treeLevels(kThreadsPerBlock) - treeLevels(kWarpSize) - 1 == 2,
"only the first two of those levels need a block-wide barrier; "
"the other six are one warp talking to itself, which is what v5 "
"moves into registers");
static_assert(kElems % kTotalThreads == 0,
"every thread must take the same number of elements, or the "
"exactness argument in main() does not hold");
static_assert(kElemsPerThread <= (1ull << 24),
"a thread's running sum of ones must stay under 2^24, which is "
"where float stops counting integers exactly");
static_assert(kElems % (2 * static_cast<size_t>(kThreadsPerBlock)) == 0,
"day 24's v4 covers exactly two elements per thread, so its "
"grid has to divide the problem exactly");
static_assert(kPartialElems >= static_cast<size_t>(kBlocks),
"one partials buffer serves every row, so it is sized for the "
"row that writes the most");
static_assert(kElems % 4 == 0,
"kElems is a multiple of four, so reduceVectorized's scalar "
"tail is empty on this input; the tail is written anyway "
"because the exercise changes kElems");
// Sums the value held by each of one warp's 32 lanes and returns the total in
// lane 0.
//
// One warp: lane L adds lane L + offset's register directly, five times. No
// shared memory, no barrier, no volatile. The mask names every lane because
// every lane of this warp is active at the call site; naming a lane that is
// not active is undefined, and that is the one way this primitive bites.
//
// Why not the volatile shared-memory tail the classic slides use: since
// Volta, threads of a warp can sit at different instructions, so a tail that
// relies on them moving in lock step is not merely fragile, it is undefined.
// snippet: warp-reduce
__device__ float warpReduceSum(float val) {
for (unsigned int offset = kWarpSize / 2; offset > 0; offset >>= 1) {
val += __shfl_down_sync(kFullMask, val, offset);
}
return val;
}
// end snippet
// v4, day 24's endpoint, rebuilt here so this page's rows are comparable to
// each other rather than to another page's problem size. One thread reads
// exactly two elements, so the grid is decided by the problem.
//
// One warp: the two loads are 32 consecutive floats each, 128 contiguous
// bytes apart by blockDim.x, so both are fully coalesced and the second is
// independent of the first, which is the latency hiding day 24 measured.
//
// Launch assumption: exactly kThreadsPerBlock threads per block and a grid of
// kTwoLoadBlocks. The guards are still here because the exercise changes
// kElems and the grid may then not divide it.
__global__ void reduceSequentialTwoLoads(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) * 2 + tid;
float acc = (i < n) ? in[i] : 0.0f;
if (i + blockDim.x < n) {
acc += in[i + blockDim.x];
}
tile[tid] = acc;
__syncthreads();
for (unsigned int half = blockDim.x / 2; half > 0; half >>= 1) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// v4c. The same tree, with the two loads replaced by a grid-stride loop over
// a fixed grid, so one thread folds 256 elements into a register before the
// block folds anything. This is the deck's own last rung and it belongs at
// the bottom of this page rather than in the middle, because every rung above
// is measured against a kernel whose load loop is long enough to matter.
//
// One thread: walks the input with a grid stride, accumulating in a register,
// then joins the block's tree.
//
// One warp: on every step of the load loop the 32 lanes read 32 consecutive
// floats, which is 128 contiguous bytes and four 32-byte sectors. The tree
// walks shared memory with a stride of one word, so no two lanes share a
// bank.
//
// Launch assumption: blockDim.x equals kThreadsPerBlock, which is what sizes
// the tile, and is a power of two, which is what the tree needs. The loop
// bound is read from blockDim.x at run time even though the tile is sized
// from a constant, because C++ needs the array bound at compile time and the
// loop does not. Closing that gap is the whole of v6.
__global__ void reduceCascaded(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid; i < n;
i += step) {
acc += in[i];
}
tile[tid] = acc;
__syncthreads();
for (unsigned int half = blockDim.x / 2; half > 0; half >>= 1) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid == 0) {
out[blockIdx.x] = tile[0];
}
}
// v5. The same kernel with the tree stopped early and the last six levels
// handed to one warp.
//
// One warp: the load loop is unchanged. The tail is warp 0 alone, reading two
// shared words each and then five registers each, while the other seven warps
// have nothing left to do.
//
// Launch assumption: as v4, plus blockDim.x >= 64. The loop leaves 64 live
// values and the tail folds them in pairs before its first shuffle, so a
// one-warp block would read tile entries nothing wrote.
__global__ void reduceShuffleTail(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid; i < n;
i += step) {
acc += in[i];
}
tile[tid] = acc;
__syncthreads();
// snippet: v5-tail
// Stop at 64 live values. Below that the survivors are one warp, and a
// warp does not need a block-wide barrier to talk to itself.
for (unsigned int half = blockDim.x / 2; half > kWarpSize; half >>= 1) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid < kWarpSize) {
const float val = warpReduceSum(tile[tid] + tile[tid + kWarpSize]);
if (tid == 0) {
out[blockIdx.x] = val;
}
}
// end snippet
}
// v6. Byte for byte v5, with `blockDim.x` replaced by a template argument in
// its three uses and the tile sized from the same argument instead of from
// kThreadsPerBlock. Nothing else changes, which is the point: the diff is
// the experiment, and it is four lines.
//
// The block size is now known while the kernel is compiled, so the trip count
// of the tree loop is known too and ptxas can emit the levels as straight
// line code with no loop counter and no comparison on `half`. The classic
// slides write that out by hand as a cascade of `if (kBlockSize >= 512)`
// blocks; a loop over a compile-time bound says the same thing in five lines.
// snippet: v6-template
// The second exception to the no-templates-before-module-8 rule in
// CUDA-CODE-STYLE.md, after timeKernel. It is not a style choice: v6 of the
// canonical ladder IS template unrolling, because the block size has to be a
// compile-time constant for the compiler to erase the loop and the bounds
// tests. A runtime blockDim.x cannot do it. The comparison against a library
// reduction belongs to day 39, where CCCL is introduced.
template <unsigned int kBlockSize>
__global__ void reduceUnrolled(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kBlockSize];
const unsigned int tid = threadIdx.x;
const size_t step = gridDim.x * static_cast<size_t>(kBlockSize);
// end snippet
float acc = 0.0f;
for (size_t i = blockIdx.x * static_cast<size_t>(kBlockSize) + tid; i < n;
i += step) {
acc += in[i];
}
tile[tid] = acc;
__syncthreads();
for (unsigned int half = kBlockSize / 2; half > kWarpSize; half >>= 1) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid < kWarpSize) {
const float val = warpReduceSum(tile[tid] + tile[tid + kWarpSize]);
if (tid == 0) {
out[blockIdx.x] = val;
}
}
}
// v7. v6 with the load loop reading four floats per instruction.
//
// One warp: a float4 is 16 bytes, so the 32 lanes now cover 512 contiguous
// bytes per load instead of 128, and the loop issues a quarter as many load
// instructions for the same bytes.
//
// Launch assumption: as v6. `in` must be 16-byte aligned, which every pointer
// from cudaMalloc is; a pointer into the middle of an array is not, and that
// is the failure this kernel invites.
//
// What it measures and what it does not: the tail below handles an n that is
// not a multiple of four, and on this file's kElems it never runs. It is here
// because the exercise changes kElems, and a vectorised kernel that silently
// drops up to three elements is the easiest wrong answer on this page.
template <unsigned int kBlockSize>
__global__ void reduceVectorized(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
__shared__ float tile[kBlockSize];
const unsigned int tid = threadIdx.x;
const size_t t = blockIdx.x * static_cast<size_t>(kBlockSize) + tid;
const size_t step = gridDim.x * static_cast<size_t>(kBlockSize);
// snippet: v7-load
const size_t n4 = n / 4;
const float4* in4 = reinterpret_cast<const float4*>(in);
float acc = 0.0f;
for (size_t i = t; i < n4; i += step) {
const float4 v = in4[i];
acc += v.x + v.y + v.z + v.w;
}
// The 0 to 3 elements that do not fill a float4, one each to the first
// few threads of the grid.
const size_t tailStart = n4 * 4;
if (tailStart + t < n) {
acc += in[tailStart + t];
}
// end snippet
tile[tid] = acc;
__syncthreads();
for (unsigned int half = kBlockSize / 2; half > kWarpSize; half >>= 1) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid < kWarpSize) {
const float val = warpReduceSum(tile[tid] + tile[tid + kWarpSize]);
if (tid == 0) {
out[blockIdx.x] = val;
}
}
}
// The ceiling. The same bytes in that the reductions read, the same bytes out
// again, no tree and no shared memory.
//
// One warp: 32 lanes read 32 consecutive floats and write 32 consecutive
// floats, four sectors each way, which is the best this access pattern can
// do. Day 11 timed the same shape on this card.
__global__ void copyFloats(const float* __restrict__ in,
float* __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];
}
}
// CPU reference. Written for obvious correctness, not speed: one plain loop,
// no blocking and no OpenMP. It accumulates in double even though every
// kernel accumulates in float, because the reference's job is to be right
// rather than to match bit for bit. It never allocates.
static double sumCpu(const float* x, size_t n) {
double total = 0.0;
for (size_t i = 0; i < n; ++i) {
total += static_cast<double>(x[i]);
}
return total;
}
// Reads `count` floats back and sums them on the host in double. Only ever
// called outside a timed region.
static double readBackSum(const float* d_values, size_t count) {
std::vector<float> h_values(count);
CUDA_CHECK(cudaMemcpy(h_values.data(), d_values, count * sizeof(float),
cudaMemcpyDeviceToHost));
double total = 0.0;
for (size_t i = 0; i < count; ++i) {
total += static_cast<double>(h_values[i]);
}
return total;
}
// Returns 1 when the total is wrong, so the caller can add it to a failure
// count and keep going. The comparison is exact and that is deliberate: see
// the input comment in main().
//
// A real branch, not an assert(). CI builds Release, Release defines NDEBUG,
// and NDEBUG deletes assert(), so a gate written that way is missing from
// exactly the build that matters.
static int reportTotal(const char* name, double got, double want) {
if (got != want) {
std::fprintf(stderr, "%s summed %.17g, want %.17g, off by %.17g\n",
name, got, want, got - want);
return 1;
}
return 0;
}
// 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 reduction row reads kElems floats and writes at most kBlocks, so the
// figure comes from one constant rather than from each kernel, where the rows
// could drift apart.
static double reduceGBs(float ms) {
const double bytes = static_cast<double>(kElems) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// The copy row moves those bytes in and the same bytes out, so its figure is
// twice the reduction's for the same time. The two columns are not the same
// quantity and the footer on the table says so.
static double copyGBs(float ms) {
return 2.0 * reduceGBs(ms);
}
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, %zu bytes of shared memory per block\n",
prop.multiProcessorCount, prop.sharedMemPerBlock);
// Every failure below records itself and falls through to the one cleanup
// block at the bottom, so no path returns with device memory allocated.
int failures = 0;
float* d_in = nullptr;
float* d_out = nullptr;
float* d_partials = nullptr;
float* d_sum = nullptr;
const size_t bytes = kElems * sizeof(float);
const size_t partialBytes = kPartialElems * sizeof(float);
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMalloc(&d_partials, partialBytes));
CUDA_CHECK(cudaMalloc(&d_sum, sizeof(float)));
// Every element is 1.0f, so the answer is kElems and the sum counts how
// many accumulations happened. A kernel that drops one element or takes
// one twice is off by exactly one and the gate below sees it, which a
// tolerance over random data could not: one element in 2^26 is a relative
// error of 1.5e-8, far inside any tolerance worth writing.
//
// The sum is exact in float at every level. Each thread adds
// kElemsPerThread ones, which the static_assert above keeps under 2^24.
// Every value entering the tree is then equal to every other, and adding
// two equal floats only moves the exponent, so no level rounds.
//
// What this input cannot catch on its own: a kernel that reads the right
// number of elements from the wrong addresses. The bandwidth column is
// the other half of that gate, because a kernel rereading one cached
// element 2^26 times would post a figure the memory system cannot deliver.
std::vector<float> h_in(kElems, 1.0f);
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
const double want = sumCpu(h_in.data(), kElems);
// Correctness first, with no clock running. Every row gets one untimed
// launch and its total is compared against the CPU reference exactly. A
// kernel that cannot produce the right answer never reaches the table.
reduceSequentialTwoLoads<<<kTwoLoadBlocks, kThreadsPerBlock>>>(
d_in, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures += reportTotal("v4 reduceSequentialTwoLoads",
readBackSum(d_partials, kTwoLoadBlocks), want);
reduceCascaded<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures += reportTotal("v4c reduceCascaded",
readBackSum(d_partials, kBlocks), want);
reduceShuffleTail<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures += reportTotal("v5 reduceShuffleTail",
readBackSum(d_partials, kBlocks), want);
reduceUnrolled<kThreadsPerBlock>
<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures += reportTotal("v6 reduceUnrolled",
readBackSum(d_partials, kBlocks), want);
reduceVectorized<kThreadsPerBlock>
<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures += reportTotal("v7 reduceVectorized",
readBackSum(d_partials, kBlocks), want);
// v8 is the same v7 launch followed by one block over its partials, so
// the answer is a single float that never left the device.
// snippet: v8-two-passes
reduceVectorized<kThreadsPerBlock>
<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
reduceUnrolled<kThreadsPerBlock><<<1, kThreadsPerBlock>>>(
d_partials, d_sum, static_cast<size_t>(kBlocks));
CUDA_CHECK(cudaGetLastError());
// end snippet
CUDA_CHECK(cudaDeviceSynchronize());
failures += reportTotal("v8 two passes", readBackSum(d_sum, 1), want);
// The copy's output is checked by reducing it with v4c, which the block
// above has already compared against the CPU reference. On an input where
// every element is equal, a total of kElems over the copy's output is the
// same statement as "every element arrived".
copyFloats<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
reduceCascaded<<<kBlocks, kThreadsPerBlock>>>(d_out, d_partials, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
failures +=
reportTotal("copyFloats", readBackSum(d_partials, kBlocks), want);
std::printf("\nn = %zu floats, %.0f MiB read per reduction\n", kElems,
static_cast<double>(bytes) / (1024.0 * 1024.0));
std::printf(
"v4 runs %d blocks of %d threads, two elements each. Every row below "
"it\nruns %d blocks of %d threads, %zu elements each.\n",
kTwoLoadBlocks, kThreadsPerBlock, kBlocks, kThreadsPerBlock,
kElemsPerThread);
std::printf("%-24s %11s %10s %11s %s\n", "kernel", "time (ms)", "GB/s",
"% of copy", "produces");
std::printf("%-24s %11s %10s %11s %s\n", "------------------------",
"---------", "--------", "---------", "--------");
int rows = 0;
const float copyMs = timeKernel([&] {
copyFloats<<<kBlocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
});
const double ceiling = copyGBs(copyMs);
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "copy (ceiling)", copyMs,
ceiling, 100.0, "n floats");
const float v4Ms = timeKernel([&] {
reduceSequentialTwoLoads<<<kTwoLoadBlocks, kThreadsPerBlock>>>(
d_in, d_partials, kElems);
});
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "v4 two loads/thread", v4Ms,
reduceGBs(v4Ms), 100.0 * reduceGBs(v4Ms) / ceiling, "partials");
const float v4cMs = timeKernel([&] {
reduceCascaded<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
});
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "v4c cascaded loads", v4cMs,
reduceGBs(v4cMs), 100.0 * reduceGBs(v4cMs) / ceiling,
"partials");
const float v5Ms = timeKernel([&] {
reduceShuffleTail<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials,
kElems);
});
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "v5 shuffle tail", v5Ms,
reduceGBs(v5Ms), 100.0 * reduceGBs(v5Ms) / ceiling, "partials");
const float v6Ms = timeKernel([&] {
reduceUnrolled<kThreadsPerBlock>
<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
});
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "v6 template unrolled",
v6Ms, reduceGBs(v6Ms), 100.0 * reduceGBs(v6Ms) / ceiling,
"partials");
const float v7Ms = timeKernel([&] {
reduceVectorized<kThreadsPerBlock>
<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
});
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "v7 float4 loads", v7Ms,
reduceGBs(v7Ms), 100.0 * reduceGBs(v7Ms) / ceiling, "partials");
const float v8Ms = timeKernel([&] {
reduceVectorized<kThreadsPerBlock>
<<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
reduceUnrolled<kThreadsPerBlock><<<1, kThreadsPerBlock>>>(
d_partials, d_sum, static_cast<size_t>(kBlocks));
});
++rows;
std::printf("%-24s %11.3f %10.1f %11.1f %s\n", "v8 two passes", v8Ms,
reduceGBs(v8Ms), 100.0 * reduceGBs(v8Ms) / ceiling,
"one float");
std::printf(
"\nThe copy row's GB/s counts the bytes it reads and the bytes it\n"
"writes. Every reduction row counts only the %.0f MiB it reads,\n"
"because that is all a reduction moves. A reduction above 100 is not\n"
"a broken measurement, it is a kernel that does not have to write.\n",
static_cast<double>(bytes) / (1024.0 * 1024.0));
std::printf(
"v8 is the only row that leaves a single number on the device. v4\n"
"leaves %d partials and the four rows above v8 leave %d, for\n"
"somebody else to finish. What v8 costs is that somebody.\n",
kTwoLoadBlocks, kBlocks);
// The row count is checked so the lesson's table cannot drift from what
// the program prints. A real branch, not an assert(), for the reason
// reportTotal() gives.
const int kExpectedRows = 7;
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kExpectedRows);
++failures;
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_partials));
CUDA_CHECK(cudaFree(d_sum));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}