code/day82-math-libs/math_libs.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 82: cuRAND, cuSPARSE, cuFFT and cuSOLVER, one program, organised by
// what each library replaces rather than by API surface.
//
// Part 1, cuRAND, replaces a random number generator you wrote yourself.
// The same Monte Carlo estimate of pi is computed three ways: a hand-written
// 64-bit LCG, the cuRAND device API with one generator state per thread, and
// the cuRAND host API filling a buffer that a plain kernel then reads. The
// state initialization is timed on its own, because it is the cost people do
// not budget for.
//
// Part 2, cuSPARSE, replaces the SpMV kernel of day 37. Day 37's two
// matrices and both of its kernels are here unchanged, so cusparseSpMV is
// measured against them on the same data in the same process. The skewed
// matrix is the one where day 37 measured 3.8 percent of the issued lane
// slots doing work.
//
// Part 3, cuFFT, replaces a transform nobody writes by hand. It checks a
// forward transform against an analytic spectrum, checks the forward-inverse
// round trip against the input scaled by the transform size (cuFFT does not
// normalize), and times one batched plan against the same work done one
// transform at a time.
//
// Part 4, cuSOLVER, replaces a blocked pivoted LU factorization. There is no
// hand-written comparison here on purpose: that is the point of the row.
//
// Every library call is status checked. The four libraries return four
// different status enums, so each gets its own check macro in the shape of
// CUDA_CHECK. Only cuSPARSE has a function that turns a status into a
// string; the other three print the numeric code, which is the whole reason
// the macro prints the call text as well.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o math_libs math_libs.cu \
// -lcurand -lcufft -lcusparse -lcusolver
// Run: ./math_libs
//
// Verified on a Tesla T4 with CUDA 12.6 on 2026-09-02. The first attempt
// exposed a DC-bin bug in the analytic oracle; both transcripts are retained.
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cuda_runtime.h>
#include <cufft.h>
#include <curand.h>
#include <curand_kernel.h>
#include <cusolverDn.h>
#include <cusparse.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)
// Four more, one per library, in the same shape and for the same reason.
// CUDA_CHECK cannot wrap any of them: a cusparseStatus_t is not a
// cudaError_t, and a status quietly assigned to a variable nobody reads is
// how day 44 shipped a discarded cuBLAS status inside a timing lambda.
//
// cuSPARSE is the only one of the four with a documented status-to-string
// function (https://docs.nvidia.com/cuda/cusparse/index.html section 4.2.4,
// checked 2026-09-01). cuRAND, cuFFT and cuSOLVER have none, so those three
// print the integer and rely on the stringized call for context.
// snippet: library-check-macros
#define CURAND_CHECK(call) \
do { \
curandStatus_t st_ = (call); \
if (st_ != CURAND_STATUS_SUCCESS) { \
std::fprintf(stderr, "cuRAND error %s:%d: %s: status %d\n", \
__FILE__, __LINE__, #call, static_cast<int>(st_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
#define CUFFT_CHECK(call) \
do { \
cufftResult st_ = (call); \
if (st_ != CUFFT_SUCCESS) { \
std::fprintf(stderr, "cuFFT error %s:%d: %s: status %d\n", \
__FILE__, __LINE__, #call, static_cast<int>(st_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
#define CUSPARSE_CHECK(call) \
do { \
cusparseStatus_t st_ = (call); \
if (st_ != CUSPARSE_STATUS_SUCCESS) { \
std::fprintf(stderr, "cuSPARSE error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cusparseGetErrorString(st_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
#define CUSOLVER_CHECK(call) \
do { \
cusolverStatus_t st_ = (call); \
if (st_ != CUSOLVER_STATUS_SUCCESS) { \
std::fprintf(stderr, "cuSOLVER error %s:%d: %s: status %d\n", \
__FILE__, __LINE__, #call, static_cast<int>(st_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// end snippet
// 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;
// Part 1, Monte Carlo. 2^24 points, two draws each. One draw pair per loop
// iteration and 256 iterations per thread, so the generator state is loaded
// once and amortised over 512 draws.
constexpr size_t kMcSamples = 16777216;
constexpr int kMcBlocks = 256;
constexpr size_t kMcThreads = static_cast<size_t>(kMcBlocks) * kThreadsPerBlock;
constexpr size_t kMcPerThread = kMcSamples / kMcThreads;
constexpr unsigned long long kMcSeed = 20260901ull;
constexpr double kPi = 3.14159265358979323846;
// The estimator counts a Bernoulli variable with p = pi/4, and 4 * phat has
// standard error sqrt(pi * (4 - pi) / N). Five of those is the gate.
constexpr double kSigmaGate = 5.0;
static_assert(kMcSamples % kMcThreads == 0,
"every thread must draw the same number of samples, or the "
"estimate is not an average over equal-sized blocks");
static_assert(kThreadsPerBlock % kWarpSize == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"the block reduction halves the stride, so the block size must "
"be a power of two");
// Part 2, SpMV. These four constants and the two matrices they describe are
// day 37's, unchanged, so the timings on this page can be read against
// day 37's own transcript. kHeavyRowLen is derived rather than typed.
constexpr size_t kRows = 524288;
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);
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(kRows % kWarpSize == 0,
"the skewed matrix is built one aligned group of 32 rows at a "
"time");
// Part 3, FFT. 1024 transforms of 1024 complex points, which is 8 MiB in and
// 8 MiB out, small enough that the per-call cost is visible against the
// arithmetic. kFftBin is the frequency of batch 0; batch b uses bin
// (kFftBin + b) % (kFftSize / 2), so a plan that ignores its batch index
// fails the spectrum check instead of passing it.
constexpr int kFftSize = 1024;
constexpr int kFftBatch = 1024;
constexpr int kFftBin = 37;
constexpr size_t kFftPoints =
static_cast<size_t>(kFftSize) * static_cast<size_t>(kFftBatch);
static_assert(kFftBin < kFftSize / 2,
"the base bin must be in the half-spectrum before batching");
// Part 4, dense solve. 1024 by 1024 in double, which is 8 MiB, and the
// f64 row of EXERCISE-DESIGN.md's tolerance table (rtol 1e-12, atol 1e-14)
// is what grades it.
constexpr int kSolveN = 1024;
constexpr double kSolveRelTol = 1e-12;
constexpr double kSolveAbsTol = 1e-14;
// 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;
}
// ---------------------------------------------------------------------
// Part 1: cuRAND against a generator you wrote yourself.
// ---------------------------------------------------------------------
// One block's hits, summed into a single device counter.
//
// Launch requirement: exactly kThreadsPerBlock threads per block, because the
// shared array is sized from that constant and the halving loop counts down
// from it. Every thread reaches every barrier; the guard in the callers
// covers the loads, never the barrier.
__device__ inline void addBlockHits(unsigned long long hits,
unsigned long long* inside) {
__shared__ unsigned long long tile[kThreadsPerBlock];
const unsigned int tid = threadIdx.x;
tile[tid] = hits;
__syncthreads();
for (unsigned int half = kThreadsPerBlock / 2; half > 0; half /= 2) {
if (tid < half) {
tile[tid] += tile[tid + half];
}
__syncthreads();
}
if (tid == 0) {
atomicAdd(inside, tile[0]);
}
}
// One generator state per thread: the same seed for the whole run, and the
// thread's own index as the sequence number.
//
// The documentation is explicit about which knob is which: "For the highest
// quality parallel pseudorandom number generation, each experiment should be
// assigned a unique seed. Within an experiment, each thread of computation
// should be assigned a unique sequence number"
// (https://docs.nvidia.com/cuda/curand/device-api-overview.html , checked
// 2026-09-01). Giving each thread a different seed instead is the common
// mistake: nothing in the algorithm then keeps two threads' streams apart.
//
// One thread does one curand_init. One warp: 32 lanes write 32 adjacent
// states, and a Philox state is 64 bytes, so the warp stores 2 KiB
// contiguously. That is why this is a separate kernel and separately timed:
// it is a store-bound pass over the whole state array.
//
// Launch assumption: gridDim.x * blockDim.x >= nStates.
// snippet: curand-init
__global__ void initSampleStates(curandStatePhilox4_32_10_t* states,
unsigned long long seed, size_t nStates) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (t < nStates) {
curand_init(seed, t, 0, &states[t]);
}
}
// end snippet
// Counts points inside the quarter circle, drawing from the cuRAND device
// API. One thread owns one state and draws kMcPerThread pairs from it.
//
// The state is copied into a register-resident local, used, and written back,
// so the loop does not touch global memory at all. The documentation's
// advice is the same: "saving and restoring generator state in global memory
// between kernel launches" beats recomputing it.
//
// One warp: no global loads inside the loop, so the only warp-level memory
// behaviour is the state load and store either side of it.
//
// Launch assumption: exactly kMcThreads threads, blockDim.x ==
// kThreadsPerBlock, because addBlockHits needs both.
__global__ void countInsideCurand(curandStatePhilox4_32_10_t* states,
unsigned long long* inside,
size_t perThread) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
curandStatePhilox4_32_10_t state = states[t];
unsigned long long hits = 0;
for (size_t k = 0; k < perThread; ++k) {
const float x = curand_uniform(&state);
const float y = curand_uniform(&state);
if (x * x + y * y <= 1.0f) {
++hits;
}
}
states[t] = state;
addBlockHits(hits, inside);
}
// A 64-bit linear congruential generator, the thing people write when they
// need "just a random number" in a kernel. The multiplier and increment are
// Knuth's MMIX values. Only the top bits of an LCG have the full period, so
// the float comes from bits 63 to 40 and nothing below them.
__device__ inline unsigned long long lcgNext(unsigned long long* state) {
*state = *state * 6364136223846793005ull + 1442695040888963407ull;
return *state;
}
__device__ inline float lcgUniform(unsigned long long* state) {
const unsigned long long bits = lcgNext(state) >> 40;
return static_cast<float>(bits) * (1.0f / 16777216.0f);
}
// The same count with the hand-written generator, seeded by thread index
// through a mixing constant so two threads do not start on the same value.
//
// What this row measures and what it does not: it measures how fast a
// two-instruction generator draws, and it does not measure whether the draws
// are any good. Nothing on this page gates the estimate this kernel
// produces, because a hand-rolled LCG failing a five-sigma test is a finding
// rather than a build failure.
//
// Launch assumption: as countInsideCurand.
__global__ void countInsideLcg(unsigned long long seed,
unsigned long long* inside, size_t perThread) {
const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
unsigned long long state = seed + t * 0x9e3779b97f4a7c15ull;
unsigned long long hits = 0;
for (size_t k = 0; k < perThread; ++k) {
const float x = lcgUniform(&state);
const float y = lcgUniform(&state);
if (x * x + y * y <= 1.0f) {
++hits;
}
}
addBlockHits(hits, inside);
}
// The host-API path: the numbers already exist in two device buffers, so
// this kernel only reads them.
//
// One warp: 32 lanes read 32 consecutive floats from each of xs and ys, four
// 32-byte sectors each, which is day 11's coalesced case. The cost this row
// carries that the device-API row does not is 128 MiB of traffic: the
// generator wrote it and this kernel reads it back.
//
// Launch assumption: blockDim.x == kThreadsPerBlock. The grid-stride loop
// condition is the bounds check.
__global__ void countInsideBuffer(const float* __restrict__ xs,
const float* __restrict__ ys,
unsigned long long* inside, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
unsigned long long hits = 0;
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
const float x = xs[i];
const float y = ys[i];
if (x * x + y * y <= 1.0f) {
hits += 1;
}
}
addBlockHits(hits, inside);
}
// Reads the counter back and turns it into an estimate of pi.
static double piFromHits(unsigned long long hits, size_t n) {
return 4.0 * static_cast<double>(hits) / static_cast<double>(n);
}
static int runMonteCarlo() {
int failures = 0;
const double se =
std::sqrt(kPi * (4.0 - kPi) / static_cast<double>(kMcSamples));
const double gate = kSigmaGate * se;
std::printf("\nPart 1: pi from %zu points, %zu threads, %zu pairs each\n",
kMcSamples, kMcThreads, kMcPerThread);
std::printf("gate: |estimate - pi| <= %.1f * sqrt(pi*(4-pi)/N) = %.3e\n",
kSigmaGate, gate);
unsigned long long* d_inside = nullptr;
curandStatePhilox4_32_10_t* d_states = nullptr;
float* d_xs = nullptr;
float* d_ys = nullptr;
const size_t stateBytes = kMcThreads * sizeof(curandStatePhilox4_32_10_t);
const size_t drawBytes = kMcSamples * sizeof(float);
CUDA_CHECK(cudaMalloc(&d_inside, sizeof(unsigned long long)));
CUDA_CHECK(cudaMalloc(&d_states, stateBytes));
CUDA_CHECK(cudaMalloc(&d_xs, drawBytes));
CUDA_CHECK(cudaMalloc(&d_ys, drawBytes));
std::printf("cuRAND device state array: %zu bytes for %zu threads\n",
stateBytes, kMcThreads);
std::printf("cuRAND host draw buffers: %zu bytes\n", 2 * drawBytes);
const int stateBlocks = static_cast<int>(
(kMcThreads + kThreadsPerBlock - 1) / kThreadsPerBlock);
// Row 1: the hand-written generator.
CUDA_CHECK(cudaMemset(d_inside, 0, sizeof(unsigned long long)));
countInsideLcg<<<kMcBlocks, kThreadsPerBlock>>>(kMcSeed, d_inside,
kMcPerThread);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
unsigned long long h_hits = 0;
CUDA_CHECK(cudaMemcpy(&h_hits, d_inside, sizeof(unsigned long long),
cudaMemcpyDeviceToHost));
const double lcgPi = piFromHits(h_hits, kMcSamples);
const float lcgMs = timeKernel([&] {
countInsideLcg<<<kMcBlocks, kThreadsPerBlock>>>(kMcSeed, d_inside,
kMcPerThread);
});
// Row 2: the cuRAND device API, with the state setup timed on its own.
// curand_init is deterministic in its three arguments, so the last of
// timeKernel's launches leaves exactly the states the sampling kernel
// below wants; there is no separate untimed init call to get out of step
// with this one.
const float initMs = timeKernel([&] {
initSampleStates<<<stateBlocks, kThreadsPerBlock>>>(d_states, kMcSeed,
kMcThreads);
});
CUDA_CHECK(cudaMemset(d_inside, 0, sizeof(unsigned long long)));
countInsideCurand<<<kMcBlocks, kThreadsPerBlock>>>(d_states, d_inside,
kMcPerThread);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(&h_hits, d_inside, sizeof(unsigned long long),
cudaMemcpyDeviceToHost));
const double devicePi = piFromHits(h_hits, kMcSamples);
// The states advance on every timed run, so this row times the same work
// on a different part of the stream, never the same draws twice.
const float deviceMs = timeKernel([&] {
countInsideCurand<<<kMcBlocks, kThreadsPerBlock>>>(d_states, d_inside,
kMcPerThread);
});
// Row 3: the cuRAND host API. The generator is created once and the
// creation is outside every timed region, like CUB's temp storage on
// day 39.
curandGenerator_t gen = nullptr;
CURAND_CHECK(curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_PHILOX4_32_10));
CURAND_CHECK(curandSetPseudoRandomGeneratorSeed(gen, kMcSeed));
CURAND_CHECK(curandGenerateUniform(gen, d_xs, kMcSamples));
CURAND_CHECK(curandGenerateUniform(gen, d_ys, kMcSamples));
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemset(d_inside, 0, sizeof(unsigned long long)));
countInsideBuffer<<<kMcBlocks, kThreadsPerBlock>>>(d_xs, d_ys, d_inside,
kMcSamples);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(&h_hits, d_inside, sizeof(unsigned long long),
cudaMemcpyDeviceToHost));
const double hostPi = piFromHits(h_hits, kMcSamples);
const float genMs = timeKernel([&] {
CURAND_CHECK(curandGenerateUniform(gen, d_xs, kMcSamples));
CURAND_CHECK(curandGenerateUniform(gen, d_ys, kMcSamples));
});
const float countMs = timeKernel([&] {
countInsideBuffer<<<kMcBlocks, kThreadsPerBlock>>>(d_xs, d_ys, d_inside,
kMcSamples);
});
const double draws = 2.0 * static_cast<double>(kMcSamples);
std::printf("%-22s %9s %9s %9s %10s %10s %8s\n", "generator", "setup ms",
"draw ms", "total ms", "Mdraw/s", "pi", "sigma");
std::printf("%-22s %9s %9s %9s %10s %10s %8s\n", "----------------------",
"---------", "---------", "---------", "----------",
"----------", "--------");
std::printf("%-22s %9s %9.3f %9.3f %10.1f %10.6f %8.2f\n", "hand LCG", "-",
lcgMs, lcgMs, draws / (lcgMs * 1.0e-3) / 1.0e6, lcgPi,
(lcgPi - kPi) / se);
std::printf("%-22s %9.3f %9.3f %9.3f %10.1f %10.6f %8.2f\n",
"cuRAND device API", initMs, deviceMs, initMs + deviceMs,
draws / (deviceMs * 1.0e-3) / 1.0e6, devicePi,
(devicePi - kPi) / se);
std::printf("%-22s %9.3f %9.3f %9.3f %10.1f %10.6f %8.2f\n",
"cuRAND host API", genMs, countMs, genMs + countMs,
draws / (genMs * 1.0e-3) / 1.0e6, hostPi, (hostPi - kPi) / se);
std::printf(
"setup is curand_init for the device row and nothing for the LCG;\n"
"the host row's setup column is the generate call itself, and its\n"
"draw column is the kernel that reads the %zu MiB back.\n",
(2 * drawBytes) >> 20);
// Both cuRAND rows are gated. The LCG row is reported and not gated,
// for the reason in countInsideLcg's comment.
if (std::fabs(devicePi - kPi) > gate) {
std::fprintf(stderr,
"cuRAND device API estimate %.9g is %.2f sigma from pi, "
"gate is %.1f\n",
devicePi, (devicePi - kPi) / se, kSigmaGate);
++failures;
}
if (std::fabs(hostPi - kPi) > gate) {
std::fprintf(stderr,
"cuRAND host API estimate %.9g is %.2f sigma from pi, "
"gate is %.1f\n",
hostPi, (hostPi - kPi) / se, kSigmaGate);
++failures;
}
CURAND_CHECK(curandDestroyGenerator(gen));
CUDA_CHECK(cudaFree(d_inside));
CUDA_CHECK(cudaFree(d_states));
CUDA_CHECK(cudaFree(d_xs));
CUDA_CHECK(cudaFree(d_ys));
return failures;
}
// ---------------------------------------------------------------------
// Part 2: cuSPARSE against day 37's two kernels, on day 37's two matrices.
// ---------------------------------------------------------------------
// One matrix in CSR on the host, from day 37. Row r owns colIdx and val over
// [rowPtr[r], rowPtr[r + 1]). The offsets are `int`, which is the 32-bit
// cuSPARSE index path this program then asks for by name.
struct HostCsr {
std::vector<int> rowPtr;
std::vector<int> colIdx;
std::vector<float> val;
};
struct DeviceCsr {
int* d_rowPtr;
int* d_colIdx;
float* d_val;
};
// Row lengths for one of the two matrices, from day 37. 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;
}
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 directly in CSR, from day 37's column rule: row r holds
// `len` consecutive columns starting at column r, shifted left only where the
// row would run past the last column. Sorted inside the row, which is what
// cuSPARSE's CSR SpMV expects and what day 37's kernels assumed anyway.
static void buildCsr(const std::vector<int>& lengths, HostCsr* csr) {
csr->rowPtr.assign(kRows + 1, 0);
for (size_t r = 0; r < kRows; ++r) {
csr->rowPtr[r + 1] = csr->rowPtr[r] + lengths[r];
}
csr->colIdx.assign(kNnz, 0);
csr->val.assign(kNnz, 0.0f);
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);
}
size_t slot = static_cast<size_t>(csr->rowPtr[r]);
for (int j = 0; j < len; ++j) {
csr->colIdx[slot] = static_cast<int>(start) + j;
csr->val[slot] = static_cast<float>(j % 7 + 1);
++slot;
}
}
}
// Lane-slots the thread-per-row kernel issues, from day 37: a warp runs until
// its longest row is finished, so one group of 32 rows costs 32 * max(row
// length). Host arithmetic, printed before any timing, because it is the
// prediction the timings get 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;
}
// y[row] = sum over row's nonzeros of val[j] * x[colIdx[j]], one thread per
// row. Day 37's kernel, unchanged.
//
// One warp: 32 lanes hold 32 consecutive rows, so at step k lane L reads
// val[rowPtr[r0 + L] + k]. On the even matrix those addresses are 16 floats
// apart and the warp touches 32 sectors for one request. On the skewed
// matrix the loop bound is per lane, so the warp keeps issuing until its
// longest row is done.
//
// Launch assumption: gridDim.x * blockDim.x >= nRows, and no barrier, so the
// single guard is the whole bounds check.
__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;
}
}
// The same product with the row given to a warp. Day 37's kernel, unchanged.
//
// One warp: at step k the 32 lanes read val[begin + 32k + lane], 128
// contiguous bytes, four sectors. The imbalance moves from between the lanes
// of one warp to between warps.
//
// Launch assumption: blockDim.x is a whole number of warps, so row = t / 32
// is warp-uniform and every lane named in the shuffle mask reaches the
// shuffle.
__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;
}
}
}
// CPU reference, from day 37. Written for obvious correctness, not speed. 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;
}
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));
}
// The whole cuSPARSE SpMV path. alpha and beta are host scalars because the
// handle is left in its default host pointer mode; passing device pointers
// without switching the mode is the classic first failure with this API.
//
// The descriptors carry what the kernels above carried in their argument
// list and their comments: the index type, the index base, and the value
// type. Nothing about the row-length distribution appears anywhere, which is
// the thing this page is about.
// snippet: cusparse-spmv
static void spmvCusparse(cusparseHandle_t handle, cusparseSpMatDescr_t matA,
cusparseDnVecDescr_t vecX, cusparseDnVecDescr_t vecY,
cusparseSpMVAlg_t alg, void* d_buffer) {
const float alpha = 1.0f;
const float beta = 0.0f;
CUSPARSE_CHECK(cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
&alpha, matA, vecX, &beta, vecY, CUDA_R_32F,
alg, d_buffer));
}
static size_t spmvBufferBytes(cusparseHandle_t handle,
cusparseSpMatDescr_t matA,
cusparseDnVecDescr_t vecX,
cusparseDnVecDescr_t vecY,
cusparseSpMVAlg_t alg) {
const float alpha = 1.0f;
const float beta = 0.0f;
size_t bytes = 0;
CUSPARSE_CHECK(cusparseSpMV_bufferSize(
handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA, vecX, &beta,
vecY, CUDA_R_32F, alg, &bytes));
return bytes;
}
// end snippet
static int runSpmv(cusparseHandle_t handle) {
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 2: y = A x, %zu rows, %zu nonzeros in each matrix\n",
kRows, kNnz);
std::printf("%-8s %8s %8s %14s %8s\n", "matrix", "min row", "max row",
"lane-slots", "useful");
std::printf("%-8s %8s %8s %14s %8s\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;
}
const size_t slots = threadPerRowLaneSlots(lengths[m]);
std::printf(
"%-8s %8d %8d %14zu %7.1f%%\n", kNames[m], shortest, longest, slots,
100.0 * static_cast<double>(kNnz) / static_cast<double>(slots));
}
// Values are 1 to 7 and x is 1 to 13, so every product is at most 91 and
// the longest row sums 419 of them: every partial sum is an integer a
// float holds exactly, whatever order it is summed in. So the comparison
// below is not measuring rounding, it is looking for an indexing bug, and
// the four paths must agree to the last bit. The flip side is that this
// matrix cannot show what CSR_ALG2's determinism buys; real values can.
const size_t rowBytes = kRows * 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_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;
DeviceCsr d_csr[2];
CUDA_CHECK(cudaMalloc(&d_x, rowBytes));
CUDA_CHECK(cudaMalloc(&d_y, rowBytes));
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);
cusparseDnVecDescr_t vecX = nullptr;
cusparseDnVecDescr_t vecY = nullptr;
CUSPARSE_CHECK(cusparseCreateDnVec(&vecX, static_cast<int64_t>(kRows), d_x,
CUDA_R_32F));
CUSPARSE_CHECK(cusparseCreateDnVec(&vecY, static_cast<int64_t>(kRows), d_y,
CUDA_R_32F));
const cusparseSpMVAlg_t kAlgs[2] = {CUSPARSE_SPMV_CSR_ALG1,
CUSPARSE_SPMV_CSR_ALG2};
const char* kAlgNames[2] = {"cusparseSpMV ALG1", "cusparseSpMV ALG2"};
std::printf("\n%-8s %-22s %11s %11s %13s\n", "matrix", "path", "time (ms)",
"GFLOP/s", "buffer bytes");
std::printf("%-8s %-22s %11s %11s %13s\n", "--------",
"----------------------", "-----------", "-----------",
"-------------");
const double flops = 2.0 * static_cast<double>(kNnz);
int rows = 0;
for (int m = 0; m < 2; ++m) {
cusparseSpMatDescr_t matA = nullptr;
CUSPARSE_CHECK(cusparseCreateCsr(
&matA, static_cast<int64_t>(kRows), static_cast<int64_t>(kRows),
static_cast<int64_t>(kNnz), d_csr[m].d_rowPtr, d_csr[m].d_colIdx,
d_csr[m].d_val, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I,
CUSPARSE_INDEX_BASE_ZERO, CUDA_R_32F));
// Day 37's two kernels first, checked then timed.
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 %13s\n", kNames[m],
"thread per row", threadMs,
flops / (threadMs * 1.0e-3) / 1.0e9, "-");
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 %13s\n", kNames[m],
"warp per row", warpMs, flops / (warpMs * 1.0e-3) / 1.0e9,
"-");
for (int a = 0; a < 2; ++a) {
const size_t bufferBytes =
spmvBufferBytes(handle, matA, vecX, vecY, kAlgs[a]);
void* d_buffer = nullptr;
// cudaMalloc(&p, 0) is legal and returns a null pointer, which is
// also what cusparseSpMV wants when it asked for no scratch.
if (bufferBytes > 0) {
CUDA_CHECK(cudaMalloc(&d_buffer, bufferBytes));
}
CUDA_CHECK(cudaMemset(d_y, 0, rowBytes));
spmvCusparse(handle, matA, vecX, vecY, kAlgs[a], d_buffer);
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 %s wrong at row %zu: got %.9g, want %.9g\n",
kNames[m], kAlgNames[a], bad, h_y[bad], h_want[m][bad]);
++failures;
}
const float libMs = timeKernel([&] {
spmvCusparse(handle, matA, vecX, vecY, kAlgs[a], d_buffer);
});
++rows;
std::printf("%-8s %-22s %11.3f %11.1f %13zu\n", kNames[m],
kAlgNames[a], libMs, flops / (libMs * 1.0e-3) / 1.0e9,
bufferBytes);
if (d_buffer != nullptr) {
CUDA_CHECK(cudaFree(d_buffer));
}
}
CUSPARSE_CHECK(cusparseDestroySpMat(matA));
}
std::printf(
"\nAll eight rows above ran the same %.0f flops over the same %zu\n"
"values and %zu column indices. Only the row lengths differ between\n"
"the two blocks, and nothing differs inside a block but the code.\n",
flops, 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 = 8;
if (rows != kExpectedRows) {
std::fprintf(stderr,
"printed %d SpMV rows, expected %d; the lesson's table "
"and this program disagree\n",
rows, kExpectedRows);
++failures;
}
CUSPARSE_CHECK(cusparseDestroyDnVec(vecX));
CUSPARSE_CHECK(cusparseDestroyDnVec(vecY));
for (int m = 0; m < 2; ++m) {
freeCsr(&d_csr[m]);
}
CUDA_CHECK(cudaFree(d_x));
CUDA_CHECK(cudaFree(d_y));
return failures;
}
// ---------------------------------------------------------------------
// Part 3: cuFFT, the plan, the batch and the missing 1/N.
// ---------------------------------------------------------------------
// The frequency batch b carries. The range excludes the Nyquist bin; DC is
// included, where the positive and negative spectral lines coincide.
static int fftBinFor(int batch) {
return (kFftBin + batch) % (kFftSize / 2);
}
static int runFft() {
int failures = 0;
const size_t complexBytes = kFftPoints * sizeof(cufftComplex);
std::vector<cufftComplex> h_signal(kFftPoints);
for (int b = 0; b < kFftBatch; ++b) {
const double k0 = static_cast<double>(fftBinFor(b));
const size_t base = static_cast<size_t>(b) * kFftSize;
for (int t = 0; t < kFftSize; ++t) {
const double phase =
2.0 * kPi * k0 * static_cast<double>(t) / kFftSize;
cufftComplex& z = h_signal[base + static_cast<size_t>(t)];
z.x = static_cast<float>(std::cos(phase));
z.y = 0.0f;
}
}
cufftComplex* d_signal = nullptr;
cufftComplex* d_freq = nullptr;
cufftComplex* d_back = nullptr;
CUDA_CHECK(cudaMalloc(&d_signal, complexBytes));
CUDA_CHECK(cudaMalloc(&d_freq, complexBytes));
CUDA_CHECK(cudaMalloc(&d_back, complexBytes));
CUDA_CHECK(cudaMemcpy(d_signal, h_signal.data(), complexBytes,
cudaMemcpyHostToDevice));
// The fourth argument of cufftPlan1d is the batch count. That is the
// whole difference between the two plans below.
// snippet: fft-plans
cufftHandle planBatched = 0;
cufftHandle planOne = 0;
CUFFT_CHECK(cufftPlan1d(&planBatched, kFftSize, CUFFT_C2C, kFftBatch));
CUFFT_CHECK(cufftPlan1d(&planOne, kFftSize, CUFFT_C2C, 1));
size_t batchedWork = 0;
size_t oneWork = 0;
CUFFT_CHECK(cufftGetSize(planBatched, &batchedWork));
CUFFT_CHECK(cufftGetSize(planOne, &oneWork));
// end snippet
std::printf("\nPart 3: %d transforms of %d complex points\n", kFftBatch,
kFftSize);
std::printf("plan work area: batched %zu bytes, single %zu bytes\n",
batchedWork, oneWork);
CUFFT_CHECK(cufftExecC2C(planBatched, d_signal, d_freq, CUFFT_FORWARD));
CUFFT_CHECK(cufftExecC2C(planBatched, d_freq, d_back, CUFFT_INVERSE));
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<cufftComplex> h_freq(kFftPoints);
std::vector<cufftComplex> h_back(kFftPoints);
CUDA_CHECK(cudaMemcpy(h_freq.data(), d_freq, complexBytes,
cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(h_back.data(), d_back, complexBytes,
cudaMemcpyDeviceToHost));
// The tolerance, with its arithmetic, per EXERCISE-DESIGN.md: an f32
// output that accumulates K terms is graded at
// rtol_K = max(1e-5, 4 * 2^-23 * sqrt(K)) and atol_K = rtol_K * max|ref|.
// Here K is the transform size, so sqrt(K) = 32. A radix-2 FFT's error
// actually grows like log2(K), which is smaller, so the bound has room.
const double rtolK =
std::fmax(1.0e-5, 4.0 * std::pow(2.0, -23.0) * std::sqrt(kFftSize));
const double specAtol = rtolK * (kFftSize / 2.0);
const double roundAtol = rtolK * 1.0;
std::printf("tolerance: rtol_K = max(1e-5, 4*2^-23*sqrt(%d)) = %.3e\n",
kFftSize, rtolK);
std::printf(" spectrum |got-exp| <= %.3e, round trip <= %.3e\n",
specAtol, roundAtol);
// A cosine of amplitude 1 at bin k has spectral lines of N/2 at k and at
// (-k mod N). At DC those two lines coincide and add to N. Nothing here
// compares cuFFT with cuFFT.
double worstSpec = 0.0;
double worstRound = 0.0;
for (int b = 0; b < kFftBatch; ++b) {
const int k0 = fftBinFor(b);
for (int k = 0; k < kFftSize; ++k) {
const size_t i =
static_cast<size_t>(b) * kFftSize + static_cast<size_t>(k);
const int mirror = (kFftSize - k0) % kFftSize;
const double wantRe =
(k == k0 ? kFftSize / 2.0 : 0.0) +
(k == mirror ? kFftSize / 2.0 : 0.0);
const double dRe = std::fabs(h_freq[i].x - wantRe);
const double dIm = std::fabs(h_freq[i].y);
worstSpec = std::fmax(worstSpec, std::fmax(dRe, dIm));
// The round trip comes back scaled by the transform size. Divide
// by N or every amplitude on this page is 1024 times too large.
const double backRe = h_back[i].x / kFftSize;
const double backIm = h_back[i].y / kFftSize;
worstRound = std::fmax(
worstRound, std::fmax(std::fabs(backRe - h_signal[i].x),
std::fabs(backIm - h_signal[i].y)));
}
}
// The scale itself, stated as a number rather than as a belief: the mean
// magnitude of the un-divided round trip against the input.
double sumIn = 0.0;
double sumBack = 0.0;
for (size_t i = 0; i < kFftPoints; ++i) {
sumIn += std::fabs(static_cast<double>(h_signal[i].x));
sumBack += std::fabs(static_cast<double>(h_back[i].x));
}
const double measuredScale = (sumIn > 0.0) ? sumBack / sumIn : 0.0;
std::printf("round trip scale, |ifft(fft(x))| / |x|: %.4f (N is %d)\n",
measuredScale, kFftSize);
if (worstSpec > specAtol) {
std::fprintf(stderr,
"forward transform off the analytic spectrum by %.6g, "
"gate %.6g\n",
worstSpec, specAtol);
++failures;
}
if (worstRound > roundAtol) {
std::fprintf(stderr,
"round trip off the input by %.6g after dividing by %d, "
"gate %.6g\n",
worstRound, kFftSize, roundAtol);
++failures;
}
const float batchedMs = timeKernel([&] {
CUFFT_CHECK(cufftExecC2C(planBatched, d_signal, d_freq, CUFFT_FORWARD));
});
const float loopMs = timeKernel([&] {
for (int b = 0; b < kFftBatch; ++b) {
const size_t off = static_cast<size_t>(b) * kFftSize;
CUFFT_CHECK(cufftExecC2C(planOne, d_signal + off, d_freq + off,
CUFFT_FORWARD));
}
});
std::printf("\n%-26s %11s %13s %9s\n", "forward transform", "time (ms)",
"us per xform", "vs batch");
std::printf("%-26s %11s %13s %9s\n", "--------------------------",
"-----------", "-------------", "---------");
std::printf("%-26s %11.3f %13.3f %9.2f\n", "one batched plan", batchedMs,
batchedMs * 1000.0 / kFftBatch, 1.0);
std::printf("%-26s %11.3f %13.3f %9.2f\n", "batch-1 plan, called 1024x",
loopMs, loopMs * 1000.0 / kFftBatch, loopMs / batchedMs);
std::printf("worst spectrum error %.6g, worst round-trip error %.6g\n",
worstSpec, worstRound);
CUFFT_CHECK(cufftDestroy(planBatched));
CUFFT_CHECK(cufftDestroy(planOne));
CUDA_CHECK(cudaFree(d_signal));
CUDA_CHECK(cudaFree(d_freq));
CUDA_CHECK(cudaFree(d_back));
return failures;
}
// ---------------------------------------------------------------------
// Part 4: cuSOLVER, the factorization nobody should write by hand.
// ---------------------------------------------------------------------
static int runSolve(cusolverDnHandle_t handle) {
int failures = 0;
const size_t n = static_cast<size_t>(kSolveN);
const size_t matrixBytes = n * n * sizeof(double);
const size_t vecBytes = n * sizeof(double);
// Column-major, because that is what every LAPACK-shaped API takes:
// A(i, j) lives at A[i + j * n]. Strictly diagonally dominant, so the
// system is well conditioned and the forward error below is a statement
// about the solver rather than about the matrix.
std::vector<double> h_a(n * n);
std::vector<double> h_xTrue(n);
std::vector<double> h_b(n);
for (size_t j = 0; j < n; ++j) {
for (size_t i = 0; i < n; ++i) {
h_a[i + j * n] = static_cast<double>((i * 7 + j * 13) % 11) / 11.0;
}
}
for (size_t i = 0; i < n; ++i) {
h_a[i + i * n] = static_cast<double>(n);
h_xTrue[i] = 1.0 + 0.25 * static_cast<double>(i % 5);
}
// b = A x, in double on the host, so the exact answer is known and the
// check is not "does the solver agree with itself".
for (size_t i = 0; i < n; ++i) {
double sum = 0.0;
for (size_t j = 0; j < n; ++j) {
sum += h_a[i + j * n] * h_xTrue[j];
}
h_b[i] = sum;
}
double* d_a = nullptr;
double* d_b = nullptr;
double* d_work = nullptr;
int* d_ipiv = nullptr;
int* d_info = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, matrixBytes));
CUDA_CHECK(cudaMalloc(&d_b, vecBytes));
CUDA_CHECK(cudaMalloc(&d_ipiv, n * sizeof(int)));
CUDA_CHECK(cudaMalloc(&d_info, sizeof(int)));
int lwork = 0;
CUSOLVER_CHECK(cusolverDnDgetrf_bufferSize(handle, kSolveN, kSolveN, d_a,
kSolveN, &lwork));
CUDA_CHECK(
cudaMalloc(&d_work, static_cast<size_t>(lwork) * sizeof(double)));
std::printf("\nPart 4: LU factorize and solve a %d by %d double system\n",
kSolveN, kSolveN);
std::printf("cusolverDnDgetrf workspace: %d doubles, %zu bytes\n", lwork,
static_cast<size_t>(lwork) * sizeof(double));
// getrf overwrites A, so it cannot be run in a timing loop without a
// reset, and the reset would sit inside the measurement. One warm call,
// one fresh upload, then one timed call. At this size the factorization
// is milliseconds, far above the half-microsecond event resolution, so a
// single timed call is honest.
CUDA_CHECK(
cudaMemcpy(d_a, h_a.data(), matrixBytes, cudaMemcpyHostToDevice));
CUSOLVER_CHECK(cusolverDnDgetrf(handle, kSolveN, kSolveN, d_a, kSolveN,
d_work, d_ipiv, d_info));
CUDA_CHECK(cudaDeviceSynchronize());
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(
cudaMemcpy(d_a, h_a.data(), matrixBytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaEventRecord(start));
CUSOLVER_CHECK(cusolverDnDgetrf(handle, kSolveN, kSolveN, d_a, kSolveN,
d_work, d_ipiv, d_info));
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float factorMs = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&factorMs, start, stop));
int h_info = 0;
CUDA_CHECK(
cudaMemcpy(&h_info, d_info, sizeof(int), cudaMemcpyDeviceToHost));
if (h_info < 0) {
std::fprintf(stderr,
"cusolverDnDgetrf devInfo = %d: argument %d is wrong\n",
h_info, -h_info);
++failures;
} else if (h_info > 0) {
std::fprintf(stderr,
"cusolverDnDgetrf devInfo = %d: U(%d,%d) is exactly "
"zero, so the matrix is singular to working precision\n",
h_info, h_info, h_info);
++failures;
}
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), vecBytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaEventRecord(start));
CUSOLVER_CHECK(cusolverDnDgetrs(handle, CUBLAS_OP_N, kSolveN, 1, d_a,
kSolveN, d_ipiv, d_b, kSolveN, d_info));
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float solveMs = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&solveMs, start, stop));
CUDA_CHECK(
cudaMemcpy(&h_info, d_info, sizeof(int), cudaMemcpyDeviceToHost));
if (h_info != 0) {
std::fprintf(stderr, "cusolverDnDgetrs devInfo = %d\n", h_info);
++failures;
}
std::vector<double> h_x(n);
CUDA_CHECK(cudaMemcpy(h_x.data(), d_b, vecBytes, cudaMemcpyDeviceToHost));
// Two different questions. The forward error asks whether the answer is
// the one we planted; the residual asks whether it satisfies the system
// we handed over. A solver can look good on one and bad on the other.
double worstFwd = 0.0;
double maxTrue = 0.0;
for (size_t i = 0; i < n; ++i) {
worstFwd = std::fmax(worstFwd, std::fabs(h_x[i] - h_xTrue[i]));
maxTrue = std::fmax(maxTrue, std::fabs(h_xTrue[i]));
}
double worstRes = 0.0;
double maxB = 0.0;
for (size_t i = 0; i < n; ++i) {
double sum = 0.0;
for (size_t j = 0; j < n; ++j) {
sum += h_a[i + j * n] * h_x[j];
}
worstRes = std::fmax(worstRes, std::fabs(sum - h_b[i]));
maxB = std::fmax(maxB, std::fabs(h_b[i]));
}
// The f64 row of EXERCISE-DESIGN.md's table, scaled the way that file
// scales a K-term accumulation: rtol_K = max(1e-12, 4 * 2^-52 * sqrt(n)),
// whose second term is 2.8e-14 at n = 1024, so the table's floor binds.
const double rtolK =
std::fmax(kSolveRelTol, 4.0 * std::pow(2.0, -52.0) * std::sqrt(n));
const double fwdGate = kSolveAbsTol + rtolK * maxTrue;
const double resGate = kSolveAbsTol + rtolK * maxB;
std::printf("tolerance: rtol_K = max(1e-12, 4*2^-52*sqrt(%zu)) = %.3e\n", n,
rtolK);
std::printf("%-26s %11s %14s %14s\n", "stage", "time (ms)", "worst err",
"gate");
std::printf("%-26s %11s %14s %14s\n", "--------------------------",
"-----------", "--------------", "--------------");
std::printf("%-26s %11.3f %14.3e %14.3e\n", "cusolverDnDgetrf", factorMs,
worstFwd, fwdGate);
std::printf("%-26s %11.3f %14.3e %14.3e\n", "cusolverDnDgetrs", solveMs,
worstRes, resGate);
std::printf(
"the getrf row's error column is |x - x_true|, the getrs row's is\n"
"|A x - b|; both come from the one solve those two calls produce.\n");
if (worstFwd > fwdGate) {
std::fprintf(stderr,
"solve is %.6g from the planted answer, gate %.6g\n",
worstFwd, fwdGate);
++failures;
}
if (worstRes > resGate) {
std::fprintf(stderr, "residual |A x - b| is %.6g, gate %.6g\n",
worstRes, resGate);
++failures;
}
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_work));
CUDA_CHECK(cudaFree(d_ipiv));
CUDA_CHECK(cudaFree(d_info));
return failures;
}
// Prints one library's version the way every one of the four reports it, as
// three separate property queries. The transcript then records which shared
// object was actually loaded, which is the check day 44 wishes it had had
// when it linked a binary against the wrong cuBLAS major by accident.
static void printLibVersion(const char* name, int major, int minor, int patch) {
std::printf(" %-10s %d.%d.%d\n", name, major, minor, patch);
}
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);
int runtimeVersion = 0;
int driverVersion = 0;
CUDA_CHECK(cudaRuntimeGetVersion(&runtimeVersion));
CUDA_CHECK(cudaDriverGetVersion(&driverVersion));
std::printf("CUDA runtime %d, driver reports %d\n", runtimeVersion,
driverVersion);
std::printf("libraries linked into this binary:\n");
int major = 0;
int minor = 0;
int patch = 0;
CURAND_CHECK(curandGetProperty(MAJOR_VERSION, &major));
CURAND_CHECK(curandGetProperty(MINOR_VERSION, &minor));
CURAND_CHECK(curandGetProperty(PATCH_LEVEL, &patch));
printLibVersion("cuRAND", major, minor, patch);
CUFFT_CHECK(cufftGetProperty(MAJOR_VERSION, &major));
CUFFT_CHECK(cufftGetProperty(MINOR_VERSION, &minor));
CUFFT_CHECK(cufftGetProperty(PATCH_LEVEL, &patch));
printLibVersion("cuFFT", major, minor, patch);
CUSPARSE_CHECK(cusparseGetProperty(MAJOR_VERSION, &major));
CUSPARSE_CHECK(cusparseGetProperty(MINOR_VERSION, &minor));
CUSPARSE_CHECK(cusparseGetProperty(PATCH_LEVEL, &patch));
printLibVersion("cuSPARSE", major, minor, patch);
CUSOLVER_CHECK(cusolverGetProperty(MAJOR_VERSION, &major));
CUSOLVER_CHECK(cusolverGetProperty(MINOR_VERSION, &minor));
CUSOLVER_CHECK(cusolverGetProperty(PATCH_LEVEL, &patch));
printLibVersion("cuSOLVER", major, minor, patch);
cusparseHandle_t sparseHandle = nullptr;
cusolverDnHandle_t solverHandle = nullptr;
CUSPARSE_CHECK(cusparseCreate(&sparseHandle));
CUSOLVER_CHECK(cusolverDnCreate(&solverHandle));
// Every part records its own failures and falls through to the cleanup
// below, so no path returns with a handle or a device buffer still open.
int failures = 0;
failures += runMonteCarlo();
failures += runSpmv(sparseHandle);
failures += runFft();
failures += runSolve(solverHandle);
CUSPARSE_CHECK(cusparseDestroy(sparseHandle));
CUSOLVER_CHECK(cusolverDnDestroy(solverHandle));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
std::printf("\nall checks passed\n");
return EXIT_SUCCESS;
}