code/day81-cublas/cublas_lt.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 81: one FP32 GEMM plus a bias, four ways, so the cost of the bias is
// a row in a table instead of an argument.
//
// gemmEx + bias kernel cublasGemmEx, then addBias as its own launch
// Lt + bias kernel cublasLtMatmul with the default epilogue, then
// the same addBias launch
// Lt bias epilogue cublasLtMatmul with CUBLASLT_EPILOGUE_BIAS, one
// launch, no second pass over C
// bias kernel alone addBias with no GEMM, so the fused path's saving
// has a measured price tag next to it
//
// All four paths do the same arithmetic: CUDA_R_32F in and out,
// CUBLAS_COMPUTE_32F, math mode pinned to CUBLAS_DEFAULT_MATH. That is the
// comparison day 44 fixed for this course, and it is the only way the
// percentages on the page mean anything.
//
// The sizes are day 44's, 512, 1024 and 2048, so the %cuBLAS column here and
// the one on that page describe the same three problems on the same card.
//
// Column major: cuBLAS reads every matrix column major, the host buffers here
// are row major, and no buffer is transposed anywhere in this file. The whole
// trick is in which pointer goes first. See the comment above callGemmEx.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o cublas_lt cublas_lt.cu \
// -lcublas -lcublasLt
// Run: ./cublas_lt
//
// VERIFIED: Tesla T4, driver 580.173.02, CUDA 12.6, 2026-09-02.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cublasLt.h>
#include <cublas_v2.h>
#include <cuda_runtime.h>
// The one error macro, byte identical to day 5's.
#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)
// cuBLAS and cuBLASLt return cublasStatus_t, which CUDA_CHECK cannot take.
// Keep that conversion in a function instead of adding another error macro.
// snippet: cublas-check
static void checkCublas(cublasStatus_t status) {
if (status != CUBLAS_STATUS_SUCCESS) {
std::fprintf(stderr, "cuBLAS error: %s\n",
cublasGetStatusString(status));
std::exit(EXIT_FAILURE);
}
}
// end snippet
constexpr int kCaseCount = 3;
constexpr int kCases[kCaseCount] = {512, 1024, 2048};
constexpr int kThreadsPerBlock = 256;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// How many candidate algorithms to ask the heuristic for. The library
// returns them in increasing estimated compute time, so index 0 is its
// answer and the rest are the search space it ranked below that.
constexpr int kMaxAlgos = 8;
// cuBLASLt takes the workspace as an argument rather than owning one. The
// documentation's floor is 16 KiB to avoid CUBLAS_STATUS_ALLOC_FAILED; 32
// MiB is the size the cuBLASLt sample uses and it is small next to the 48
// MiB of matrices at the largest case.
constexpr size_t kWorkspaceBytes = 32u << 20;
constexpr int kPathCount = 4;
static const char* const kPathNames[kPathCount] = {
"gemmEx + bias kernel", "Lt + bias kernel", "Lt bias epilogue",
"bias kernel alone"};
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kWorkspaceBytes >= 16u << 10,
"the cuBLAS documentation's floor for avoiding "
"CUBLAS_STATUS_ALLOC_FAILED is 16 KiB");
// out[row][col] = in[row][col] + bias[col], row major, out of place.
//
// One thread owns one element. A warp's 32 threads take 32 consecutive
// columns of one row, so the read of `in`, the write of `out` and the read of
// `bias` are each 32 consecutive floats: four 32-byte sectors per warp per
// access, which is day 11's coalesced case.
//
// Out of place on purpose. Timing an in-place update in a loop of ten runs
// would either accumulate the bias ten times or need a reset inside the
// timed region, and the byte count is the same either way: one read of C and
// one write of C.
//
// Launch assumption: gridDim.x * blockDim.x >= size, gridDim.y == size.
// snippet: add-bias
__global__ void addBias(const float* __restrict__ in,
const float* __restrict__ bias, float* __restrict__ out,
int size) {
const int col = blockIdx.x * blockDim.x + threadIdx.x;
const int row = blockIdx.y;
if (col < size) {
const size_t elem = static_cast<size_t>(row) * size + col;
out[elem] = in[elem] + bias[col];
}
}
// end snippet
static void launchAddBias(const float* d_in, const float* d_bias, float* d_out,
int size, cudaStream_t stream) {
const int blocksX = (size + kThreadsPerBlock - 1) / kThreadsPerBlock;
const dim3 grid(static_cast<unsigned int>(blocksX),
static_cast<unsigned int>(size));
addBias<<<grid, kThreadsPerBlock, 0, stream>>>(d_in, d_bias, d_out, size);
}
// Inputs are deliberately not representable in binary floating point. Small
// integers held as floats would make every partial sum exact on both
// processors, the comparison below would collapse to equality, and the
// tolerance would be decoration. These values are 0.25 + 0.1 * k, and 0.1 is
// a repeating binary fraction, so the GPU's summation order and the host's
// genuinely disagree in the last bits and the bound has to hold.
//
// Every value is positive, so no partial sum cancels. That is what makes a
// relative tolerance meaningful: |expected| is the same size as the sum of
// the magnitudes that produced it.
static void fillInputs(float* a, float* b, float* bias, int size) {
const size_t n = static_cast<size_t>(size) * size;
for (size_t i = 0; i < n; ++i) {
a[i] = 0.25f + 0.1f * static_cast<float>((i * 31u + 7u) % 9u);
b[i] = 0.25f + 0.1f * static_cast<float>((i * 17u + 3u) % 11u);
}
for (int j = 0; j < size; ++j) {
bias[j] = -1.0f + 0.1f * static_cast<float>(j % 21);
}
}
// C = A * B + bias, row major, accumulated in double.
//
// No Kahan summation here, unlike the harness references in
// code/reference/. Plain double accumulation over K = 2048 terms carries a
// relative error near K * 2^-53, about 2e-13, which is seven orders of
// magnitude under the f32 bound this reference is used to enforce. Kahan
// would cost another pass over 8.6 billion terms and buy nothing.
static void matmulBiasCpu(const float* a, const float* b, const float* bias,
double* out, int size) {
for (int row = 0; row < size; ++row) {
for (int col = 0; col < size; ++col) {
double sum = 0.0;
for (int k = 0; k < size; ++k) {
sum +=
static_cast<double>(
a[row * static_cast<size_t>(size) + k]) *
static_cast<double>(b[k * static_cast<size_t>(size) + col]);
}
out[row * static_cast<size_t>(size) + col] =
sum + static_cast<double>(bias[col]);
}
}
}
// The tolerance from EXERCISE-DESIGN.md, scaled by the reduction depth:
//
// rtol_K = max(1e-5, 4 * 2^-23 * sqrt(K))
// atol_K = rtol_K * max|expected|
//
// sqrt(K) and not K because the rounding errors in a blocked GEMM are not
// correlated; the worst-case linear bound would pass a genuinely broken
// kernel. The 4 is slack for a different but valid summation order, and a
// split-K GEMM is exactly that. At K = 2048 this is 2.16e-5, about 180 ulp
// at 1.0, and the program prints the arithmetic next to the number.
static double relTolerance(int size) {
const double eps = 1.0 / 8388608.0; // 2^-23, the f32 spacing at 1.0
const double scaled = 4.0 * eps * std::sqrt(static_cast<double>(size));
return (scaled > 1e-5) ? scaled : 1e-5;
}
// Returns the index of the first element outside the tolerance, or n.
static size_t firstMismatch(const float* got, const double* want, size_t n,
double rtol, double atol) {
for (size_t i = 0; i < n; ++i) {
const double diff = std::fabs(static_cast<double>(got[i]) - want[i]);
if (!(diff <= atol + rtol * std::fabs(want[i]))) {
return i;
}
}
return n;
}
static double maxAbs(const double* v, size_t n) {
double m = 0.0;
for (size_t i = 0; i < n; ++i) {
if (std::fabs(v[i]) > m) {
m = std::fabs(v[i]);
}
}
return m;
}
// Day 9's timeKernel with one change, and the change is the whole reason the
// stream exists: the events are recorded on the stream the work is queued
// on, not on the legacy default stream. Everything this file launches, the
// library calls included, runs on `stream`.
template <typename LaunchFn>
static float timeOnStream(LaunchFn launch, cudaStream_t stream) {
cudaEvent_t start;
cudaEvent_t stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this path, not just the first path in the program. Lazy module
// loading has been the default since CUDA 12.2 on Linux, and both
// libraries pick and cache a kernel on the first call, so an unwarmed
// comparison flatters whichever path ran second.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaStreamSynchronize(stream));
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start, stream));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop, stream));
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;
}
// C = A * B with row-major A, B and C, through a column-major library, with
// nothing transposed and nothing copied.
//
// Read a row-major size x size buffer as column major and you are holding
// its transpose. So the buffers on the device are A^T, B^T and C^T as far as
// cuBLAS is concerned, and (A*B)^T = B^T * A^T. Passing d_b first and d_a
// second with both operations CUBLAS_OP_N asks for exactly that product, and
// the answer it writes, C^T in column major, is C in row major. The swap is
// free; a real transpose would cost two passes over the data.
//
// snippet: gemm-ex-call
static cublasStatus_t callGemmEx(cublasHandle_t handle, const float* d_a,
const float* d_b, float* d_c, int size) {
const float alpha = 1.0f;
const float beta = 0.0f;
return cublasGemmEx(handle, CUBLAS_OP_N, CUBLAS_OP_N, size, size, size,
&alpha, d_b, CUDA_R_32F, size, d_a, CUDA_R_32F, size,
&beta, d_c, CUDA_R_32F, size, CUBLAS_COMPUTE_32F,
CUBLAS_GEMM_DEFAULT);
}
// end snippet
// Everything cuBLASLt needs to describe one problem. cuBLAS packs all of
// this into one call and one handle; Lt makes you say it, which is what buys
// the epilogue and the algorithm list.
struct LtProblem {
cublasLtMatmulDesc_t desc;
cublasLtMatrixLayout_t layoutA; // holds d_b, per the swap above
cublasLtMatrixLayout_t layoutB; // holds d_a
cublasLtMatrixLayout_t layoutD; // C, and C == D with beta = 0
cublasLtMatmulPreference_t pref;
};
// Builds the descriptors for one case at one epilogue setting.
//
// The bias vector's length "must match matrix D rows" and it is "broadcast to
// all columns" (cuBLAS documentation, cublasLtEpilogue_t). Under the swap
// above, D's rows are the row-major C's columns, so a bias per output column,
// which is the bias every linear layer in every network has, lands correctly
// with no reshaping. The convention that trips everyone on the way in pays
// for itself here.
static cublasStatus_t buildLtProblem(LtProblem* p, int size,
cublasLtEpilogue_t epilogue,
const float* d_bias) {
const cublasOperation_t opN = CUBLAS_OP_N;
const size_t workspace = kWorkspaceBytes;
cublasStatus_t st =
cublasLtMatmulDescCreate(&p->desc, CUBLAS_COMPUTE_32F, CUDA_R_32F);
if (st != CUBLAS_STATUS_SUCCESS) {
return st;
}
st = cublasLtMatmulDescSetAttribute(p->desc, CUBLASLT_MATMUL_DESC_TRANSA,
&opN, sizeof(opN));
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatmulDescSetAttribute(
p->desc, CUBLASLT_MATMUL_DESC_TRANSB, &opN, sizeof(opN));
}
// snippet: lt-epilogue-attrs
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatmulDescSetAttribute(p->desc,
CUBLASLT_MATMUL_DESC_EPILOGUE,
&epilogue, sizeof(epilogue));
}
if (st == CUBLAS_STATUS_SUCCESS && epilogue == CUBLASLT_EPILOGUE_BIAS) {
st = cublasLtMatmulDescSetAttribute(p->desc,
CUBLASLT_MATMUL_DESC_BIAS_POINTER,
&d_bias, sizeof(d_bias));
}
// end snippet
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatrixLayoutCreate(&p->layoutA, CUDA_R_32F, size, size,
size);
}
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatrixLayoutCreate(&p->layoutB, CUDA_R_32F, size, size,
size);
}
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatrixLayoutCreate(&p->layoutD, CUDA_R_32F, size, size,
size);
}
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatmulPreferenceCreate(&p->pref);
}
if (st == CUBLAS_STATUS_SUCCESS) {
st = cublasLtMatmulPreferenceSetAttribute(
p->pref, CUBLASLT_MATMUL_PREF_MAX_WORKSPACE_BYTES, &workspace,
sizeof(workspace));
}
return st;
}
static void destroyLtProblem(LtProblem* p) {
if (p->pref != nullptr) {
checkCublas(cublasLtMatmulPreferenceDestroy(p->pref));
p->pref = nullptr;
}
if (p->layoutD != nullptr) {
checkCublas(cublasLtMatrixLayoutDestroy(p->layoutD));
p->layoutD = nullptr;
}
if (p->layoutB != nullptr) {
checkCublas(cublasLtMatrixLayoutDestroy(p->layoutB));
p->layoutB = nullptr;
}
if (p->layoutA != nullptr) {
checkCublas(cublasLtMatrixLayoutDestroy(p->layoutA));
p->layoutA = nullptr;
}
if (p->desc != nullptr) {
checkCublas(cublasLtMatmulDescDestroy(p->desc));
p->desc = nullptr;
}
}
// One cublasLtMatmul call. Same operand swap as callGemmEx: d_b is the
// library's A. beta is zero, so C is never read and C and D can be the same
// buffer with the same layout; the documentation asks for zero beta and host
// pointer mode for the best bias-epilogue path, and host pointer mode is the
// default.
// snippet: lt-matmul-call
static cublasStatus_t callLtMatmul(cublasLtHandle_t lt, const LtProblem* p,
const cublasLtMatmulAlgo_t* algo,
const float* d_a, const float* d_b,
float* d_out, void* workspace,
cudaStream_t stream) {
const float alpha = 1.0f;
const float beta = 0.0f;
return cublasLtMatmul(lt, p->desc, &alpha, d_b, p->layoutA, d_a, p->layoutB,
&beta, d_out, p->layoutD, d_out, p->layoutD, algo,
workspace, kWorkspaceBytes, stream);
}
// end snippet
// Copies the device result back and compares it to the CPU reference.
// Returns EXIT_FAILURE and names the first bad element, so a broken path
// cannot post a time.
static int checkPath(const char* name, const float* d_out, float* h_out,
const double* want, int size, double rtol, double atol) {
const size_t n = static_cast<size_t>(size) * size;
CUDA_CHECK(
cudaMemcpy(h_out, d_out, n * sizeof(float), cudaMemcpyDeviceToHost));
const size_t bad = firstMismatch(h_out, want, n, rtol, atol);
if (bad != n) {
std::fprintf(stderr,
"%s wrong at size %d, row %zu col %zu: got %.9g, "
"want %.9g, allowed %.3g\n",
name, size, bad / static_cast<size_t>(size),
bad % static_cast<size_t>(size),
static_cast<double>(h_out[bad]), want[bad],
atol + rtol * std::fabs(want[bad]));
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
// Runs and times every path at one size, plus the algorithm list.
static int runOneCase(cublasHandle_t handle, cublasLtHandle_t lt, int size,
const float* d_a, const float* d_b, const float* d_bias,
float* d_c, float* d_out, void* d_workspace,
const float* h_a, const float* h_b, const float* h_bias,
float* h_out, double* h_ref, cudaStream_t stream) {
const size_t n = static_cast<size_t>(size) * size;
const double gflop = 2.0 * static_cast<double>(size) * size * size / 1.0e9;
const double biasBytes = 2.0 * static_cast<double>(n) * sizeof(float);
matmulBiasCpu(h_a, h_b, h_bias, h_ref, size);
const double rtol = relTolerance(size);
const double atol = rtol * maxAbs(h_ref, n);
LtProblem plain = {nullptr, nullptr, nullptr, nullptr, nullptr};
LtProblem fused = {nullptr, nullptr, nullptr, nullptr, nullptr};
int status = EXIT_SUCCESS;
cublasStatus_t st =
buildLtProblem(&plain, size, CUBLASLT_EPILOGUE_DEFAULT, nullptr);
if (st != CUBLAS_STATUS_SUCCESS) {
std::fprintf(stderr, "plain Lt descriptors failed: %s\n",
cublasGetStatusString(st));
status = EXIT_FAILURE;
}
if (status == EXIT_SUCCESS) {
st = buildLtProblem(&fused, size, CUBLASLT_EPILOGUE_BIAS, d_bias);
if (st != CUBLAS_STATUS_SUCCESS) {
std::fprintf(stderr, "bias-epilogue descriptors failed: %s\n",
cublasGetStatusString(st));
status = EXIT_FAILURE;
}
}
// The heuristic, for the fused problem. The library returns candidates in
// increasing estimated compute time, so results[0] is its answer. A
// return of zero candidates is the honest failure mode for an epilogue
// this card or this toolkit does not implement, and it is a refusal, not
// a crash.
cublasLtMatmulHeuristicResult_t results[kMaxAlgos] = {};
int algoCount = 0;
if (status == EXIT_SUCCESS) {
st = cublasLtMatmulAlgoGetHeuristic(
lt, fused.desc, fused.layoutA, fused.layoutB, fused.layoutD,
fused.layoutD, fused.pref, kMaxAlgos, results, &algoCount);
if (st != CUBLAS_STATUS_SUCCESS || algoCount == 0) {
std::fprintf(stderr,
"no bias-epilogue algorithm at size %d: %s, "
"%d candidates\n",
size, cublasGetStatusString(st), algoCount);
status = EXIT_FAILURE;
}
}
cublasLtMatmulHeuristicResult_t plainResults[kMaxAlgos] = {};
int plainCount = 0;
if (status == EXIT_SUCCESS) {
st = cublasLtMatmulAlgoGetHeuristic(
lt, plain.desc, plain.layoutA, plain.layoutB, plain.layoutD,
plain.layoutD, plain.pref, kMaxAlgos, plainResults, &plainCount);
if (st != CUBLAS_STATUS_SUCCESS || plainCount == 0) {
std::fprintf(stderr,
"no default-epilogue algorithm at size %d: %s, "
"%d candidates\n",
size, cublasGetStatusString(st), plainCount);
status = EXIT_FAILURE;
}
}
// Correctness for all four paths before any of them is timed.
if (status == EXIT_SUCCESS) {
checkCublas(callGemmEx(handle, d_a, d_b, d_c, size));
launchAddBias(d_c, d_bias, d_out, size, stream);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaStreamSynchronize(stream));
status =
checkPath(kPathNames[0], d_out, h_out, h_ref, size, rtol, atol);
}
if (status == EXIT_SUCCESS) {
checkCublas(callLtMatmul(lt, &plain, &plainResults[0].algo, d_a, d_b,
d_c, d_workspace, stream));
launchAddBias(d_c, d_bias, d_out, size, stream);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaStreamSynchronize(stream));
status =
checkPath(kPathNames[1], d_out, h_out, h_ref, size, rtol, atol);
}
if (status == EXIT_SUCCESS) {
// Zero the output first, so a fused path that writes nothing is
// caught by the comparison instead of being rescued by the answer
// the previous path left in the same memory.
CUDA_CHECK(cudaMemsetAsync(d_out, 0, n * sizeof(float), stream));
checkCublas(callLtMatmul(lt, &fused, &results[0].algo, d_a, d_b, d_out,
d_workspace, stream));
CUDA_CHECK(cudaStreamSynchronize(stream));
status =
checkPath(kPathNames[2], d_out, h_out, h_ref, size, rtol, atol);
}
double ms[kPathCount] = {0.0, 0.0, 0.0, 0.0};
if (status == EXIT_SUCCESS) {
ms[0] = static_cast<double>(timeOnStream(
[&] {
checkCublas(callGemmEx(handle, d_a, d_b, d_c, size));
launchAddBias(d_c, d_bias, d_out, size, stream);
},
stream));
ms[1] = static_cast<double>(timeOnStream(
[&] {
checkCublas(callLtMatmul(lt, &plain, &plainResults[0].algo,
d_a, d_b, d_c, d_workspace, stream));
launchAddBias(d_c, d_bias, d_out, size, stream);
},
stream));
ms[2] = static_cast<double>(timeOnStream(
[&] {
checkCublas(callLtMatmul(lt, &fused, &results[0].algo, d_a,
d_b, d_out, d_workspace, stream));
},
stream));
ms[3] = static_cast<double>(timeOnStream(
[&] { launchAddBias(d_c, d_bias, d_out, size, stream); }, stream));
}
if (status == EXIT_SUCCESS) {
std::printf(
"\nsize %d, %.3f GFLOP per GEMM, mean of %d runs after "
"%d warm-ups\n",
size, gflop, kTimedRuns, kWarmupRuns);
std::printf("tolerance |got-exp| <= %.3g + %.3g*|exp|\n", atol, rtol);
std::printf(
" reduction over K=%d, "
"rtol = 4*2^-23*sqrt(K) = %.3g\n",
size, rtol);
std::printf("path ms GFLOP/s vs fused\n");
std::printf("--------------------- --------- --------- ----------\n");
for (int p = 0; p < kPathCount; ++p) {
if (p == kPathCount - 1) {
// No GFLOP/s for the bias row. It does none of the GEMM's
// arithmetic, so that column would divide the matmul's
// flops by a time that never computed them.
std::printf("%-21s %9.4f %9s %10.3f\n", kPathNames[p], ms[p],
"-", ms[p] / ms[2]);
} else {
std::printf("%-21s %9.4f %9.1f %10.3f\n", kPathNames[p], ms[p],
gflop * 1.0e3 / ms[p], ms[p] / ms[2]);
}
}
std::printf("the bias kernel alone moves %.1f MB at %.1f GB/s\n",
biasBytes / 1.0e6, biasBytes / (ms[3] * 1.0e6));
}
// The algorithm list, which is the other half of what Lt sells. wavesCount
// is the library's own device-utilization estimate: 1.0 means the launch
// fills the GPU exactly once.
if (status == EXIT_SUCCESS) {
std::printf(
"\nbias-epilogue candidates at size %d, "
"heuristic order\n",
size);
std::printf("index ms waves workspace B\n");
std::printf("----- ------- --------- ------------\n");
for (int i = 0; i < algoCount; ++i) {
if (results[i].state != CUBLAS_STATUS_SUCCESS) {
std::printf("%5d - - - (%s)\n", i,
cublasGetStatusString(results[i].state));
continue;
}
const double algoMs = static_cast<double>(timeOnStream(
[&] {
checkCublas(callLtMatmul(lt, &fused, &results[i].algo, d_a,
d_b, d_out, d_workspace, stream));
},
stream));
std::printf("%5d %7.4f %9.3f %12zu\n", i, algoMs,
static_cast<double>(results[i].wavesCount),
results[i].workspaceSize);
}
}
destroyLtProblem(&fused);
destroyLtProblem(&plain);
return status;
}
int main() {
int device = 0;
CUDA_CHECK(cudaGetDevice(&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);
const int maxSize = kCases[kCaseCount - 1];
const size_t maxElems = static_cast<size_t>(maxSize) * maxSize;
const size_t maxBytes = maxElems * sizeof(float);
std::vector<float> h_a(maxElems);
std::vector<float> h_b(maxElems);
std::vector<float> h_bias(static_cast<size_t>(maxSize));
std::vector<float> h_out(maxElems);
std::vector<double> h_ref(maxElems);
float* d_a = nullptr;
float* d_b = nullptr;
float* d_bias = nullptr;
float* d_c = nullptr;
float* d_out = nullptr;
void* d_workspace = nullptr;
CUDA_CHECK(cudaMalloc(&d_a, maxBytes));
CUDA_CHECK(cudaMalloc(&d_b, maxBytes));
CUDA_CHECK(
cudaMalloc(&d_bias, static_cast<size_t>(maxSize) * sizeof(float)));
CUDA_CHECK(cudaMalloc(&d_c, maxBytes));
CUDA_CHECK(cudaMalloc(&d_out, maxBytes));
CUDA_CHECK(cudaMalloc(&d_workspace, kWorkspaceBytes));
cudaStream_t stream = nullptr;
CUDA_CHECK(cudaStreamCreate(&stream));
cublasHandle_t handle = nullptr;
checkCublas(cublasCreate(&handle));
// The handle carries the stream, and every subsequent call on this
// handle runs on it. cuBLASLt takes the stream as an argument instead,
// which is why callLtMatmul has one and callGemmEx does not.
//
// Order matters: cublasSetStream "unconditionally resets the cuBLAS
// library workspace back to the default workspace pool", so a
// cublasSetWorkspace call has to come after this one, not before.
checkCublas(cublasSetStream(handle, stream));
// FP32 all the way through, stated in code rather than only in prose. On
// a T4 there is no TF32 path to fall into; on an Ampere or newer card
// this line is the difference between a fair comparison and a rout.
checkCublas(cublasSetMathMode(handle, CUBLAS_DEFAULT_MATH));
// A cuBLAS handle "encapsulates a cuBLASLt handle. Any valid
// cublasHandle_t can be used in place of cublasLtHandle_t with a simple
// cast" (cuBLAS documentation, cublasLtHandle_t). One context, two APIs.
cublasLtHandle_t lt = reinterpret_cast<cublasLtHandle_t>(handle);
// Printed, not asserted: the library that answers is the one that was
// linked, and day 44 shipped a binary linked against the wrong major
// by accident.
int version = 0;
checkCublas(cublasGetVersion(handle, &version));
std::printf("cuBLAS version %d, cuBLASLt version %zu\n", version,
cublasLtGetVersion());
std::printf("cuBLASLt workspace: %zu bytes\n", kWorkspaceBytes);
int status = EXIT_SUCCESS;
for (int c = 0; c < kCaseCount && status == EXIT_SUCCESS; ++c) {
const int size = kCases[c];
const size_t elems = static_cast<size_t>(size) * size;
fillInputs(h_a.data(), h_b.data(), h_bias.data(), size);
CUDA_CHECK(cudaMemcpy(d_a, h_a.data(), elems * sizeof(float),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), elems * sizeof(float),
cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_bias, h_bias.data(),
static_cast<size_t>(size) * sizeof(float),
cudaMemcpyHostToDevice));
status = runOneCase(handle, lt, size, d_a, d_b, d_bias, d_c, d_out,
d_workspace, h_a.data(), h_b.data(), h_bias.data(),
h_out.data(), h_ref.data(), stream);
}
if (status == EXIT_SUCCESS) {
std::printf("\nall paths match the CPU reference at every size\n");
}
checkCublas(cublasDestroy(handle));
CUDA_CHECK(cudaStreamDestroy(stream));
CUDA_CHECK(cudaFree(d_workspace));
CUDA_CHECK(cudaFree(d_out));
CUDA_CHECK(cudaFree(d_c));
CUDA_CHECK(cudaFree(d_bias));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_a));
return status;
}