Split View: cuDNN 완전 해부: 딥러닝 연산이 GPU에서 빛처럼 빠른 이유

Naive CUDA convolution:     ~5 TFLOPS achieved (0.5% of H100 theoretical)
cuDNN optimized convolution: ~900 TFLOPS achieved (90% of H100 theoretical)
Difference: ~180x

import torch

# This code calls cuDNN internally
conv = torch.nn.Conv2d(64, 128, kernel_size=3, padding=1).cuda()
x = torch.randn(32, 64, 56, 56, device='cuda')
y = conv(x)  # → calls cudnnConvolutionForward()

# Verify cuDNN is in use
print(torch.backends.cudnn.enabled)    # True
print(torch.backends.cudnn.version())  # e.g., 8906

For each output pixel (i, j):
  for c_out in range(C_out):
    for c_in in range(C_in):
      for kh in range(K):
        for kw in range(K):
          output[c_out, i, j] += input[c_in, i+kh, j+kw] * filter[c_out, c_in, kh, kw]

Total ops: N × C_out × C_in × K^2 × H_out × W_out

Input feature map [3x3 image, 3x3 kernel]:

Original input:        im2col result:
┌─────────┐            ┌────────────────────────────┐
│ 1  2  3 │  im2col    │ 1  2  4  5  (patch 0, pos 0,0) │
│ 4  5  6 │ ─────────→ │ 2  3  5  6  (patch 1, pos 0,1) │
│ 7  8  9 │            │ 4  5  7  8  (patch 2, pos 1,0) │
└─────────┘            │ 5  6  8  9  (patch 3, pos 1,1) │
                       └────────────────────────────┘
Each row = all pixel values in one receptive field

Then: output = filter_matrix × im2col_matrix
      (standard GEMM → Tensor Core fully utilized!)

# Python illustration of im2col behavior
import numpy as np

def im2col(input, kernel_h, kernel_w, stride=1, pad=0):
    N, C, H, W = input.shape
    out_h = (H + 2*pad - kernel_h) // stride + 1
    out_w = (W + 2*pad - kernel_w) // stride + 1

    # Apply padding
    img = np.pad(input, [(0,0),(0,0),(pad,pad),(pad,pad)], mode='constant')

    # Build im2col matrix: (C*kernel_h*kernel_w) x (N*out_h*out_w)
    col = np.zeros((N, C, kernel_h, kernel_w, out_h, out_w))
    for j in range(kernel_h):
        jj = j + stride * np.arange(out_h)
        for i in range(kernel_w):
            ii = i + stride * np.arange(out_w)
            col[:, :, j, i, :, :] = img[:, :, jj[:, None], ii[None, :]]

    col = col.transpose(0, 4, 5, 1, 2, 3).reshape(N*out_h*out_w, -1)
    return col
# cuDNN's actual im2col is far more optimized than this

Standard 3x3 convolution (2x2 output):
- Input patch: 4x4 = 16 elements
- Kernel: 3x3 = 9 elements
- Output: 2x2 = 4 elements
- Multiplications required (naive): 4 x 9 = 36

After Winograd F(2x2, 3x3) transformation:
- Transformed input: 4x4 = 16 elements (linear transform, no multiplications)
- Element-wise multiplication: 16 operations (4x4)
- Inverse transform: 4x4 → 2x2 (linear transform, no multiplications)
- Multiplications required: 16
- Reduction: 36 → 16 = 2.25x fewer multiplications

Y = A^T [(G × g × G^T) element-wise-multiply (B^T × d × B)] A

Where:
d = input tile (4x4)
g = 3x3 kernel
B, G, A = fixed transform matrices (precomputed constants)
element-wise-multiply = Hadamard product

Key insight: G × g × G^T is computed once per kernel
             (can be precomputed at inference time)
             B^T × d × B is computed per input tile
             Element-wise product gives 16 multiplications
             instead of 36

Spatial domain convolution:    O(N × K^2)   (N = output size, K = kernel size)
Frequency domain convolution:  O(N × log N) (FFT + element-wise multiply)

The larger K is, the more FFT wins:
K=3:  9 ops vs log(N) ≈ 7  → marginal difference
K=11: 121 ops vs log(N) ≈ 7 → FFT uses 17x fewer multiplications

# What does this single line actually do?
torch.backends.cudnn.benchmark = True

# See the benchmark effect directly
import torch
import time

torch.backends.cudnn.benchmark = False
model = resnet50().cuda()
x = torch.randn(32, 3, 224, 224).cuda()

t0 = time.time()
for _ in range(10): y = model(x)
torch.cuda.synchronize()
print(f"benchmark=False: {(time.time()-t0)*1000:.0f}ms")

torch.backends.cudnn.benchmark = True
# First run includes benchmarking (SLOW!)
t0 = time.time()
y = model(x)  # This single call takes several seconds
torch.cuda.synchronize()
print(f"First run (includes benchmarking): {(time.time()-t0)*1000:.0f}ms")

# Subsequent runs use cached optimal algorithm
t0 = time.time()
for _ in range(10): y = model(x)
torch.cuda.synchronize()
print(f"benchmark=True (after warmup): {(time.time()-t0)*1000:.0f}ms")
# Typically 20-40% faster

y = gamma * (x - mean(x)) / sqrt(var(x) + eps) + beta

Naive implementation memory round trips:
1. Load x from HBM
2. Compute mean → store to HBM
3. Load x + load mean (2 HBM reads)
4. Compute variance → store to HBM
5. Load x + load mean + load var (3 HBM reads)
6. Compute normalization → store to HBM
7. ReLU: load from HBM → compute → store to HBM

Total: ~10 HBM accesses

Fused BN+ReLU kernel:
1. Load x into registers/shared memory (1 HBM read)
2. Compute mean, variance in registers (zero HBM access)
3. Normalize in registers (zero HBM access)
4. ReLU in registers (zero HBM access)
5. Store result to HBM (1 HBM write)

Total: 2 HBM accesses (5x reduction!)

# Using fused BN+ReLU in PyTorch
import torch.nn as nn

# Naive: Conv → BN → ReLU (3 separate kernel launches)
class NaiveBlock(nn.Module):
    def __init__(self, c):
        super().__init__()
        self.conv = nn.Conv2d(c, c, 3, padding=1)
        self.bn = nn.BatchNorm2d(c)
        self.relu = nn.ReLU()

    def forward(self, x):
        x = self.conv(x)
        x = self.bn(x)    # Separate cuDNN kernel
        x = self.relu(x)  # Separate elementwise kernel
        return x

# With torch.compile(), the compiler automatically fuses these
# through cuDNN Graph API and nvFuser
model = torch.compile(model)  # Auto-fuses conv+bn+relu into one kernel

# Standard Attention — memory inefficient
def standard_attention(Q, K, V, scale):
    # S = Q × K^T: creates N×N attention score matrix
    S = torch.matmul(Q, K.transpose(-2, -1)) * scale  # O(N^2) memory!
    # Softmax
    P = torch.softmax(S, dim=-1)
    # Multiply by V
    O = torch.matmul(P, V)
    return O

# Problem: S requires O(N^2) memory
# N=8192 (document length), FP16:
# S size = 8192 * 8192 * 2 bytes = 128MB per head!

Standard attention memory flow:
Q, K, V → read from HBM
S = QK^T → write to HBM (N^2 size!)
softmax(S) → write to HBM
P × V → write to HBM
Total HBM accesses: O(N^2)

FlashAttention memory flow:
Load Q, K, V in tiles → into shared memory
Compute partial attention sums in shared memory
  (never creates full N^2 matrix)
Write only the final output to HBM
Total HBM accesses: O(N)  ← from N^2 down to N!

FlashAttention tiling:

Partition Q into blocks Q_1, Q_2, ... Q_Tc
Partition K, V into blocks K_1, V_1, ... K_Tr, V_Tr

for i in range(Tc):
  Load Q_i into SRAM
  O_i = 0, l_i = 0, m_i = -inf  (softmax statistics init)

  for j in range(Tr):
    Load K_j, V_j into SRAM

    # Compute entirely in SRAM (zero HBM access)
    S_ij = Q_i × K_j^T  (tile-sized attention scores)
    m_ij = max(m_i, rowmax(S_ij))  (numerical stability)
    P_ij = exp(S_ij - m_ij)

    # Online softmax update (numerically stable!)
    O_i = diag(exp(m_i - m_ij)) × O_i + P_ij × V_j
    l_i = exp(m_i - m_ij) × l_i + rowsum(P_ij)
    m_i = m_ij

  # Final normalization, then write to HBM (1 write)
  O_i = diag(l_i)^(-1) × O_i → store to HBM

NCHW layout (batch × channel × height × width):
Batch=1, channels=3 (RGB), 4x4 image:

Memory arrangement:
[R00 R01 R02 R03 | R10 R11 ... | R30 R31 R32 R33 |
 G00 G01 G02 G03 | G10 G11 ... | G30 G31 G32 G33 |
 B00 B01 B02 B03 | B10 B11 ... | B30 B31 B32 B33]

→ Pixels within the same channel are contiguous

NHWC layout (batch × height × width × channel):
Memory arrangement:
[R00 G00 B00 | R01 G01 B01 | R02 G02 B02 | R03 G03 B03 |
 R10 G10 B10 | R11 G11 B11 | ...
 R30 G30 B30 | R31 G31 B31 | R32 G32 B32 | R33 G33 B33]

→ All channels at the same spatial position are contiguous

# Using NHWC (channels_last) in PyTorch
x_nchw = torch.randn(32, 64, 56, 56, device='cuda')

# Convert to NHWC (channels_last format)
x_nhwc = x_nchw.to(memory_format=torch.channels_last)

# Convert model to channels_last format
model = model.to(memory_format=torch.channels_last)

# Forward pass will auto-select NHWC kernels in cuDNN
output = model(x_nhwc)

# Benchmark comparison
# Measured on H100:
# NCHW: ~12ms/batch
# NHWC: ~9ms/batch  (~25% faster)

Original model (ONNX/PyTorch)
      ↓
  Graph analysis
      ↓
  Layer fusion
  ┌─────────────────────────────────────────────┐
  │ Conv → BN → ReLU → Conv → BN → ReLU        │
  │           ↓ after fusion                    │
  │     single optimized kernel                 │
  └─────────────────────────────────────────────┘
      ↓
  Precision selection (FP32 → FP16 → INT8)
      ↓
  Kernel auto-selection (benchmarked per input size)
      ↓
  Optimized execution engine

import tensorrt as trt

# INT8 calibrator setup
class MyCalibrator(trt.IInt8EntropyCalibrator2):
    def __init__(self, data_loader, cache_file):
        super().__init__()
        self.data_loader = iter(data_loader)
        self.cache_file = cache_file
        self.batch_allocation = None

    def get_batch_size(self):
        return 32

    def get_batch(self, names):
        try:
            batch = next(self.data_loader)[0].numpy()
            # Calibrate activation ranges using representative dataset
            # TensorRT maps FP32 range of each layer to INT8
            if self.batch_allocation is None:
                self.batch_allocation = cuda.mem_alloc(batch.nbytes)
            cuda.memcpy_htod(self.batch_allocation, batch)
            return [int(self.batch_allocation)]
        except StopIteration:
            return None

# Builder configuration
builder = trt.Builder(logger)
config = builder.create_builder_config()
config.set_flag(trt.BuilderFlag.INT8)  # Enable INT8
config.int8_calibrator = MyCalibrator(calibration_loader, 'cache.bin')

# Trace LLM operations using PyTorch Profiler
import torch
from torch.profiler import profile, ProfilerActivity

model = GPT2Model.from_pretrained('gpt2').cuda()
input_ids = torch.randint(0, 50257, (1, 512)).cuda()

with profile(
    activities=[ProfilerActivity.CPU, ProfilerActivity.CUDA],
    with_stack=True
) as prof:
    output = model(input_ids)

# Print top CUDA kernels
print(prof.key_averages().table(sort_by="cuda_time_total", row_limit=10))

GPT-2 Forward Pass Internal CUDA Kernels (512 tokens, batch=1):

Operation       Kernel                         Time   Share
──────────────────────────────────────────────────────────
Linear (QKV)    cublasSgemm (Tensor Core)      2.1ms  35%
Attention       flash_attn_fwd_kernel          1.4ms  23%
Linear (output) cublasSgemm (Tensor Core)      0.9ms  15%
LayerNorm       layer_norm_kernel              0.3ms   5%
GELU            vectorized_elementwise_kernel  0.2ms   3%
Residual Add    vectorized_elementwise_kernel  0.1ms   2%
Embedding       Embedding_cuda                 0.1ms   2%
Other                                         0.9ms  15%

Total: ~6ms (on H100)