code/day37-spmv/spmv.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 37: sparse matrices, COO to CSR, and SpMV two ways.
//
// The variable this program isolates is the shape of the data, not the code.
// Two matrices are built with the same number of rows, the same number of
// columns and exactly the same number of nonzeros. They differ in one thing,
// how those nonzeros are spread over the rows:
//
// even every row holds kEvenRowLen nonzeros.
// skewed each aligned group of 32 rows holds 31 rows of kLightRowLen and
// one row of kHeavyRowLen. A group of 32 rows is exactly one warp of
// the thread-per-row kernel, so the imbalance lands inside a warp.
//
// Both matrices are multiplied by the same vector with the same two kernels,
// so every row of the timing table performs 2 * kNnz flops and streams the
// same kNnz values and kNnz column indices. Anything that moves between the
// rows moved because of the shape.
//
// A third kernel, probeNonzeroStream, reads the same two arrays with no row
// structure at all, and is the floor the two SpMV kernels are reported
// against. It is called a probe because its output is not a matrix-vector
// product: it is one partial sum per thread over an arbitrary slice of the
// nonzeros. What it measures honestly is the time to stream the matrix once
// with a perfect access pattern.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o spmv spmv.cu
// Run: ./spmv
//
// 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 <cmath>
#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. The built-in `warpSize` is a runtime
// value, so it cannot size an array or appear in a static_assert.
constexpr int kWarpSize = 32;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
// A square matrix, 2^19 rows. x and y are 2 MiB each, which is under the
// T4's 4 MiB of L2, so the gather into x is cheap here. A matrix whose x
// does not fit adds a cost this program does not measure; the README says so.
constexpr size_t kRows = 524288;
// The three row lengths. kHeavyRowLen is derived rather than typed, so the
// two matrices cannot drift into holding different numbers of nonzeros.
constexpr int kEvenRowLen = 16;
constexpr int kLightRowLen = 3;
constexpr int kHeavyRowLen =
kWarpSize * kEvenRowLen - (kWarpSize - 1) * kLightRowLen;
constexpr size_t kNnz = kRows * static_cast<size_t>(kEvenRowLen);
// The probe walks the nonzeros with a fixed grid, 8 blocks per SM on the
// T4's 40 SMs, and one partial sum per thread.
constexpr int kProbeBlocks = 320;
constexpr size_t kProbePartials =
static_cast<size_t>(kProbeBlocks) * kThreadsPerBlock;
static_assert(kThreadsPerBlock % kWarpSize == 0,
"the warp-per-row kernel derives its row from t / 32, which is "
"only warp-uniform when the block is a whole number of warps");
static_assert(kRows % kWarpSize == 0,
"the skewed matrix is built one aligned group of 32 rows at a "
"time");
static_assert((kWarpSize - 1) * kLightRowLen + kHeavyRowLen ==
kWarpSize * kEvenRowLen,
"the two matrices must hold the same number of nonzeros, or "
"nothing in the timing table is comparable");
static_assert(kHeavyRowLen > 0 && static_cast<size_t>(kHeavyRowLen) <= kRows,
"the heaviest row must fit inside the matrix's columns");
// One matrix in CSR on the host. Plain data, no methods.
//
// rowPtr holds kRows + 1 offsets: row r owns colIdx and val over the half
// open range [rowPtr[r], rowPtr[r + 1]). That one array is the whole
// difference from COO, which stores a row index per nonzero instead.
//
// The offsets are `int`, which is what cuSPARSE's 32-bit path uses and what
// almost every published kernel assumes. It caps a matrix at 2^31 nonzeros.
struct HostCsr {
std::vector<int> rowPtr;
std::vector<int> colIdx;
std::vector<float> val;
};
// The same matrix on the device. The fields keep the d_ prefix so a reader
// can still see which side of the bus they name.
struct DeviceCsr {
int* d_rowPtr;
int* d_colIdx;
float* d_val;
};
// Row lengths for one of the two matrices. Deterministic, no RNG: the pair
// exists so that this function is the only difference between them.
static void fillRowLengths(bool skewed, std::vector<int>* lengths) {
for (size_t r = 0; r < kRows; ++r) {
(*lengths)[r] = kEvenRowLen;
}
if (!skewed) {
return;
}
// One heavy row per aligned group of 32, and 31 light rows beside it.
// The group is one warp of the thread-per-row kernel, so 31 lanes stop
// after kLightRowLen steps and wait for the lane that has kHeavyRowLen.
for (size_t g = 0; g < kRows / kWarpSize; ++g) {
const size_t first = g * kWarpSize;
for (int k = 0; k < kWarpSize - 1; ++k) {
(*lengths)[first + static_cast<size_t>(k)] = kLightRowLen;
}
(*lengths)[first + static_cast<size_t>(kWarpSize - 1)] = kHeavyRowLen;
}
}
// Builds the matrix in COO, then converts it to CSR, because that conversion
// is the concrete difference between the two formats.
//
// Columns: row r holds `len` consecutive columns starting at column r, shifted
// left only where the row would otherwise run past the last column. Sorted
// inside the row, and the same rule for both matrices, so the gather into x is
// held fixed while the row lengths change. A real matrix scatters those reads;
// this one does not, on purpose, so the gather cannot explain the difference
// between the two timing blocks.
//
// Values depend on the position inside the row and x depends on the column,
// so a kernel that reads val at the wrong offset and a kernel that reads the
// wrong column produce different wrong answers instead of the same one.
static void buildCsr(const std::vector<int>& lengths, HostCsr* csr) {
std::vector<int> cooRow;
std::vector<int> cooCol;
std::vector<float> cooVal;
cooRow.reserve(kNnz);
cooCol.reserve(kNnz);
cooVal.reserve(kNnz);
for (size_t r = 0; r < kRows; ++r) {
const int len = lengths[r];
size_t start = r;
if (start + static_cast<size_t>(len) > kRows) {
start = kRows - static_cast<size_t>(len);
}
for (int j = 0; j < len; ++j) {
cooRow.push_back(static_cast<int>(r));
cooCol.push_back(static_cast<int>(start) + j);
cooVal.push_back(static_cast<float>(j % 7 + 1));
}
}
// Count each row, prefix sum the counts into offsets, then one pass that
// drops every nonzero into its row's slice. O(nnz) with one cursor array,
// and it works whatever order the COO entries arrive in.
// snippet: coo-to-csr
csr->rowPtr.assign(kRows + 1, 0);
for (size_t e = 0; e < cooRow.size(); ++e) {
++csr->rowPtr[static_cast<size_t>(cooRow[e]) + 1];
}
for (size_t r = 0; r < kRows; ++r) {
csr->rowPtr[r + 1] += csr->rowPtr[r];
}
std::vector<int> cursor(csr->rowPtr.begin(), csr->rowPtr.end() - 1);
csr->colIdx.assign(cooCol.size(), 0);
csr->val.assign(cooVal.size(), 0.0f);
for (size_t e = 0; e < cooRow.size(); ++e) {
const size_t row = static_cast<size_t>(cooRow[e]);
const size_t slot = static_cast<size_t>(cursor[row]++);
csr->colIdx[slot] = cooCol[e];
csr->val[slot] = cooVal[e];
}
// end snippet
}
// Lane-slots the thread-per-row kernel issues: a warp runs until its longest
// row is finished, so one group of 32 rows costs 32 * max(row length).
// Arithmetic on the row lengths, not a measurement, and it is the number the
// timing table is judged against.
static size_t threadPerRowLaneSlots(const std::vector<int>& lengths) {
size_t slots = 0;
for (size_t first = 0; first < kRows; first += kWarpSize) {
int longest = 0;
for (int k = 0; k < kWarpSize; ++k) {
const int len = lengths[first + static_cast<size_t>(k)];
if (len > longest) {
longest = len;
}
}
slots += static_cast<size_t>(longest) * kWarpSize;
}
return slots;
}
// Lane-slots the warp-per-row kernel issues: one row costs 32 lanes times the
// number of 32-wide passes it takes, so a row of 3 costs a whole pass and a
// row of 419 costs fourteen.
static size_t warpPerRowLaneSlots(const std::vector<int>& lengths) {
size_t slots = 0;
for (size_t r = 0; r < kRows; ++r) {
const int passes = (lengths[r] + kWarpSize - 1) / kWarpSize;
slots += static_cast<size_t>(passes) * kWarpSize;
}
return slots;
}
// y[row] = sum over row's nonzeros of val[j] * x[colIdx[j]].
//
// One thread owns one row and walks its slice of colIdx and val.
//
// One warp: 32 lanes hold 32 consecutive rows. At step k, lane L reads
// val[rowPtr[r0 + L] + k]. On the even matrix those addresses are 16 floats
// apart, so the 32 lanes land in 32 different 32-byte sectors and one request
// costs 32 of them, which is day 11's floor. On the skewed matrix the
// addresses are not even a fixed stride apart, and worse, the loop bound is
// per lane, so the warp keeps issuing until its longest row is done.
//
// Launch assumption: gridDim.x * blockDim.x >= nRows. There is no barrier in
// this kernel, so the single guard on the index is the whole bounds check.
// snippet: thread-per-row
__global__ void spmvCsrThreadPerRow(const int* __restrict__ rowPtr,
const int* __restrict__ colIdx,
const float* __restrict__ val,
const float* __restrict__ x,
float* __restrict__ y, size_t nRows) {
const size_t row =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (row < nRows) {
float sum = 0.0f;
for (int j = rowPtr[row]; j < rowPtr[row + 1]; ++j) {
sum += val[j] * x[colIdx[j]];
}
y[row] = sum;
}
}
// end snippet
// The same product, one warp per row.
//
// One thread owns every 32nd nonzero of its warp's row, then the warp adds
// its 32 partial sums with the shuffle reduction from day 23 and lane 0
// writes the answer.
//
// One warp: at step k the 32 lanes read val[begin + 32k + lane], which is 32
// consecutive floats, 128 contiguous bytes, four 32-byte sectors. colIdx is
// read the same way. The imbalance does not disappear, it moves from between
// the lanes of one warp to between warps, where the scheduler can hide it
// behind other resident warps.
//
// Launch assumption: blockDim.x is a whole number of warps, which the
// static_assert on kThreadsPerBlock fixes, and gridDim.x * blockDim.x >=
// nRows * 32. Because a warp's 32 threads hold 32 consecutive values of t,
// row = t / 32 is the same for all of them, so the guard below is
// warp-uniform and every lane named in the shuffle mask reaches the shuffle.
// snippet: warp-per-row
__global__ void spmvCsrWarpPerRow(const int* __restrict__ rowPtr,
const int* __restrict__ colIdx,
const float* __restrict__ val,
const float* __restrict__ x,
float* __restrict__ y, size_t nRows) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row = t / kWarpSize;
const unsigned int lane = threadIdx.x % kWarpSize;
if (row < nRows) {
const int end = rowPtr[row + 1];
float sum = 0.0f;
for (int j = rowPtr[row] + static_cast<int>(lane); j < end;
j += kWarpSize) {
sum += val[j] * x[colIdx[j]];
}
for (int offset = kWarpSize / 2; offset > 0; offset /= 2) {
sum += __shfl_down_sync(0xffffffffu, sum, offset);
}
if (lane == 0) {
y[row] = sum;
}
}
}
// end snippet
// The floor: the same multiply and the same gather over every nonzero, with
// no row structure at all.
//
// One thread walks a grid-stride slice of the nonzeros and keeps a partial
// sum. One warp: 32 lanes read 32 consecutive values and 32 consecutive
// column indices, which is the best access pattern this data admits.
//
// What it measures and what it does not: the time to stream val and colIdx
// once and gather x, which is the ceiling neither SpMV kernel can beat. Its
// output is not a matrix-vector product, so main() checks it against the sum
// of the whole reference vector instead of element by element.
//
// Launch assumption: partials has exactly gridDim.x * blockDim.x entries, so
// the write below needs no guard and the loop condition is the only bounds
// check.
__global__ void probeNonzeroStream(const int* __restrict__ colIdx,
const float* __restrict__ val,
const float* __restrict__ x,
float* __restrict__ partials, size_t nnz) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
float sum = 0.0f;
for (size_t j = t; j < nnz; j += step) {
sum += val[j] * x[colIdx[j]];
}
partials[t] = sum;
}
// CPU reference. Written for obvious correctness, not speed: one row at a
// time, no blocking, no OpenMP. It accumulates in double even though the
// kernels accumulate in float, and it never allocates.
static void spmvCsrCpu(const int* rowPtr, const int* colIdx, const float* val,
const float* x, float* y, size_t nRows) {
for (size_t r = 0; r < nRows; ++r) {
double sum = 0.0;
for (int j = rowPtr[r]; j < rowPtr[r + 1]; ++j) {
sum += static_cast<double>(val[j]) *
static_cast<double>(x[static_cast<size_t>(colIdx[j])]);
}
y[r] = static_cast<float>(sum);
}
}
// Returns the first index where got and want differ by more than the relative
// tolerance, or n if they agree everywhere. Returning the index rather than a
// bool is the point: "wrong at 524256" names the warp, "wrong" does not.
static size_t firstMismatch(const float* got, const float* want, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const float scale = (want[i] == 0.0f) ? 1.0f : std::fabs(want[i]);
if (std::fabs(got[i] - want[i]) > relTolerance * scale) {
return i;
}
}
return n;
}
// Sums in double, because the totals here run past what float counts exactly.
static double totalOf(const float* v, size_t n) {
double total = 0.0;
for (size_t i = 0; i < n; ++i) {
total += static_cast<double>(v[i]);
}
return total;
}
// 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 performs one multiply and one add per nonzero, so the
// flop count comes from one constant rather than from each kernel, where the
// two could drift apart.
static double spmvGFlops(float ms) {
const double flops = 2.0 * static_cast<double>(kNnz);
return flops / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Uploads one CSR matrix. Every cudaMalloc here has its cudaFree in
// freeCsr(), and main() calls that from its single cleanup block.
static void uploadCsr(const HostCsr& h, DeviceCsr* d) {
const size_t ptrBytes = (kRows + 1) * sizeof(int);
const size_t idxBytes = kNnz * sizeof(int);
const size_t valBytes = kNnz * sizeof(float);
CUDA_CHECK(cudaMalloc(&d->d_rowPtr, ptrBytes));
CUDA_CHECK(cudaMalloc(&d->d_colIdx, idxBytes));
CUDA_CHECK(cudaMalloc(&d->d_val, valBytes));
CUDA_CHECK(cudaMemcpy(d->d_rowPtr, h.rowPtr.data(), ptrBytes,
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d->d_colIdx, h.colIdx.data(), idxBytes,
cudaMemcpyHostToDevice));
CUDA_CHECK(
cudaMemcpy(d->d_val, h.val.data(), valBytes, cudaMemcpyHostToDevice));
}
static void freeCsr(DeviceCsr* d) {
CUDA_CHECK(cudaFree(d->d_rowPtr));
CUDA_CHECK(cudaFree(d->d_colIdx));
CUDA_CHECK(cudaFree(d->d_val));
}
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)\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
// 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;
const char* kNames[2] = {"even", "skewed"};
HostCsr h_csr[2];
std::vector<int> lengths[2];
for (int m = 0; m < 2; ++m) {
lengths[m].assign(kRows, 0);
fillRowLengths(m == 1, &lengths[m]);
buildCsr(lengths[m], &h_csr[m]);
}
std::printf("\nPart 1: two matrices, %zu rows, %zu nonzeros each\n", kRows,
kNnz);
std::printf("%-8s %10s %8s %8s %9s\n", "matrix", "nonzeros", "min row",
"max row", "mean row");
std::printf("%-8s %10s %8s %8s %9s\n", "--------", "----------", "--------",
"--------", "---------");
for (int m = 0; m < 2; ++m) {
int shortest = lengths[m][0];
int longest = lengths[m][0];
for (size_t r = 0; r < kRows; ++r) {
if (lengths[m][r] < shortest) {
shortest = lengths[m][r];
}
if (lengths[m][r] > longest) {
longest = lengths[m][r];
}
}
const size_t built = static_cast<size_t>(h_csr[m].rowPtr[kRows]);
if (built != kNnz) {
std::fprintf(stderr,
"%s matrix holds %zu nonzeros, expected %zu; the two "
"matrices are not comparable\n",
kNames[m], built, kNnz);
++failures;
}
std::printf("%-8s %10zu %8d %8d %9.2f\n", kNames[m], built, shortest,
longest,
static_cast<double>(kNnz) / static_cast<double>(kRows));
}
// Lane-slots are host arithmetic over the row lengths, printed before any
// kernel runs, because they are the prediction the timings are judged
// against rather than a result.
std::printf("\nLane-slots issued, and the share with a nonzero to do\n");
std::printf("%-8s %-16s %12s %8s\n", "matrix", "kernel", "lane-slots",
"useful");
std::printf("%-8s %-16s %12s %8s\n", "--------", "----------------",
"------------", "--------");
for (int m = 0; m < 2; ++m) {
const size_t threadSlots = threadPerRowLaneSlots(lengths[m]);
const size_t warpSlots = warpPerRowLaneSlots(lengths[m]);
std::printf("%-8s %-16s %12zu %7.1f%%\n", kNames[m], "thread per row",
threadSlots,
100.0 * static_cast<double>(kNnz) /
static_cast<double>(threadSlots));
std::printf(
"%-8s %-16s %12zu %7.1f%%\n", kNames[m], "warp per row", warpSlots,
100.0 * static_cast<double>(kNnz) / static_cast<double>(warpSlots));
}
const size_t rowBytes = kRows * sizeof(float);
const size_t partialBytes = kProbePartials * sizeof(float);
std::vector<float> h_x(kRows);
for (size_t c = 0; c < kRows; ++c) {
h_x[c] = static_cast<float>(c % 13 + 1);
}
std::vector<float> h_y(kRows);
std::vector<float> h_partials(kProbePartials);
std::vector<float> h_want[2];
for (int m = 0; m < 2; ++m) {
h_want[m].assign(kRows, 0.0f);
spmvCsrCpu(h_csr[m].rowPtr.data(), h_csr[m].colIdx.data(),
h_csr[m].val.data(), h_x.data(), h_want[m].data(), kRows);
}
float* d_x = nullptr;
float* d_y = nullptr;
float* d_partials = nullptr;
DeviceCsr d_csr[2];
CUDA_CHECK(cudaMalloc(&d_x, rowBytes));
CUDA_CHECK(cudaMalloc(&d_y, rowBytes));
CUDA_CHECK(cudaMalloc(&d_partials, partialBytes));
CUDA_CHECK(cudaMemcpy(d_x, h_x.data(), rowBytes, cudaMemcpyHostToDevice));
for (int m = 0; m < 2; ++m) {
uploadCsr(h_csr[m], &d_csr[m]);
}
const int threadBlocks =
static_cast<int>((kRows + kThreadsPerBlock - 1) / kThreadsPerBlock);
const int warpBlocks = static_cast<int>(
(kRows * kWarpSize + kThreadsPerBlock - 1) / kThreadsPerBlock);
std::printf(
"\nPart 2: y = A x, %d timed runs after %d warm-ups, %zu flops per "
"run\n",
kTimedRuns, kWarmupRuns, 2 * kNnz);
std::printf("%-8s %-22s %11s %11s %9s\n", "matrix", "kernel", "time (ms)",
"GFLOP/s", "x floor");
std::printf("%-8s %-22s %11s %11s %9s\n", "--------",
"----------------------", "-----------", "-----------",
"---------");
int rows = 0;
for (int m = 0; m < 2; ++m) {
// The probe first, because it is the floor the other two are read
// against, and it runs on the same card in the same process.
probeNonzeroStream<<<kProbeBlocks, kThreadsPerBlock>>>(
d_csr[m].d_colIdx, d_csr[m].d_val, d_x, d_partials, kNnz);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_partials.data(), d_partials, partialBytes,
cudaMemcpyDeviceToHost));
const double probeTotal = totalOf(h_partials.data(), kProbePartials);
const double wantTotal = totalOf(h_want[m].data(), kRows);
if (std::fabs(probeTotal - wantTotal) > 1.0e-6 * std::fabs(wantTotal)) {
std::fprintf(stderr,
"%s probe total %.17g, reference total %.17g\n",
kNames[m], probeTotal, wantTotal);
++failures;
}
const float probeMs = timeKernel([&] {
probeNonzeroStream<<<kProbeBlocks, kThreadsPerBlock>>>(
d_csr[m].d_colIdx, d_csr[m].d_val, d_x, d_partials, kNnz);
});
++rows;
std::printf("%-8s %-22s %11.3f %11.1f %9.2f\n", kNames[m],
"stream probe (floor)", probeMs, spmvGFlops(probeMs), 1.0);
// One untimed, checked launch per kernel before any number from it
// reaches the table. d_y is cleared first so a row the kernel never
// writes shows up as a mismatch rather than as last kernel's answer.
CUDA_CHECK(cudaMemset(d_y, 0, rowBytes));
spmvCsrThreadPerRow<<<threadBlocks, kThreadsPerBlock>>>(
d_csr[m].d_rowPtr, d_csr[m].d_colIdx, d_csr[m].d_val, d_x, d_y,
kRows);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_y.data(), d_y, rowBytes, cudaMemcpyDeviceToHost));
size_t bad =
firstMismatch(h_y.data(), h_want[m].data(), kRows, kRelTolerance);
if (bad != kRows) {
std::fprintf(stderr,
"%s thread per row wrong at row %zu: got %.9g, want "
"%.9g\n",
kNames[m], bad, h_y[bad], h_want[m][bad]);
++failures;
}
const float threadMs = timeKernel([&] {
spmvCsrThreadPerRow<<<threadBlocks, kThreadsPerBlock>>>(
d_csr[m].d_rowPtr, d_csr[m].d_colIdx, d_csr[m].d_val, d_x, d_y,
kRows);
});
++rows;
std::printf("%-8s %-22s %11.3f %11.1f %9.2f\n", kNames[m],
"thread per row", threadMs, spmvGFlops(threadMs),
threadMs / probeMs);
CUDA_CHECK(cudaMemset(d_y, 0, rowBytes));
spmvCsrWarpPerRow<<<warpBlocks, kThreadsPerBlock>>>(
d_csr[m].d_rowPtr, d_csr[m].d_colIdx, d_csr[m].d_val, d_x, d_y,
kRows);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_y.data(), d_y, rowBytes, cudaMemcpyDeviceToHost));
bad = firstMismatch(h_y.data(), h_want[m].data(), kRows, kRelTolerance);
if (bad != kRows) {
std::fprintf(stderr,
"%s warp per row wrong at row %zu: got %.9g, want "
"%.9g\n",
kNames[m], bad, h_y[bad], h_want[m][bad]);
++failures;
}
const float warpMs = timeKernel([&] {
spmvCsrWarpPerRow<<<warpBlocks, kThreadsPerBlock>>>(
d_csr[m].d_rowPtr, d_csr[m].d_colIdx, d_csr[m].d_val, d_x, d_y,
kRows);
});
++rows;
std::printf("%-8s %-22s %11.3f %11.1f %9.2f\n", kNames[m],
"warp per row", warpMs, spmvGFlops(warpMs),
warpMs / probeMs);
}
std::printf(
"\nAll six rows above ran the same %zu flops over the same %zu\n"
"values and %zu column indices. The two matrices hold the same\n"
"nonzeros in a different set of rows, and nothing else differs.\n",
2 * kNnz, kNnz, kNnz);
// The row count is checked so the lesson's table 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 = 6;
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kExpectedRows);
++failures;
}
for (int m = 0; m < 2; ++m) {
freeCsr(&d_csr[m]);
}
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
CUDA_CHECK(cudaFree(d_partials));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}