100 DAYS OF GPU Course preview
Day 81Module 9
in-technical-review

cuBLAS and cuBLASLt

Day 44 ended with a hand-written float4 matmul at 67.7 percent of cuBLAS's cublasGemmEx on a Tesla T4. cuBLAS is a set of kernels plus an algorithm chooser, not one kernel.

This lesson handles row-major matrices by swapping two arguments, with no copy. It then compares cublasGemmEx with cublasLtMatmul and measures a fused bias against the bias kernel it replaces.

What the handle owns, and what it reads column major

cublasCreate gives you a handle, and the handle holds everything that is not an argument: the device, the math mode, the workspace, and the stream every subsequent call runs on. cublasSetStream(handle, stream) is how this library's work joins the rest of your pipeline, and it has an edge worth knowing before you hit it. The call "unconditionally resets the cuBLAS library workspace back to the default workspace pool" (https://docs.nvidia.com/cuda/cublas/index.html#cublassetstream , checked 2026-09-01), so a cublasSetWorkspace written above it is silently undone.

cuBLAS "uses column-major storage, and 1-based indexing" (https://docs.nvidia.com/cuda/cublas/index.html#data-layout , checked 2026-09-01) and your C++ arrays are row major. The usual first reaction is to transpose the inputs, two passes over global memory before any arithmetic starts.

Do not transpose the buffers. A row-major buffer read as column major is its transpose, and matrix algebra gives the identity (A*B)^T = B^T * A^T. So hand the library B first and A second, leave both operations at CUBLAS_OP_N, and ask for C^T.

The answer it writes, column major, is the answer you wanted, row major, in the buffer you already had.

Diagram: one buffer, two readings. Left: a 3 by 4 row-major array drawn as a flat run of 12 boxes with the row-major reading traced over it, and the same run with the column-major reading traced over it, arriving at a 4 by 3 grid. Right: the call, with d_b entering the A slot and d_a entering the B slot, and the output labelled C^T on the library's side and C on yours. Caption on the left pair: "same 12 floats, two shapes." Caption on the right: "two arguments swapped, zero bytes copied." Alt text: A row-major array read column major is its own transpose, so passing B before A and asking for the transposed product returns the row-major answer in place. Twelve floats, two readings, no copy.

Lt is not a newer spelling of cublasGemmEx

cuBLASLt's cublasLtMatmul adds a descriptor, three layouts, a preference object, caller-owned workspace, and an explicit algorithm query. Those objects expose options that cublasGemmEx cannot.

The extra objects are the product. Two things you cannot express in cublasGemmEx become arguments in cublasLtMatmul. The first is the epilogue: CUBLASLT_EPILOGUE_BIAS adds a bias vector to the product inside the same kernel, and the ReLU and GELU variants apply an activation there too, so the pointwise pass you were going to launch afterwards never runs.

That is day 48's fusion argument, handed to you by the library instead of written by you. The second is the algorithm: cublasLtMatmulAlgoGetHeuristic hands back a ranked list of candidate kernels with the workspace and the wave count each one wants, and you pick. cublasGemmEx makes that choice for you and does not say what it rejected.

The epilogue also pays back the previous section. The bias vector's length "must match matrix D rows" and it "is broadcast to all columns" (https://docs.nvidia.com/cuda/cublas/index.html#cublasltepilogue-t , checked 2026-09-01). After the swap, the library's D rows are your row-major output's columns, so one bias per output feature, which is what every linear layer wants, lands with no reshaping at all.

Four paths, and the bias with its own price tag

Full program in cublas_lt.cu, with the install recipe in the README beside it.

The four paths are one GEMM plus one bias, arranged so no two of them differ in more than one thing. gemmEx + bias kernel calls cublasGemmEx and then launches addBias. Lt + bias kernel swaps in cublasLtMatmul with the default epilogue and keeps the same addBias, so its gap to the row above is the price of the interface and nothing else.

Lt bias epilogue moves the bias inside. bias kernel alone runs addBias on its own, which turns the third row's saving from a claim into a subtraction.

Sizes match day 44, so the two pages use the same three problems. Timing uses CUDA events on the stream where work is queued, three warm-ups and the mean of ten, because the first call to any cuBLAS entry point picks and caches an algorithm. Every path is CUDA_R_32F in and out under CUBLAS_COMPUTE_32F, with CUBLAS_DEFAULT_MATH set in code.

On Ampere and newer GPUs, the compute type and math mode determine whether a run uses FP32 or TF32 tensor cores, so report both settings.

The library's status type is not cudaError_t, so CUDA_CHECK cannot wrap it. A small function checks every library call without adding a second error macro:

static void checkCublas(cublasStatus_t status) {
    if (status != CUBLAS_STATUS_SUCCESS) {
        std::fprintf(stderr, "cuBLAS error: %s\n",
                     cublasGetStatusString(status));
        std::exit(EXIT_FAILURE);
    }
}

Day 44 wrote those branches out by hand and one of them, inside a timing lambda, dropped the status on the floor for a while. A discarded library status is a wrong number with nothing on stderr.

The swap the whole page turns on is two argument positions:

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);
}

And the epilogue is two attributes on the descriptor, one naming the epilogue and one handing over the vector:

    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));
    }

What the fusion should save is arithmetic you can do before the run, the way day 48 did it. At 2048 the output is 16.8 MB; addBias reads it once and writes it once, so 33.6 MB crosses the bus, and day 11 measured 244.7 GB/s of coalesced copy bandwidth. See day 11 for the measurement.

Call it 140 microseconds plus one launch, which day 48 measured at 2.668 microseconds empty, against the 2.858 ms day 44 timed for cublasGemmEx at this size.

Roughly a twentieth, and a larger fraction as the matrices shrink, since the GEMM costs N cubed and the bias N squared.

Results

Verified on a Tesla T4 (compute capability 7.5), driver 580.173.02 and CUDA 12.6. The full build and run transcript is code/day81-cublas/evidence/run-2026-09-02.txt; timings are means of ten runs after three warm-ups.

size path ms GFLOP/s vs fused
512 gemmEx + bias kernel 0.1394 1925.3 1.048
512 Lt + bias kernel 0.1391 1929.9 1.046
512 Lt bias epilogue 0.1330 2018.3 1.000
512 bias kernel alone 0.0099 - 0.075
1024 gemmEx + bias kernel 0.8689 2471.4 1.049
1024 Lt + bias kernel 0.8682 2473.5 1.048
1024 Lt bias epilogue 0.8284 2592.3 1.000
1024 bias kernel alone 0.0410 - 0.049
2048 gemmEx + bias kernel 6.1782 2780.7 1.301
2048 Lt + bias kernel 6.1812 2779.4 1.302
2048 Lt bias epilogue 4.7476 3618.6 1.000
2048 bias kernel alone 0.1357 - 0.029
Five bets, in the order I would be most surprised to lose them.
  1. Held. The two separate-bias API rows differed by 0.22 percent at 512, 0.08 percent at 1024 and 0.05 percent at 2048, all inside 2 percent.
  2. Refuted. Subtracting the 0.1357 ms bias row leaves 6.0425 ms at 2048, nowhere near day 44's 2.858 ms within 5 percent.
  3. Refuted. Fusion saved 1.4306 ms, or 23.2 percent of the separate path, outside the predicted 3 to 7 percent. That saving was over ten times the standalone bias kernel's 0.1357 ms, not within 30 percent.
  4. Partly held, narrowly refuted overall. The bias kernel took 0.1357 ms at 247.4 GB/s, inside both direct bands. That rate is 100.8 percent of day 11's 245.348 GB/s copy row, just outside the predicted 80 to 100 percent ceiling band.
  5. Refuted. Candidate 0 was fastest at 512 and within 0.3 percent at
    1. At 2048 it took 3.8400 ms against candidate 7's 3.5681 ms, 7.6 percent slower and outside the 5 percent band.

The gate under all five held: the heuristic returned eight candidates at every size. All paths matched the CPU reference and the process exited 0.

CUDA 13.0 (V13.0.88) re-verification on the same T4 also compiled, ran all three native cuBLAS/cuBLASLt paths and exited 0. It reported cuBLAS and cuBLASLt version 130002, and every output still matched the CPU reference. The 512 and 1024 timings stayed close to the CUDA 12.6 results.

At 2048, however, the selected fused epilogue moved from 4.7476 to 6.0334 ms, so its advantage over the separate path shrank from 1.301x to 1.024x. Candidate 0 showed the same change, from 3.8400 to 6.0353 ms, while candidate 7 stayed close at 3.5681 versus 3.5873 ms. This is material algorithm-selection drift, not a correctness change.

Both transcripts remain in front matter.

The epilogue saving is the bias kernel's bytes over the GEMM's time, and both parts vary with GPU and size. The milliseconds in this table apply only to the recorded Tesla T4 run at these three sizes.

Run it yourself

Build the program for a GPU with compute capability 7.5 or newer. It does not need a profiler or performance counters. Link both cuBLAS libraries:

nvcc -std=c++17 -O3 -arch=sm_75 -o cublas_lt cublas_lt.cu -lcublas -lcublasLt

cuBLASLt is a separate shared object even though it ships inside the cuBLAS package, so linking only -lcublas without -lcublasLt fails at the link step. Compiler Explorer is out: the 2048 case's CPU reference is 8.6 billion multiply-adds, which walks past the 20 second run cap on its own.

Exercise

Take your best hand-written GEMM, day 44's float4 kernel or the tensor-core version from capstone 4, and replace it with a cuBLAS call at the same dtype, compute type and math mode. Then add the bias to both, once as a separate kernel and once as an Lt epilogue, and report the four numbers.

Time: 30 to 45 minutes. Submit: the four times at N = 2048, your percentage of cublasGemmEx, and one sentence on which bias placement moved that percentage more, and why.

Check: every path is compared against the double CPU reference at every element before anything is timed, and a failure prints the row, column, value, expectation and the bound that was broken, so a mismatched transpose convention fails instead of posting a fast wrong answer. The tolerance is the length-scaled 4 * 2^-23 * sqrt(K), printed with its arithmetic so you can see it was not invented to fit the result.

Hint 1 Both placements move the same bytes and one of them also launches something. Which of your two implementations spends the larger fraction of its time on the bias?
Hint 2 Adding a bias to both sides of your percentage changes them by the same number of milliseconds, not the same proportion. What does adding a fixed constant to the top and bottom of a fraction below one do to it?
Solution Adding the same bias cost to both sides pulls your percentage toward 100, because a fixed addition to the top and bottom of a fraction under one raises it. The separate-bias comparison flatters you, and the slower your kernel the more it flatters you. Fusing the bias into the library's GEMM and leaving yours as a second launch is the honest version, and the one that gets worse for you.

A comparison against a library is only as honest as the work you made both sides do. Every stage you leave outside the timed region is a stage the library is allowed to fuse and you are not.

Pitfalls

Your result is the transpose of what you wanted, or nonsense at the edges. You passed A then B and the library read both as column major. On a square problem that compiles, runs, and returns B*A. Swap the pointers rather than transposing the data, and check a non-square case, since square problems hide the leading dimensions that would have caught you.

Your percentage of cuBLAS looks much better than everyone else's. The first call to any cuBLAS entry point picks and caches an algorithm, so an unwarmed library row is timing a decision as well as a multiply. Warm every path, this program does three runs. The other half of this trap is the compute type: CUBLAS_COMPUTE_32F_FAST_TF32 on an Ampere or newer card is tensor cores, not FP32, and it is one enum away from what you meant.

cublasLtMatmul returns CUBLAS_STATUS_NOT_SUPPORTED and will not say why. The descriptor, the four layouts, the epilogue and the workspace all have to agree, and the call reports the disagreement without naming it. Run cublasLtMatmulAlgoGetHeuristic first and read the count: zero candidates means the combination has no kernel, and moving one attribute at a time until the count changes is faster than reading the support matrix.

The call fails with CUBLAS_STATUS_ALLOC_FAILED, or your user workspace vanished. Both are the workspace. On size, the documentation says "workspace size equal to or larger than 16KiB is enough to prevent CUBLAS_STATUS_ALLOC_FAILED error, while a larger workspace can provide performance benefits for some routines" (https://docs.nvidia.com/cuda/cublas/index.html#cublassetworkspace , checked 2026-09-01). On disappearance, cublasSetStream resets the cuBLAS workspace to the default pool unconditionally and says nothing, so the order is stream first, workspace second.

cuBLASLt sidesteps both by taking the workspace as an argument on every call.

Go deeper

Next

Day 82 does this to four more libraries at once, and one of them gets a real fight: cuSPARSE's SpMV against day 37's two hand-written kernels, the one comparison in this module where the input picks the winner. Day 84 opens CUTLASS, where the kernels the heuristic was choosing between come from, and it needs compute capability 8.0. The occupancy and launch overhead this page leaned on are days 45 and 48, and the fusion argument returns on day 96, where the epilogue you cannot get from a library is the one you write yourself.