code/day43-matmul-1/matmul_steps.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 43: the first three steps of a matmul optimisation ladder, with exactly
// one thing moving per step.
//
// Step 1, matmulStridedRows. One thread per output element, and threadIdx.x
// chooses the row. A warp is then sixteen consecutive rows of C and two
// consecutive columns, so its load of A is sixteen addresses 4N bytes apart.
//
// Step 2, matmulCoalescedCols. The same kernel with threadIdx.x choosing the
// column instead. Two lines move. The launch, the arithmetic, the guard and
// the answer are identical.
//
// Step 3, matmulTiled. Step 2 with both inputs staged through shared memory,
// which cuts global loads per thread from 2N to 2 * ceil(N / 16).
//
// Steps 2 and 3 are day 16's two kernels with the three dimensions collapsed
// to one, because every case here is square. The N = 512 rows below and day
// 16's 512 rows are therefore two independent measurements of one pair, and
// they have to agree.
//
// Both cases divide by the tile on every axis, so the edge guard never fires.
// Day 16 owns the ragged size; holding it fixed here leaves the memory access
// pattern as the only variable in the table.
//
// The two sizes are picked against the T4's 4096 KiB L2. Three 512 x 512
// float matrices are 3 MiB and fit in it. Three 1024 x 1024 are 12 MiB and do
// not. That is the whole reason there are two cases rather than one.
//
// All three correctness launches happen before any timing, so the first three
// kernel launches of the process are one of each step. That is what makes
// `ncu --launch-count 3` collect the set this lesson quotes. The program
// prints the launch index of the second case so the second ncu command's
// --launch-skip is read off the run rather than counted by hand.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o matmul_steps \
// matmul_steps.cu
// Run: ./matmul_steps
//
#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)
// A 16 x 16 tile is one 256-thread block, the course default, and two tiles
// of floats is 2 KiB of shared memory against the 48 KiB a Turing block gets
// without an opt-in call. Day 44 is where the tile shape becomes the subject.
constexpr int kTileDim = 16;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
constexpr int kNumCases = 2;
constexpr int kNumKernels = 3;
constexpr size_t kCaseDim[kNumCases] = {512, 1024};
constexpr size_t kMaxDim = 1024;
// The T4's L2 is 4096 KiB, from FACT-SHEET.md section 4b. main() reads the
// real figure off the device and prints a calibration note when it differs,
// because the two case sizes were chosen against this number and against
// nothing else. On a card with a bigger L2 the second case stops being the
// out-of-cache case and the page says so.
constexpr size_t kReferenceL2Bytes = 4ull * 1024ull * 1024ull;
// Sectors one warp's global load request covers. Arithmetic from the block
// shape, not a measurement: the ncu report is what checks it.
//
// A 16 x 16 block linearises as threadIdx.x + 16 * threadIdx.y, so one warp
// is sixteen consecutive x and two consecutive y. Global memory is served in
// 32-byte sectors, which hold eight floats.
//
// Step 1, x is the row. A's sixteen rows are 4N bytes apart, so sixteen
// sectors, and B's two columns are adjacent floats inside one sector.
// (16 + 1) / 2 = 8.5.
// Step 2, x is the column. A's two rows are one sector each, and B's sixteen
// consecutive floats are 64 bytes, which is two sectors. (2 + 2) / 2 = 2.0.
// Step 3. Both tile loads are two rows of sixteen consecutive floats, four
// sectors each. (4 + 4) / 2 = 4.0. Higher than step 2, on one sixteenth as
// many requests, which is the point of the table.
constexpr double kSectorsPerRequest[kNumKernels] = {8.5, 2.0, 4.0};
constexpr const char* kKernelNames[kNumKernels] = {
"matmulStridedRows", "matmulCoalescedCols", "matmulTiled"};
static_assert(kTileDim == 16,
"kSectorsPerRequest above is derived for a 16 x 16 block. "
"Changing the tile means re-deriving all three figures, which "
"is day 44's job and not this file's");
static_assert((kTileDim * kTileDim) % 32 == 0,
"block size must be a whole number of warps");
static_assert(2 * kTileDim * kTileDim * sizeof(float) <= 48 * 1024,
"the two tiles must fit the 48 KiB a block gets on a T4 "
"without the cudaFuncSetAttribute opt-in");
static_assert(kCaseDim[0] % kTileDim == 0 && kCaseDim[1] % kTileDim == 0,
"both cases must divide by the tile, so the edge guard never "
"fires and the access pattern is the only variable");
static_assert(3 * kCaseDim[0] * kCaseDim[0] * sizeof(float) <=
kReferenceL2Bytes,
"case 0 exists because three of its matrices fit the T4's L2");
static_assert(3 * kCaseDim[1] * kCaseDim[1] * sizeof(float) > kReferenceL2Bytes,
"case 1 exists because three of its matrices do not");
static_assert(kCaseDim[1] == kMaxDim,
"every buffer is allocated once at kMaxDim x kMaxDim, so the "
"largest case has to be that size");
// Integer ceiling division. constexpr so one function sizes a grid at run
// time and appears in a static_assert. Anything that follows only from the
// constants above has to be a static_assert rather than an assert: CI builds
// Release, Release defines NDEBUG, and NDEBUG deletes assert().
static constexpr size_t ceilDiv(size_t a, size_t b) {
return (a + b - 1) / b;
}
static_assert(ceilDiv(kMaxDim, static_cast<size_t>(kTileDim)) <= 65535,
"gridDim.y stops at 65535, unlike gridDim.x which goes to "
"2^31 - 1");
// C = A * B, square N x N. One thread owns one element of C and reads a row
// of A and a column of B out of global memory: 2N loads for one store.
//
// Memory: threadIdx.x chooses the row here, so a warp's 32 lanes sit on
// sixteen consecutive rows of C and two consecutive columns. a[row * n + p]
// is sixteen addresses 4N bytes apart, one 32-byte sector each, and
// b[p * n + col] is two adjacent floats inside a single sector. Sixteen of
// the seventeen sectors that pair of requests pulls carry four useful bytes.
// The store spreads over sixteen sectors for the same reason.
//
// This kernel breaks the course's own index rule on purpose. CUDA-CODE-STYLE
// pins `row` to y and `col` to x, and that rule exists because of this exact
// kernel. Step 2 puts it back and the table prices the difference.
//
// Launch assumption: a grid that covers N on both axes. Because N is a
// multiple of the block dimension in every case here, the guard never fires,
// but it stays so the kernel is the same shape as day 16's. No barrier, which
// is the only reason a guard may wrap the whole body.
// snippet: strided-kernel
__global__ void matmulStridedRows(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ c, size_t n) {
const size_t row =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t col =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
if (row < n && col < n) {
float acc = 0.0f;
for (size_t p = 0; p < n; ++p) {
acc += a[row * n + p] * b[p * n + col];
}
c[row * n + col] = acc;
}
}
// end snippet
// The same product, the same 2N loads per thread, the same launch. Two lines
// differ from matmulStridedRows: threadIdx.x now chooses the column.
//
// Memory: a warp is two consecutive rows of C and sixteen consecutive
// columns. a[row * n + p] is two addresses, one sector each. b[p * n + col]
// is sixteen consecutive floats starting on a 64-byte boundary, which is two
// sectors, and both halves of the warp read the same two. Four sectors for
// the pair of requests against step 1's seventeen.
//
// Launch assumption: identical to step 1's, and the caller passes the same
// grid and block, because both cases are square. Nothing about the launch
// distinguishes these two kernels.
__global__ void matmulCoalescedCols(const float* __restrict__ a,
const float* __restrict__ b,
float* __restrict__ c, size_t n) {
// snippet: coalesced-index
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
// end snippet
if (row < n && col < n) {
float acc = 0.0f;
for (size_t p = 0; p < n; ++p) {
acc += a[row * n + p] * b[p * n + col];
}
c[row * n + col] = acc;
}
}
// Step 2 with both inputs staged through shared memory. One block owns one
// 16 x 16 tile of C and walks ceil(N / 16) tile steps. On each step every
// thread loads one element of A and one of B, the block synchronises, and
// every thread does sixteen multiply-adds out of shared memory. Global loads
// per thread fall from 2N to 2 * ceil(N / 16), a factor of the tile width.
//
// Memory: both global loads are two rows of sixteen consecutive floats, so
// four sectors per request, which is more per request than step 2 and one
// sixteenth as many requests. In shared memory a warp reads
// tileA[threadIdx.y][p], two addresses its halves share, and
// tileB[p][threadIdx.x], sixteen consecutive words read by two lanes each.
// Both are broadcasts or distinct banks, so neither conflicts; day 15 is
// where that stops being free.
//
// The ternary guards are dead code at these sizes, because both cases divide
// by the tile. They stay because removing them would make this a different
// kernel from day 16's and the 512 rows would stop being comparable.
//
// Launch assumption: exactly kTileDim x kTileDim threads per block. Every
// thread reaches both barriers; the guard covers the loads and the store and
// never the __syncthreads(), because a barrier only part of a block arrives
// at is undefined behaviour. Day 14 has the rule.
__global__ void matmulTiled(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ c,
size_t n) {
__shared__ float tileA[kTileDim][kTileDim];
__shared__ float tileB[kTileDim][kTileDim];
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
// ceilDiv is a host function, so the same arithmetic is written out here.
const size_t tiles = (n + kTileDim - 1) / kTileDim;
float acc = 0.0f;
for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
const size_t aCol = tileIdx * kTileDim + threadIdx.x;
const size_t bRow = tileIdx * kTileDim + threadIdx.y;
tileA[threadIdx.y][threadIdx.x] =
(row < n && aCol < n) ? a[row * n + aCol] : 0.0f;
tileB[threadIdx.y][threadIdx.x] =
(bRow < n && col < n) ? b[bRow * n + col] : 0.0f;
__syncthreads();
for (int p = 0; p < kTileDim; ++p) {
acc += tileA[threadIdx.y][p] * tileB[p][threadIdx.x];
}
__syncthreads();
}
if (row < n && col < n) {
c[row * n + col] = acc;
}
}
// One launch, chosen by step index, so the correctness pass and the timing
// pass cannot disagree about which kernel they ran. The grid and the block
// are the same for all three: the cases are square, so covering N rows and N
// columns is the same launch whichever axis a kernel reads them from.
static void launchStep(int step, const float* d_a, const float* d_b, float* d_c,
size_t n, dim3 grid, dim3 block) {
if (step == 0) {
matmulStridedRows<<<grid, block>>>(d_a, d_b, d_c, n);
} else if (step == 1) {
matmulCoalescedCols<<<grid, block>>>(d_a, d_b, d_c, n);
} else {
matmulTiled<<<grid, block>>>(d_a, d_b, d_c, n);
}
}
// CPU reference. Three nested loops in the obvious order, written for obvious
// correctness rather than speed: no blocking, no OpenMP, no intrinsics. It
// never allocates; the caller owns every buffer.
//
// It accumulates in double where the kernels accumulate in float, because a
// reference exists to be right rather than to match bit for bit. On these
// inputs both are exact anyway: makeInputs keeps every partial sum an integer
// well under 2^24.
//
// At N = 1024 this is 1.07 billion multiply-adds with poor locality on b, so
// it takes several seconds. That is the slowest part of the program and it is
// the price of checking every element of both cases.
static void matmulCpu(const float* a, const float* b, float* c, size_t n) {
for (size_t row = 0; row < n; ++row) {
for (size_t col = 0; col < n; ++col) {
double total = 0.0;
for (size_t p = 0; p < n; ++p) {
total += static_cast<double>(a[row * n + p]) *
static_cast<double>(b[p * n + col]);
}
c[row * n + col] = static_cast<float>(total);
}
}
}
// Fills both inputs with small integers held as floats: A in [-3, 3] and B in
// [-2, 2], the same generator day 16 uses. Every product is at most 6, so the
// largest dot product either case can produce is 6 * 1024 = 6,144, far below
// the 2^24 above which a float stops representing consecutive integers. Every
// intermediate on both processors is exact, so a mismatch in main() is an
// indexing bug and can be nothing else.
//
// The two periods, 7 and 5, divide neither 512 nor 1024, so a kernel that
// transposes an index reads a different value rather than the same one back.
static void makeInputs(std::vector<float>* h_a, std::vector<float>* h_b,
size_t elems) {
for (size_t e = 0; e < elems; ++e) {
(*h_a)[e] = static_cast<float>(static_cast<int>(e % 7) - 3);
(*h_b)[e] = static_cast<float>(static_cast<int>(e % 5) - 2);
}
}
// 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: main() turns it back into a row and a column, and "wrong
// at row 16, column 0" names the second block row. "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;
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// Copy it verbatim; the alternative is a second copy of the event
// boilerplate, which is how a warm-up goes missing from one of them.
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;
}
// 2 * N^3 flops for a square matrix multiply: one multiply and one add per
// term of every dot product. Computed here from the case dimension rather
// than per kernel, because all three kernels do exactly this much arithmetic
// and a second copy is a second thing to get wrong.
static double gflops(size_t n, float ms) {
const double work = 2.0 * static_cast<double>(n) * static_cast<double>(n) *
static_cast<double>(n);
return work / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Global loads a thread issues, from the constants. Steps 1 and 2 read a full
// row and a full column; step 3 reads two elements per tile step.
static size_t loadsPerThread(int step, size_t n) {
if (step == 2) {
return 2 * ceilDiv(n, static_cast<size_t>(kTileDim));
}
return 2 * n;
}
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, %d bytes of L2, %zu bytes of shared memory per "
"block\n",
prop.multiProcessorCount, prop.l2CacheSize, prop.sharedMemPerBlock);
std::printf(
"tile %d x %d, %d threads per block, %zu bytes of shared "
"memory per block\n\n",
kTileDim, kTileDim, kTileDim * kTileDim,
2 * sizeof(float) * kTileDim * kTileDim);
// The two case sizes only mean what the lesson says they mean while the
// card's L2 is the T4's. A different card gets a note rather than a
// silently wrong story.
if (static_cast<size_t>(prop.l2CacheSize) != kReferenceL2Bytes) {
std::printf(
"Note: the two case sizes were chosen against a %zu byte "
"L2 and this card has %d. Re-read the fits/does not fit "
"column below before comparing with the lesson.\n\n",
kReferenceL2Bytes, prop.l2CacheSize);
}
std::printf("Working set, three N x N float matrices against this L2\n");
for (int c = 0; c < kNumCases; ++c) {
const size_t working = 3 * kCaseDim[c] * kCaseDim[c] * sizeof(float);
std::printf(" N = %4zu %12zu bytes %s\n", kCaseDim[c], working,
(working <= static_cast<size_t>(prop.l2CacheSize))
? "fits L2"
: "does not fit L2");
}
// Both blocks below are arithmetic from the constants at the top of this
// file, not measurements. They are printed so the ncu report has
// something to be checked against.
std::printf(
"\nSectors per global load request, from the 16 x 16 block "
"shape\n");
for (int s = 0; s < kNumKernels; ++s) {
std::printf(" %-20s %5.1f\n", kKernelNames[s], kSectorsPerRequest[s]);
}
const int launchesPerCase = kNumKernels * (1 + kWarmupRuns + kTimedRuns);
std::printf(
"\nProfiler: launches 1 to %d are one of each kernel at N = "
"%zu.\n Launches %d to %d are one of each at N = "
"%zu.\n",
kNumKernels, kCaseDim[0], launchesPerCase + 1,
launchesPerCase + kNumKernels, kCaseDim[1]);
std::printf(
"\nTiming: mean of %d runs after %d warm-ups, CUDA events, "
"copies not included\n\n",
kTimedRuns, kWarmupRuns);
const size_t maxElems = kMaxDim * kMaxDim;
const size_t maxBytes = maxElems * sizeof(float);
std::vector<float> h_a(maxElems);
std::vector<float> h_b(maxElems);
std::vector<float> h_c(maxElems);
std::vector<float> h_want(maxElems);
float* d_a = nullptr;
float* d_b = nullptr;
float* d_c = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, maxBytes));
CUDA_CHECK(cudaMalloc(&d_b, maxBytes));
CUDA_CHECK(cudaMalloc(&d_c, maxBytes));
std::printf("%4s %-19s %11s %10s %10s %8s\n", "N", "kernel",
"time (ms)", "GFLOP/s", "loads/thr", "sect/req");
std::printf("%4s %-19s %11s %10s %10s %8s\n", "----",
"-------------------", "-----------", "----------",
"----------", "--------");
// Failure is recorded rather than returned, so every path falls through
// to the one cleanup block below and no cudaMalloc escapes its cudaFree.
int rows = 0;
int badCase = kNumCases;
int badStep = 0;
size_t badIndex = 0;
int badTimingCase = kNumCases;
float ms[kNumKernels] = {0.0f, 0.0f, 0.0f};
float badMs[kNumKernels] = {0.0f, 0.0f, 0.0f};
for (int c = 0; c < kNumCases && badCase == kNumCases; ++c) {
const size_t n = kCaseDim[c];
const size_t elems = n * n;
const size_t bytes = elems * sizeof(float);
makeInputs(&h_a, &h_b, elems);
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
matmulCpu(h_a.data(), h_b.data(), h_want.data(), n);
const dim3 block(kTileDim, kTileDim);
const dim3 grid(static_cast<unsigned int>(
ceilDiv(n, static_cast<size_t>(kTileDim))),
static_cast<unsigned int>(
ceilDiv(n, static_cast<size_t>(kTileDim))));
// Every kernel is checked before any kernel is timed. That ordering
// is what makes the first three launches of the process one of each
// step, which is what the ncu commands in README.md rely on.
for (int s = 0; s < kNumKernels; ++s) {
// The output is cleared first so a kernel that skips elements is
// caught by the comparison rather than by the previous step's
// answer agreeing with it.
CUDA_CHECK(cudaMemset(d_c, 0, bytes));
launchStep(s, d_a, d_b, d_c, n, grid, block);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_c.data(), d_c, bytes, cudaMemcpyDeviceToHost));
const size_t bad =
firstMismatch(h_c.data(), h_want.data(), elems, kRelTolerance);
if (bad != elems) {
badCase = c;
badStep = s;
badIndex = bad;
break;
}
}
if (badCase != kNumCases) {
break;
}
for (int s = 0; s < kNumKernels; ++s) {
ms[s] = timeKernel(
[&] { launchStep(s, d_a, d_b, d_c, n, grid, block); });
std::printf("%4zu %-19s %11.3f %10.1f %10zu %8.1f\n", n,
kKernelNames[s], ms[s], gflops(n, ms[s]),
loadsPerThread(s, n), kSectorsPerRequest[s]);
++rows;
}
// A zero or negative interval means the event pair never separated,
// and every ratio derived from it would be nonsense. Recorded rather
// than returned, so this path also falls through to the one cleanup
// block below.
if (ms[0] <= 0.0f || ms[1] <= 0.0f || ms[2] <= 0.0f) {
badTimingCase = c;
for (int s = 0; s < kNumKernels; ++s) {
badMs[s] = ms[s];
}
break;
}
std::printf(
" speedup step1->step2 %5.2fx step2->step3 %5.2fx"
" step1->step3 %5.2fx\n",
static_cast<double>(ms[0]) / ms[1],
static_cast<double>(ms[1]) / ms[2],
static_cast<double>(ms[0]) / ms[2]);
}
CUDA_CHECK(cudaFree(d_a));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_c));
if (badCase != kNumCases) {
const size_t n = kCaseDim[badCase];
std::fprintf(stderr,
"%s wrong at %zu (row %zu, col %zu) on the %zu x %zu "
"case\n",
kKernelNames[badStep], badIndex, badIndex / n,
badIndex % n, n, n);
return EXIT_FAILURE;
}
if (badTimingCase != kNumCases) {
std::fprintf(stderr,
"a timed interval was not positive at N = %zu: "
"%.6f %.6f %.6f ms\n",
kCaseDim[badTimingCase], badMs[0], badMs[1], badMs[2]);
return EXIT_FAILURE;
}
// 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. It exists so the lesson's table and
// this program cannot drift apart on how many rows there are.
if (rows != kNumCases * kNumKernels) {
std::fprintf(stderr,
"printed %d rows, expected %d; the lesson's table and "
"this program disagree\n",
rows, kNumCases * kNumKernels);
return EXIT_FAILURE;
}
std::printf(
"\nall %d kernels match the CPU reference at every element of "
"both cases\n",
kNumKernels);
return EXIT_SUCCESS;
}