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

cuFFT, cuRAND, cuSPARSE and cuSOLVER

Take a thousand samples, run them through cuFFT forwards, run the result back through cuFFT inverse, and compare with what you started with. Every value comes back a thousand times too big. Nothing is broken and no flag was missed.

The library documents it in one sentence: "cuFFT performs un-normalized FFTs; that is, performing a forward FFT on an input data set followed by an inverse FFT on the resulting set yields data that is equal to the input, scaled by the number of elements" (https://docs.nvidia.com/cuda/cufft/index.html section 3.15.3, checked 2026-09-01).

Four libraries in the toolkit replace common custom kernels. Each adds a cost beyond the launch: a missing divide, a state array, a scratch buffer, or a plan. This lesson shows which library to call and how to include that extra cost in a comparison.

Four libraries, four things you stop writing

The math libraries ship with the CUDA toolkit. Each needs its own -l flag. On day 39, CUB beat the hand-written reduction by 1.38x and the scan by 2.33x; the custom kernels still help explain the library results.

Sorted by what each one replaces:

cuRAND replaces the random number generator you were going to write in four lines. It has two halves. The host API fills a device buffer with a single call, the way you would use a Python RNG.

The device API gives you curand_init and curand_uniform inside your own kernel, so the numbers never touch memory.

cuSPARSE replaces the SpMV kernel of day 37. That day measured the thing this one has to answer: on a skewed matrix, one thread per row issues 219,676,672 lane slots and 3.8 percent of them do work, while one warp per row wins there and loses on an even matrix. You cannot pick between those two from the source, so the question for the library is whether it picks for you.

cuFFT provides fast Fourier transforms. Build a plan once and execute it many times; the plan stores the batch count.

cuSOLVER provides blocked, pivoted LU factorization. Use it instead of implementing the full factorization from CUDA kernels.

Diagram: how much of the work you still own, in four steps. One horizontal axis, left to right, from "your kernel" to "the library's algorithm". Four stacked bands, each a single row split into a host side and a device side by a vertical launch line, with the boxes you write shaded and the boxes the library owns hollow. Band 1, your own kernel: one shaded host launch and one shaded device box holding the generator, the loop and the reduction. Caption "1 launch, 0 library calls, and every bug is yours." Band 2, cuRAND device API: two shaded host launches, and inside the second device box a hollow curand_uniform sitting in your loop. Caption "2 launches, 0 host library calls, the loop is still yours." Band 3, cuSPARSE and cuFFT: four hollow host boxes (create, size the buffer or the plan, execute, destroy) and a hollow device box you never see. Caption "4 host calls, 0 kernels of yours, and a buffer you allocate." Band 4, cuSOLVER: two hollow host boxes, getrf then getrs, over one hollow device box. Caption "2 host calls for the algorithm you were never going to write." Alt text: "Four steps from owning a kernel to owning nothing. Your own kernel is one launch and every bug. cuRAND's device API keeps your loop. cuSPARSE and cuFFT take four host calls and a buffer. cuSOLVER takes two."

The call is never just the call

A library call looks like it costs one launch, because one line replaced thirty. Every one of the four charges something else on top, and in three of them the extra term is a whole pass over memory.

curand_init is a kernel of its own, and it writes one generator state per thread. The documentation says so plainly: "Calls to curand_init() are slower than calls to curand() or curand_uniform()" (https://docs.nvidia.com/cuda/curand/device-api-overview.html , checked 2026-09-01). cusparseSpMV asks for a scratch buffer whose size you have to query first, exactly as cub::DeviceScan did on day 39.

cufftPlan1d allocates a work area that can be as large as the data. And the host half of cuRAND writes its numbers to global memory so a second kernel can read them back, which is 128 MiB of traffic the device API never moves.

The program below times setup separately from the work. A benchmark that hides curand_init inside the sampling loop is measuring a different program from the one you will ship.

Holding the data still and changing who computes

Full program in code/day82-math-libs/math_libs.cu. Four parts, one binary, no arguments.

The comparison is against your own code, not against another library. Part 1's cuRAND rows run against a hand-written LCG on the same estimate of pi from the same number of samples. Part 2 keeps day 37's two matrices and both of its kernels unchanged, so cusparseSpMV is measured on the data day 37 already published numbers for.

Setup is timed on its own line. curand_init, the host generator's fill, the cuSPARSE buffer query, the two plans' work areas: each is reported beside the work rather than folded into it.

Correctness is checked before anything is timed, and never against the library itself. The transforms are checked against the analytic spectrum of a cosine, the solve against a planted answer, the four SpMV paths against a double-precision CPU reference, and the two cuRAND estimates against pi to within five standard errors of the estimator. Comparing cuFFT with cuFFT would pass on a plan that transforms the wrong batch.

The tolerances scale with the reduction depth and print their arithmetic. The FFT gate is rtol_K = max(1e-5, 4*2^-23*sqrt(K)) at K = 1024, from day 66's table; a radix-2 FFT's error grows like log2(K) rather than sqrt(K), so the bound has room and the program says so on its own output. The dense solve uses the f64 row of the same table.

Seeding is the part of cuRAND people get wrong, so it gets its own kernel:

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

One seed for the whole run, the thread index as the sequence number, offset zero. That order is the documented one: "each experiment should be assigned a unique seed. Within an experiment, each thread of computation should be assigned a unique sequence number" (same page, checked 2026-09-01).

Handing each thread a different seed instead is the common version, and it is guesswork: two seeds have no defined relationship, so nothing keeps the two streams apart.

The cuSPARSE half is four calls, and the second one exists only to tell you how much scratch to allocate:

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

CUSPARSE_CHECK is one of four macros this file carries, one per library, each in the shape of the course's CUDA_CHECK. The four libraries return four unrelated status enums, so CUDA_CHECK cannot wrap any of them, and day 44 shipped a bug by letting a cuBLAS status fall on the floor inside a timing lambda. Only cuSPARSE has a function that turns its status into a string; the other three print the integer, which is worth knowing before you are reading one at midnight.

And the cuFFT plans, where the whole batching story is one argument:

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

Note. This matrix cannot show what CUSPARSE_SPMV_CSR_ALG2 is for. Its values are integers small enough that every partial sum is exact in float, so all four SpMV paths must agree to the last bit whatever order they add in, and the two algorithms differ here only in time. On real values ALG2 is the one documented to give "deterministic (bit-wise) results for each run" (https://docs.nvidia.com/cuda/cusparse/index.html section 6.6.7, checked 2026-09-01), and run-to-run reproducibility is day 68's subject and a real reason to pick it.

Results

Verified on a Tesla T4, driver 580.173.02, CUDA 12.6. The passing rerun ended with all checks passed and exit code 0. The first attempt is retained separately: its FFT oracle expected N/2 at DC, but the positive and negative cosine lines coincide there and add to N. That made the oracle wrong by 512; after correcting it, the worst spectrum error was 9.15527e-05 against a 7.812e-03 gate.

generator setup ms draw ms total ms Mdraw/s pi sigma
hand LCG - 0.582 0.582 57669.8 3.142054 1.15
cuRAND device API 0.042 0.473 0.515 70926.9 3.142070 1.19
cuRAND host API 0.759 0.533 1.292 44212.2 3.141763 0.42
matrix path ms GFLOP/s buffer bytes
even thread per row 0.609 27.6 -
even warp per row 0.983 17.1 -
even cuSPARSE ALG1 0.390 43.1 32772
even cuSPARSE ALG2 0.379 44.2 164484
skewed thread per row 1.772 9.5 -
skewed warp per row 1.196 14.0 -
skewed cuSPARSE ALG1 0.366 45.8 32772
skewed cuSPARSE ALG2 0.380 44.2 164484
result measurement
batched FFT 0.074 ms
1024 batch-1 FFT calls 5.540 ms
FFT loop / batch 74.77x
batched / single-plan work area 8192 / 8 bytes
worst spectrum error / gate 9.15527e-05 / 7.812e-03
unnormalized round-trip scale 1023.9999
worst round-trip error 7.7486e-07
cusolverDnDgetrf 14.180 ms; error 1.177e-14, gate 2.010e-12
cusolverDnDgetrs 3.851 ms; error 1.182e-11, gate 2.747e-09
Five predictions, now settled.

  1. Mixed. cuSPARSE wins by more on the skewed matrix than on the even one: 3.27x against the best hand kernel versus 1.61x. The claim that it could not win by much on the even matrix is not verified; 1.61x is a material win. The same-process rows are the comparison; day 37's older absolute timings are not substituted into this ratio.

  2. Verified. The host half of cuRAND did not beat the device half. Its

    1.292 ms total was 2.51x the device API's 0.515 ms. curand_init was 0.042 ms, 8.9 percent of the device sampling pass. It writes 128 MiB of draws and a second kernel reads them back. The measured

    244.7 GB/s copy ceiling gives those two rows a floor of about half a millisecond each before any arithmetic. The device API's loop touches no global memory at all.

  3. Verified. The batched plan beat the loop by 74.77x. A batch-1 plan called 1024 times pays at least 1024 empty launches, and day 48 measured one at 2.668 microseconds on the same Tesla T4, so that path has a floor near three milliseconds. The batched call moves 16 MiB.

  4. Verified. The measured unnormalized round-trip scale was 1023.9999.

  5. Not verified. getrs was only 3.68x cheaper than getrf, 3.851 ms against 14.180 ms, not between 50x and 400x. The factorization is two thirds of n cubed and the solve is two n squared, but that operation-count ratio did not predict the measured library-call ratio.

Compare the cuSPARSE margin between the two matrices and the ratio between setup and work. Keep the Tesla T4 label with the milliseconds in this table.

The native CUDA 13.0 (V13.0.88) rerun also passed every gate and exited 0. It printed cuRAND 10.4.0, cuFFT 12.0.0, cuSPARSE 12.6.3 and cuSOLVER 12.0.4, versus 10.3.7, 11.3.0, 12.5.4 and 11.7.1 in the CUDA 12.6 run. The cuRAND, cuSPARSE and cuFFT results retained the same ordering and similar timings.

cuFFT's reported plan work areas changed from 8192 and 8 bytes to zero for both plans. The material timing drift was cusolverDnDgetrs, which moved from 3.851 to 55.482 ms while getrf moved from 14.180 to 15.644 ms. The original prediction verdict remains tied to the CUDA 12.6 table above; CUDA 13 reverses that relationship, with the solve 3.55x slower than the factorization in this single run.

The CUDA 13 transcript contains no ldd capture, so it establishes the printed component versions and behavior, not exact loaded SONAMEs or paths.

Run it yourself

Build for compute capability 7.5 or newer. Put all four libraries on the link line:

nvcc -std=c++17 -O3 -arch=sm_75 -o math_libs math_libs.cu \
    -lcurand -lcufft -lcusparse -lcusolver

The repo's README carries the exact conda-forge packages and versions, with the CUDA major pinned, because a binary linked against a CUDA 13 library while compiling with a 12.6 toolkit is how day 44 lost a session. There is no Compiler Explorer embed: the two CPU SpMV references walk 8,388,608 nonzeros each before a single library call runs, which is gone past the 20 second execution cap on its own.

Exercise

Add cusparseSpMV_preprocess to the timed path on both matrices and report what it changes. The call is optional, it takes the same arguments as the SpMV it precedes, and the documentation says it "may accelerate subsequent calls to cusparseSpMV()" when the sparsity pattern is reused. Report the four ratios (ALG1 and ALG2, even and skewed) against the numbers in the table above, and say in one sentence which matrix it helped and why.

Time: 30 to 45 minutes. Submit: the four ratios and the sentence.

Check: the harness reruns the same double-precision reference over all 524,288 rows of both matrices after your change, so a preprocess call whose buffer you then reused for something else fails instead of getting faster. It prints the first bad row with what it got and what it wanted, and it prints all eight timings pass or fail, because the timings are the answer and the pass is only the gate.

Hint 1 The preprocess call is not free and it is not per launch. Work out how many SpMV calls it has to amortise over before it can pay for itself, then look at which of the two matrices has an imbalance worth planning around.
Hint 2 The documentation attaches conditions to the buffer, not to the matrix: the buffer becomes "active" for that matrix, the contents must be unmodified between the two calls, and only alpha, beta, the vectors and the matrix values may change afterwards. What does that rule out in a program that shares one scratch allocation between two algorithms?
Solution The buffer is the answer. `cusparseSpMV_preprocess` writes acceleration data into the scratch you pass it and marks that buffer active for that matrix, so a program that allocates one buffer and hands it to both algorithms in turn has invalidated the first one's work before the timing starts. Give each algorithm its own allocation, or preprocess again after switching.

The ratio, whatever it comes out at, is a statement about reuse rather than about cuSPARSE. Preprocessing is a bet that the sparsity pattern outlives the values: true of a solver iterating on one matrix, false of a program that builds a new one each step.

Pitfalls

Your inverse FFT gives answers N times too large. cuFFT does not normalize, and neither direction is scaled: "Scaling either transform by the reciprocal of the size of the data set is left for the user to perform as seen fit" (https://docs.nvidia.com/cuda/cufft/index.html section 3.15.3, checked 2026-09-01). Divide once, on whichever side of the pair you choose, and write down which.

Your cuSPARSE call returns success and writes nothing useful. Like cub::DeviceScan on day 39, cusparseSpMV needs a scratch buffer whose size you query first with cusparseSpMV_bufferSize. Unlike CUB, it takes that buffer as an ordinary argument, so passing a null pointer to a non-zero request is a crash or a wrong answer rather than a "sizing" call.

You give every thread its own seed. cuRAND's parallelism lives in the sequence argument, not the seed argument. Same seed, different sequence number per thread, is the documented pattern, and it is the one the program here uses.

Your uniform draws include 1.0 and exclude 0.0. Both halves of cuRAND say so: the host generator returns "values between 0.0f and 1.0f, excluding 0.0f and including 1.0f" (https://docs.nvidia.com/cuda/curand/group__HOST.html , checked 2026-09-01), and curand_uniform "may return from 0.0 to 1.0, where 1.0 is included and 0.0 is excluded" (https://docs.nvidia.com/cuda/curand/device-api-overview.html , checked 2026-09-01). That is the opposite convention from most host RNGs, and it turns 1.0f / (1.0f - u) from a rare surprise into a division by zero.

You pass row-major data to cuSOLVER. Every LAPACK-shaped API on the GPU is column-major, so A(i, j) lives at A[i + j * n]. A symmetric test matrix hides this and a rectangular one does not.

You ignore devInfo. cusolverDnDgetrf returns a status for the call and writes a separate integer on the device for the factorization: negative means an argument is wrong, positive i means U(i,i) is exactly zero, so the matrix is singular and everything downstream is noise. The status can be CUSOLVER_STATUS_SUCCESS while devInfo says the answer is worthless.

Go deeper

Next

Day 83 is cuDNN and its graph API, where the choice stops being one call and becomes a description of an operation graph that the library picks an engine for. The two things this page measured come back there in a harder form: the setup cost is now a plan and a heuristic search, and the correctness check is against PyTorch rather than against arithmetic you can do on paper. Day 90 is the checkpoint that reads all four capstones back and asks which of their kernels should have been one of these calls.