Going Further with CUDA for Python Programmers

L003 의 마지막 예제는 단순 CUDA matmul 이었다. 결과는 정확했지만 — Jeremy 자신이 강의 끝에 인정한 사실 — “cuBLAS 와 비교하지 않았다, CUDA expert 가 아니라서.” L005 가 그 자리에서 시작한다. 같은 matmul 을 한 단계 더 빠르게 만드는 한 가지 도구 — shared memory — 를 끝까지 깐다. 그 도구를 통해 PMPP Ch.5 의 메모리 모델 (specifically tiling) 이 실전에서 어떻게 사용되는지를 보여준다.

강의가 답하는 질문 두 개.

1. L003 의 naive matmul 이 왜 느린가 — 산술적으로 정확히. (§02 의 data reuse 회계.)
2. shared memory 한 도구로 어디까지 빠르게 만들 수 있나 — Python 시점에서 먼저 흉내내고, CUDA 로 mechanical 변환.

matmul C = A·B 에서 A 가 (M×K), B 가 (K×N) 이면 C 는 (M×N). naive 커널은 한 thread = 한 출력 element 패턴. thread (i, j) 가 inner loop k 를 돌면서 A[i, k] 와 B[k, j] 를 K 번 읽는다. 이 사실을 출력 전체로 펼쳐 보면.

• A 의 row i 는 C 의 row i 의 모든 출력 element (N 개) 가 다 읽는다 → A 의 한 row 가 N 번 다시 읽힌다.
• B 의 column j 는 C 의 column j 의 모든 출력 element (M 개) 가 다 읽는다 → B 의 한 column 이 M 번 다시 읽힌다.
• 총 HBM read = M·K·N + K·N·M = 2·M·N·K, FLOPs = 2·M·N·K (한 element 당 K 번의 multiply-add). arithmetic intensity = 1 FLOP / 4 byte (FP32). roofline 위에서 far-left, memory-bound.

L003 의 naive 커널이 cuBLAS 보다 “많이” 느린 첫 이유가 여기에 있다 — HBM 의 같은 자리를 N번 다시 가져온다. data reuse 가 있는데 그걸 잡지 않은 것. 이걸 잡으려면 한 번 가져와서 on-chip 에 저장하고 N번 재사용이 답이다. 그게 shared memory 의 역할.

강의에서 Jeremy 가 깐 정의 — “shared memory 는 한 SM 위에서 도는 한 block 안의 thread 들에게만 보이는 on-chip SRAM 이다.” 이 한 줄이 모든 디테일을 따라온다.

이 정의에서 직접 따라 나오는 사용 패턴 — 한 block 의 thread 들이 협력해서 HBM 의 일부를 shared 로 가져오고, 그 위에서 inner-loop 를 돌리고, 끝나면 결과만 HBM 으로 다시 쓴다. 이게 “tiling” 의 본질.

tiled matmul 의 분해. 출력 C 를 TILE×TILE 짜리 tile 들로 자르고, 한 block 이 한 출력 tile 을 책임진다. 한 출력 tile 을 계산하려면 A 의 한 row-band 와 B 의 한 column-band 가 필요한데, 그 band 들을 한꺼번에 shared 로 가져올 순 없다 (너무 큼). 그래서 TILE 크기로 잘라 phase 를 돈다 — phase 마다 A 의 한 tile, B 의 한 tile 을 shared 로 load 하고, 그 위에서 partial sum 을 누적.

// 강의의 tiled matmul 골격 (요약)
__global__ void matmul_tiled(float* A, float* B,
                              float* C, int M, int N, int K)
{
  __shared__ float shA[TILE][TILE];
  __shared__ float shB[TILE][TILE];

  int bi = blockIdx.y, bj = blockIdx.x;
  int ti = threadIdx.y, tj = threadIdx.x;
  int row = bi * TILE + ti;
  int col = bj * TILE + tj;

  float acc = 0.0f;
  for (int p = 0; p < (K + TILE - 1) / TILE; ++p) {
    // phase: load A tile, B tile
    int aCol = p * TILE + tj;
    int bRow = p * TILE + ti;
    shA[ti][tj] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f;
    shB[ti][tj] = (bRow < K && col < N) ? B[bRow * N + col] : 0.0f;
    __syncthreads();

    // inner: dot product on shared
    for (int k = 0; k < TILE; ++k)
      acc += shA[ti][k] * shB[k][tj];
    __syncthreads();
  }
  if (row < M && col < N) C[row * N + col] = acc;
}

강의의 학습 패턴은 L003 그대로. tiled matmul 도 Python 으로 먼저 흉내낸다. shared memory 는 Python dict 또는 mutable list 로, __syncthreads 는 “일단 외부 loop 가 모든 thread 의 step 을 마친 뒤 다음으로 넘어간다” 의 구조로 모델링.

# L005 notebook (요약) — shared 를 dict 로 mocking
def matmul_tiled_pyk(blockidx_x, blockidx_y, threadidx_x, threadidx_y,
                       blockdim, A, B, C, M, N, K):
    shared_A = torch.zeros(blockdim, blockdim)        # block 별로 따로
    shared_B = torch.zeros(blockdim, blockdim)

    bi, bj = blockidx_y, blockidx_x
    ti, tj = threadidx_y, threadidx_x
    row, col = bi*blockdim + ti, bj*blockdim + tj

    acc = 0.0
    for p in range(math.ceil(K / blockdim)):
        # load phase — shared 에 자기 한 element 채움
        aCol = p*blockdim + tj
        bRow = p*blockdim + ti
        shared_A[ti, tj] = A[row, aCol] if row < M and aCol < K else 0
        shared_B[ti, tj] = B[bRow, col] if bRow < K and col < N else 0
        # 외부 loop 가 모든 thread 를 한 step 끝낸 뒤 다음으로 가는 구조 = sync

        for k in range(blockdim):
            acc += shared_A[ti, k] * shared_B[k, tj]
    if row < M and col < N:
        C[row, col] = acc

그리고 이 Python 코드를 ChatGPT 또는 손으로 CUDA C 로 옮긴다 — 이미 본 mechanical 변환. 변수명을 그대로 유지하면 거의 한 글자도 안 바뀐다 (shared_A → shA 정도). 그게 이 강의의 학습 사다리의 가치.

__syncthreads() 는 block 안의 모든 thread 가 이 자리에 도달할 때까지 기다린다는 barrier. 다른 block 의 thread 와는 무관. shared memory 에 쓰고 다른 thread 가 읽기 전, 그리고 shared 의 한 자리를 다음 phase 가 덮어쓰기 전 — 이 두 자리가 표준.

강의에서 Jeremy 가 강조하는 한 줄 — “그래서 sync 의 위치는 ‘shared 에 write 한 직후, 그리고 shared 에서 read 끝낸 직후’ 두 자리.” 이 두 자리만 외워두면 거의 모든 tiled 커널이 race-free.

강의의 tiled matmul 코드는 boundary 가 있더라도 sync 자체는 모든 thread 가 통과한다 — boundary 는 load 의 값 만 0 으로 처리. 이 패턴이 CUDA tiled 커널의 표준.

정적 shared (__shared__ float tile[16][16]) 는 컴파일 시점에 크기가 박힌다. 강의에서 Jeremy 가 추가로 깐 dynamic shared memory 는 launch 시 크기를 결정한다 — 같은 커널이 여러 tile 사이즈에 대응 가능.

// dynamic shared 의 선언과 launch
__global__ void matmul_dyn(float* A, float* B,
                            float* C, int M, int N, int K, int TILE)
{
  extern __shared__ float buf[];           // 크기는 launch 시 지정
  float* shA = &buf[0];
  float* shB = &buf[TILE * TILE];      // 한 buffer 를 둘로 split
  // ... 이전과 동일
}

// host launch — 3번째 인자가 dynamic shared 크기 (bytes)
int sharedBytes = 2 * TILE * TILE * sizeof(float);
matmul_dyn<<<blocks, threads, sharedBytes>>>(
    A, B, C, M, N, K, TILE);

이 패턴의 장점 — 같은 커널 binary 가 여러 tile 사이즈에 대응. autotune 에서 tile 사이즈를 sweep 할 때 유용. 단점은 — compiler 가 사이즈를 모르니 register 할당 등 일부 최적화가 보수적.

강의의 한 가지 부수 메시지 — “최근 NVIDIA GPU 들은 모두 같은 양의 shared 를 가진다, 그래서 dynamic shared 의 portability 이점이 옛날만큼 크지 않다.” 정적 shared 가 종종 더 깔끔.

강의 후반의 작은 운영 팁. load_inline 의 name 을 함수마다 다르게 주면 — 각 함수가 자기 cache 디렉터리에 저장되어 다른 셀에서 바뀐 다른 함수를 빌드해도 이 함수는 재컴파일 안 함. Jeremy 가 L005 에서 추가한 패턴.

# L003 의 load_cuda 가 name="inline_ext" 고정이었던 것에서 변경
def load_cuda(cuda_src, cpp_src, funcs, name=None,
              opt=False, verbose=False):
    if name is None:
        name = funcs[0]                # 함수 이름을 module 이름으로
    return load_inline(
        cuda_sources=[cuda_src],
        cpp_sources=[cpp_src],
        functions=funcs,
        extra_cuda_cflags=["-O2"] if opt else [],
        verbose=verbose,
        name=name,
    )

이 작은 변경의 효과는 — naive matmul 커널과 tiled matmul 커널이 다른 cache 디렉터리에 저장. 두 커널을 번갈아 수정해도 한 쪽 수정이 다른 쪽 cache 를 invalidate 하지 않음. notebook 의 iterate 속도가 한 단 빨라진다.

강의에서 Jeremy 가 솔직하게 인정하는 한 자리. “우리 tiled matmul 을 cuBLAS / cuSPARSE / cuTLASS 와 직접 비교하지 않았다, 정확히 얼마나 차이나는지 모른다.” 이 한 줄이 “직접 짠 커널이 어디까지 갈 수 있는가” 의 솔직한 좌표.

위 표는 (1024, 1024, 1024) FP32 matmul 의 대략적 시간 (A100 기준 추정). 정확한 수치는 자기 GPU 에서 측정 필요. 메시지는 — shared memory 한 단을 추가하는 것만으로도 6배 가까운 가속이 가능하지만, cuBLAS 까지의 100배 격차에는 한참 못 미친다는 사실.

나머지 격차의 출처들.

• register tiling — 한 thread 가 shared 의 한 element 만 읽지 않고 여러 element 를 register 에 캐시. 한 thread 가 4×4 또는 8×8 의 출력 element 를 책임진다.
• vectorized loads — float4, __ldg 등으로 한 instruction 에 여러 element load. memory bandwidth 의 효율 ↑.
• double buffering / pipeline — 다음 phase 의 load 와 현재 phase 의 compute 를 overlap. cp.async (Hopper) 또는 software pipelining.
• Tensor Core — FP32 곱셈/덧셈을 MMA (matrix-multiply-accumulate) instruction 으로 묶음. 한 instruction 에 16×16 matmul. 별개 정량 영역.

강의의 결론적인 톤은 — “production 에서는 cuBLAS / cuTLASS 를 쓴다. 직접 짜는 건 학습용이거나, cuBLAS 가 잘 안 다루는 shape (예: tall-skinny matmul, fused attention) 때문에.” 그게 L012 FlashAttention, L023 CUTLASS 의 출발점.

tiled matmul 의 정수와 — 직접 손에 박아야 하는 코드 자료.

손에 새기기 — 실습 시퀀스

1. naive vs tiled 직접 비교 — matmul_l5.ipynb 그대로 돌려 (1024, 1024, 1024) 의 시간을 측정. naive 와 tiled 의 비율이 책에서 말하는 ~6배 와 일치하는지.
2. TILE 사이즈 sweep — TILE ∈ {8, 16, 32, 64} 로 측정. shared 사용량 변화와 시간 변화를 표로. occupancy 와의 관계 추정.
3. __syncthreads 일부러 빠뜨려보기 — 작은 사이즈 (16×16) 에서는 결과가 맞을 수 있지만, 큰 사이즈에서 race 가 발생하는지 확인.
4. boundary 처리 검증 — K 가 TILE 의 배수가 아닌 경우 (예: K=1000, TILE=16) 결과가 PyTorch 와 일치하는지 torch.allclose.
5. dynamic shared 로 바꿔보기 — 정적 shared 를 dynamic 으로 변환하고 launch 의 3번째 인자 추가. 결과는 같아야 함.
6. cuBLAS 와 비교 — torch.matmul 의 시간과 자기 tiled 커널 시간을 직접 비교. 격차가 강의의 추정과 비슷한지.
7. register tiling 의 한 단계 — 한 thread 가 2×2 출력 element 를 책임지게 변경. 시간이 더 줄어드는지.
8. 한 페이지 plan — 자기가 다루는 모델의 한 GEMM 호출을 골라 “이 shape 에서는 cuBLAS 가 답인가, 직접 짠 fused kernel 이 답인가” 의 결정 기준.

tiled matmul 의 패턴 (load → sync → compute → sync) 이 GPU 시리즈의 거의 모든 후속 강의에서 다른 옷을 입고 다시 등장한다.

이 학습 사다리가 의도적으로 비워둔 자리들과 — 자기 GPU 에서 직접 봐야 손에 박히는 사실들.

• cuBLAS 와의 정확한 격차 — Jeremy 가 강의에서 직접 측정 안 함. (확인 필요) 자기 GPU 에서 (1024, 1024, 1024) FP32 의 시간을 양쪽에서 재 비교.
• bank conflict — shared 가 32 bank 로 나뉘어 있고 같은 bank 를 동시 접근하면 직렬화. 이 강의의 shA[ti][k] 패턴이 bank conflict 를 일으키는가? 답은 stride 에 따라. L008 가 이 문제의 본론.
• register tiling 의 한 단 — 강의에서는 한 thread = 한 출력 element 의 패턴만. 한 thread 가 2×2 또는 4×4 element 를 책임지면 register reuse 가 늘어 더 빨라진다. 이건 L023 또는 자기 실험.
• Tensor Core — 강의에서 일절 다루지 않음. FP16/BF16 matmul + WMMA / MMA instruction 으로 가는 것이 production 의 길. L023 CUTLASS 강의에서.
• tile 사이즈와 occupancy 의 정확한 trade-off — 강의의 16×16 = 256 thread, shared 2×16×16×4 = 2 KB 가 어떤 occupancy 를 만들어 주는지의 정량은 빠짐.
• multi-GPU GEMM — 이 강의는 single GPU. tensor parallel 의 GEMM 분배 (예: AllReduce) 는 별도 강의 (L017 등).
• cp.async 와 double buffering — Hopper 의 async copy 가 다음 phase load 와 현재 compute 의 overlap 을 가능하게. 강의에서 다루지 않음. CUTLASS / FlashAttention 에서.