Setup
GPU-BLAST 1.1 is a 2013 CUDA patch on top of NCBI BLAST+ 2.2.28 that offloads the two-hit protein sequence alignment step to GPU. I ran it against Swiss-Prot (575,503 proteins, 208 M residues) with 100 query proteins on a Tesla T4 (sm_75, 15 GiB GDDR6, 70 W TDP) running CUDA 11.8.
The build was the first hurdle — seven iterations to navigate CUDA 12 removing
the legacy texture<> API, five C++ compatibility patches against g++11, and a
manual linker fix so blastp actually linked against libgpublast.a. A 2013
codebase on a 2026 Ubuntu toolchain requires that kind of patience.
Wall time
| Mode | Time | Speedup |
|---|---|---|
blastp 2.2.28 — 1 thread |
219.8 s | 1× |
blastp 2.2.28 — 4 threads |
67.7 s | 3.2× |
GPU-BLAST T4 (-gpu T, 1 thread) |
103.4 s | 2.13× |
Hit count is 577 for all three runs — the GPU path is bit-exact with CPU. The
problem: four CPU threads beat it. And the GPU path is single-threaded only;
-num_threads 4 -gpu T crashes because the device state is shared across search
threads without locking.
What nvprof shows
GPU activities (5 queries)
89.9% 712.5 ms GPU_BLASTP_kernelTwoHit ← one launch per batch
6.0% 48.0 ms [CUDA memcpy HtoD] ×10
3.9% 30.6 ms [CUDA memcpy DtoH] ×2
Kernel config: 512 blocks × 64 threads = 32,768 threads (2 warps/block)
API overhead (first run):
307.9 ms cudaMalloc ×9 ← 198 MB GPU database
127.3 ms cudaThreadExit
One kernel. 90% of GPU time. 64 threads per block, which is 2 warps.
Roofline: the 1% problem
ncu on sm_75, 1 query:
| Counter | Value |
|---|---|
| Global load bytes | 447.93 MB |
| Integer instructions | 5.374 B |
| Kernel time | ~142.5 ms |
| SM throughput | 49.7% of peak |
Arithmetic intensity = 5.374B / 447.93 MB = 12 ops/byte
That's left of the T4 ridge point (27 ops/byte), so the kernel is nominally memory-bound. At that intensity, the bandwidth ceiling is 12 × 300 GB/s = 3,600 GOPS. The kernel achieves 37.7 GOPS — 1.05% of the ceiling.
This is not a bandwidth problem. 447 MB at peak T4 bandwidth transfers in ~1.5 ms. The kernel runs for 142.5 ms — a 95× gap. The bottleneck is access latency, not throughput.
Why: scoring matrix lookup
The inner loop of BLAST protein alignment is:
int score = score_matrix[query_aa * 20 + subject_aa];
query_aa and subject_aa are residue types from two different sequences.
They're data-dependent: adjacent threads processing adjacent database sequences
read completely different positions in the substitution table. No spatial
locality. No temporal locality within a warp's working set.
Every access to the 400-byte BLOSUM62 matrix in global memory is potentially a cache miss. On T4, a DRAM miss costs ~400 cycles. With 64 threads per block (2 warps), there aren't enough warps in flight to context-switch away from that stall. The 49.7% SM throughput — half the time stalled — makes this concrete.
The roofline sits at 3,600 GOPS. The kernel is at 37.7. The gap is pure latency overhead accumulated one scoring lookup at a time.
Why 2.13× instead of the paper's 4×
The original GPU-BLAST paper measured ~4× on a GTX 580 (Fermi, 244 W TDP, 48 warps/SM at full occupancy). The T4 is a 70 W inference chip: 32 warps/SM, no overclocking headroom, and the lower resident warp count means less latency hiding for irregular workloads. Higher peak bandwidth (300 vs 192 GB/s) doesn't help when bandwidth isn't the bottleneck.
What would fix it
Score matrix in shared memory is the single highest-value change — 4 lines of code, 400 bytes, and lookup latency drops from ~400 cycles to 4:
__shared__ int8_t s_blosum[20][20];
if (threadIdx.x < 400)
s_blosum[threadIdx.x / 20][threadIdx.x % 20] = blosum_global[threadIdx.x];
__syncthreads();
// now: score = s_blosum[query_aa][subject_aa];
Estimated gain on the inner loop alone: 10–50×. Combined with larger thread
blocks (64 → 512 threads, so 16 warps/block) and fixed thread-safety for
-num_threads 4, the GPU path would easily clear the 4-CPU-thread time of 67.7
s. Strategies, in order:
| Priority | Strategy | Lines of code | Expected gain |
|---|---|---|---|
| 1 | Score matrix → shared memory | ~10 | 5–20× inner loop |
| 2 | Thread block size 64 → 256+ | 1 | 2–4× occupancy |
| 3 | Thread-safe streams (-num_threads 4) |
Medium | 4× throughput |
| 4 | Pipelined H2D/kernel overlap | Medium | 1.5–2× |
| 5 | Coalesced database layout | High (DB rebuild) | 2–3× |
| 6 | Two-phase seeding/extension (Diamond's design) | Architecture | 3–10× |
None of this was obvious in 2013 without ncu. The paper's authors knew the kernel was memory-bound — the per-cycle latency diagnosis required hardware counters that weren't accessible until later tooling.
The build lesson
Seven iterations to compile a 2013 CUDA codebase on Ubuntu 22.04:
| Error | Fix |
|---|---|
texture<> removed in CUDA 12 |
CUDA 11.8 alongside 12.6 |
IStringDecoder* GetEncoder() type error |
Patch ncbistr.hpp return type |
Two comparator operator() missing const |
scheduler.cpp, thread_pool.cpp |
dynamic exception specs in C++17 mode |
-std=c++14 in both Makefile CXXFLAGS |
swap static_assert (hash_impl) |
Whole build via g++-9 |
Linker: undefined GPU_BLASTP_* symbols |
Add -lgpublast -lcudart to blastp Makefile |
.gpuinfo truncated on write |
Sort Swiss-Prot by length before method=2 |
The T4 added one more: nvprof metrics are blocked on sm_75+, so all hardware counter work used ncu instead.
The argument for containerizing GPU software at the version it was built on, not the version the cluster has, is right there in that table.
Takeaway
A kernel at 1% roofline efficiency in the nominally-memory-bound regime is not limited by bandwidth. It's spending most of its time stalled on DRAM latency, one irregular access at a time, with too few warps to hide the stall. The roofline tells you the ceiling; ncu's per-counter breakdown tells you why you're not reaching it.
BLAST's scoring matrix is a textbook worst case for this failure mode — small table, fully random access pattern, no locality to exploit. Four lines of shared memory code would change the answer materially. That's the real finding from the 1% number.