제3부 딥러닝 라이브러리

3.1 cuBLAS·cuBLASLt

3.1.1 이 장에서 배우는 것 / 선수 지식

이 장의 좁은 주제는 NVIDIA의 BLAS 라이브러리 — cuBLAS와 cuBLASLt — 를 직접 호출해 GEMM을 실행하고 튜닝하는 것이다. 2.7에서 직접 작성한 타일 행렬곱이 "안에서 무엇을 하는지"를 익혔다면, 이 장은 그 일을 production-grade로 대신 해 주는 라이브러리를 다룬다. PyTorch·TF의 matmul·Linear·Conv(im2col 경로)가 거의 전부 결국 이 라이브러리로 내려간다.

배우는 것:

• BLAS 레벨 1·2·3을 구분하고 왜 GEMM(레벨 3)이 GPU의 중심인지 연산-바이트 비율로 설명한다.
• cuBLAS의 열우선(column-major) 함정과 행우선을 곱하는 전치 트릭.
• cuBLAS SGEMM 호출을 작성·검증한다.
• cuBLASLt로 bias + GELU 에필로그 융합 GEMM을 작성한다.
• 핸들·스트림·워크스페이스·휴리스틱(autotuning)으로 워밍업·재사용 함정을 피한다.

선수 지식: 1.11~1.12(텐서코어·MMA), 1.13~1.14(정밀도·compute type vs data type), 2.7(타일 GEMM 골격), 2.10(스트림), 2.15(라이브러리·핸들 패턴). 환경은 0.3. 코드는 CUDA Toolkit 12.x 가정. 링크는 -lcublas(cuBLASLt는 -lcublasLt).

3.1.2 BLAS 레벨과 GEMM의 중심성

레벨 1·2는 데이터 대비 연산이 같은 차수라 비율이 상수 — HBM 대역폭에 묶인다(1.10 메모리 바운드). 레벨 3 GEMM만 데이터 $O(n^2)$ 대비 연산 $O(n^3)$ 이라 비율이 n에 비례, n이 크면 텐서코어를 포화시키는 연산 바운드가 된다. GPU 설계·라이브러리 최적화의 무게중심이 전부 GEMM에 쏠리는 이유.

거의 모든 딥러닝 연산은 GEMM이거나 GEMM으로 환원된다: Linear/FC($Y=X\cdot W^T$), 어텐션 QKᵀ·(·)V(배치 GEMM, 3.6), Conv(im2col 또는 implicit GEMM, 3.2). 그래서 "GPU를 빠르게 쓴다"는 실무적으로 "GEMM을 텐서코어로 포화시키고, 주변의 메모리 바운드 잡일(bias·activation·정규화)을 GEMM 커널에 융합한다"는 뜻이다 — 후자가 cuBLASLt 에필로그의 존재 이유.

3.1.3 cuBLAS API — 핸들·스트림·열우선

핸들과 수명주기. cuBLAS는 상태를 핸들(cublasHandle_t) 에 담는다 — 라이브러리 컨텍스트(스트림 바인딩·수학 모드·워크스페이스)를 든 대형 객체. 프로그램 시작에 한 번 만들고 끝에 한 번 파괴. 호출마다 만들면 안 된다(함정).

cublasHandle_t handle;
cublasCreate(&handle);            // 시작에 한 번
cublasSetStream(handle, stream);  // (선택) 비동기 스트림 바인딩 (2.10)
// ... 여러 번의 GEMM ...
cublasDestroy(handle);            // 끝에 한 번

열우선 레이아웃 — 이 장 최대 함정. cuBLAS는 Fortran 전통을 따라 열우선으로 행렬을 본다. C/C++·대부분의 프레임워크 텐서는 행우선이다. 행우선 $M\times N$ 배열을 cuBLAS에 그대로 넘기면 cuBLAS는 그것을 $N\times M$ 열우선, 즉 전치로 읽는다:

같은 메모리(a0…a5)를 C/PyTorch는 행우선 2×3으로, cuBLAS는 열우선 3×2(= 행우선의 전치)로 읽는다.

행우선을 곱하는 표준 트릭. 항등식 하나를 외운다: $(A\cdot B)^T = B^T\cdot A^T$. cuBLAS에 행우선 버퍼를 넘기면 전치로 읽히므로 cuBLAS가 보는 것은 $A^T$·$B^T$. 인자 순서를 B·A로 바꿔 호출하면 cuBLAS는 $B^T\cdot A^T = (A\cdot B)^T$를 열우선으로 쓰는데, 이것을 다시 행우선으로 읽으면 정확히 $A\cdot B$:

// 원하는 것(행우선): C(M×N) = A(M×K)·B(K×N)
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
    N, M, K,            // m=N, n=M, k=K (B를 첫 인자로)
    &alpha, B, N,       // B, ldb=N
    A, K, &beta,        // A, lda=K
    C, N);              // C, ldc=N → 행우선 C(M×N)

배치 GEMM. 독립적인 작은 GEMM 다발(어텐션 등)은 cublasSgemmStridedBatched(베이스 포인터 + stride, 거의 항상 선호)가 cublasSgemmBatched(포인터 배열)보다 낫다 — 포인터 목록을 따로 안 만들고 텐서가 [batch, M, K] 연속이면 stride = M*K 한 줄. 어텐션의 $Q\cdot K^T$가 전형(3.6).

3.1.4 워크스루 I — cuBLAS SGEMM (열우선 주의)

// sgemm_demo.cu  빌드: nvcc -O2 sgemm_demo.cu -lcublas -o sgemm_demo
#include <cublas_v2.h>
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define CK(x) do{if((x)!=cudaSuccess){printf("CUDA %s:%d\n",__FILE__,__LINE__);exit(1);}}while(0)
#define BK(x) do{if((x)!=CUBLAS_STATUS_SUCCESS){printf("cuBLAS %s:%d\n",__FILE__,__LINE__);exit(1);}}while(0)
int main(){
  const int M=256, K=192, N=128;           // 직사각으로 차원 혼동을 드러냄
  size_t szA=(size_t)M*K, szB=(size_t)K*N, szC=(size_t)M*N;
  float *hA=(float*)malloc(szA*4), *hB=(float*)malloc(szB*4), *hC=(float*)malloc(szC*4), *hRef=(float*)malloc(szC*4);
  for(size_t i=0;i<szA;i++) hA[i]=(float)((i%7)-3)*0.5f;
  for(size_t i=0;i<szB;i++) hB[i]=(float)((i%5)-2)*0.25f;
  for(int i=0;i<M;i++) for(int j=0;j<N;j++){ float s=0; for(int k=0;k<K;k++) s+=hA[i*K+k]*hB[k*N+j]; hRef[i*N+j]=s; }  // CPU 행우선
  float *dA,*dB,*dC; CK(cudaMalloc(&dA,szA*4)); CK(cudaMalloc(&dB,szB*4)); CK(cudaMalloc(&dC,szC*4));
  CK(cudaMemcpy(dA,hA,szA*4,cudaMemcpyHostToDevice)); CK(cudaMemcpy(dB,hB,szB*4,cudaMemcpyHostToDevice));
  cublasHandle_t h; BK(cublasCreate(&h));
  // BK(cublasSetMathMode(h, CUBLAS_TF32_TENSOR_OP_MATH)); // Ampere+: FP32→TF32 텐서코어(속도↑ 정밀도↓, 1.13)
  const float alpha=1.f, beta=0.f;
  BK(cublasSgemm(h, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, dB, N, dA, K, &beta, dC, N)); // (AB)ᵀ=BᵀAᵀ
  CK(cudaMemcpy(hC,dC,szC*4,cudaMemcpyDeviceToHost)); CK(cudaDeviceSynchronize());
  double maxerr=0; for(size_t i=0;i<szC;i++){ double e=fabs(hC[i]-hRef[i]); if(e>maxerr) maxerr=e; }
  printf("M=%d K=%d N=%d  maxerr=%.3e  %s\n", M,K,N, maxerr, maxerr<1e-3?"OK":"MISMATCH");
  cublasDestroy(h); return 0;
}

maxerr는 FP32 누적 순서 차이로 0이 아니라 작은 값(예 1e-4)이 정상. MISMATCH가 뜨면 거의 항상 레이아웃 트릭을 잘못 적용한 것이다.

일부러 틀려 보기: 트릭을 빼고 A, B 순서·M,N,K로 기본 설정 그대로 호출하면 — 직사각이면 차원 불일치 에러(다행), 정사각이면 조용히 틀린 값(MISMATCH만). 후자가 위험하다(빌드·실행은 되고 숫자만 틀림). 반드시 작은 직사각 + CPU 레퍼런스로 첫 호출을 검증.

3.1.5 cuBLASLt — 에필로그 융합

cuBLAS(legacy)는 고정 시그니처라 GEMM 본체만 한다. bias·activation·정규화는 별도 elementwise 커널을 다시 런치해야 하고, 그 커널들은 메모리 바운드라 GEMM 출력을 HBM에 썼다가 다시 읽는다 — 낭비. cuBLASLt는 디스크립터 기반으로 GEMM 직후 텐서코어 출력이 아직 레지스터·공유 메모리에 있을 때 후처리를 융합(epilogue fusion) 한다.

융합이 줄이는 것은 FLOP가 아니라 HBM 왕복이다. activation·bias는 메모리 바운드(3.1.2)라 GEMM 커널에 융합하면 사실상 무비용에 가깝다. 프레임워크 Linear+GELU가 cuBLASLt 에필로그로 내려가는 이유.

디스크립터 4종 (분리가 유연성의 근원): cublasLtMatmulDesc_t(operation — compute type·scale·에필로그 종류·bias 포인터·전치), cublasLtMatrixLayout_t×3(A·B·C/D — 타입·차원·ld·레이아웃 순서·배치 stride), cublasLtMatmulPreference_t(워크스페이스 상한), cublasLtMatmulHeuristicResult_t(휴리스틱이 고른 알고리즘·타일). 에필로그는 operation 속성: CUBLASLT_EPILOGUE_BIAS·_RELU·_GELU·_GELU_BIAS(둘 다)·_GELU_AUX(역전파 중간값) 등.

휴리스틱·autotuning. cublasLtMatmulAlgoGetHeuristic이 후보 알고리즘(타일·스테이지·split-K·텐서코어 instruction)을 점수 순으로 돌려준다. ① 1순위를 그대로 쓰거나(대부분 충분), ② 상위 K개를 실측해 가장 빠른 것을 고정(형상 고정·반복 시 캐시), ③ 직접 지정(재현성). 텐서코어 활용은 이 안에서 결정 — 타입이 TF32·FP16·BF16·FP8이고 차원이 텐서코어 타일에 맞으면(보통 K·M·N이 8·16의 배수) 텐서코어 알고리즘을 고른다. 차원이 16의 배수가 아니면 텐서코어를 못 쓸 수 있어, 모델 hidden size를 8·16 배수로 설정하는 것이 실무 관행.

3.1.6 워크스루 II — cuBLASLt로 bias+GELU 융합 GEMM

D = GELU(A·B + bias)를 한 번의 cuBLASLt matmul로. 행우선 입력을 CUBLASLT_ORDER_ROW로 처리해 전치 트릭 없이.

// ltgelu_demo.cu  빌드: nvcc -O2 ltgelu_demo.cu -lcublasLt -o ltgelu_demo  (핵심 부분만)
cublasLtHandle_t lt; cublasLtCreate(&lt);
// 1) operation: FP32 누적, 에필로그 = GELU_BIAS, bias 포인터
cublasLtMatmulDesc_t op; cublasLtMatmulDescCreate(&op, CUBLAS_COMPUTE_32F, CUDA_R_32F);
cublasLtEpilogue_t epi = CUBLASLT_EPILOGUE_GELU_BIAS;
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_EPILOGUE, &epi, sizeof(epi));
cublasLtMatmulDescSetAttribute(op, CUBLASLT_MATMUL_DESC_BIAS_POINTER, &dBias, sizeof(dBias));
// 2) 레이아웃: 행우선 A(M×K), B(K×N), D(M×N)
cublasLtMatrixLayout_t la,lb,ld;
cublasLtMatrixLayoutCreate(&la, CUDA_R_32F, M, K, K);   // ld=K (row-major)
cublasLtMatrixLayoutCreate(&lb, CUDA_R_32F, K, N, N);
cublasLtMatrixLayoutCreate(&ld, CUDA_R_32F, M, N, N);
cublasLtOrder_t row = CUBLASLT_ORDER_ROW;
for (auto L : {la,lb,ld}) cublasLtMatrixLayoutSetAttribute(L, CUBLASLT_MATRIX_LAYOUT_ORDER, &row, sizeof(row));
// 3) 워크스페이스 + 휴리스틱
size_t wsBytes = 32*1024*1024; void* ws; cudaMalloc(&ws, wsBytes);
cublasLtMatmulPreference_t pref; cublasLtMatmulPreferenceCreate(&pref);
cublasLtMatmulPreferenceSetAttribute(pref, CUBLASLT_MATMUL_PREF_MAX_WORKSPACE_BYTES, &wsBytes, sizeof(wsBytes));
cublasLtMatmulHeuristicResult_t heur[8]; int got=0;
cublasLtMatmulAlgoGetHeuristic(lt, op, la, lb, ld, ld, pref, 8, heur, &got);
// 4) 실행: D = GELU(alpha·A·B + bias)  (beta=0이라 C 미사용, C와 D를 같은 dD로)
const float alpha=1.f, beta=0.f;
cublasLtMatmul(lt, op, &alpha, dA, la, dB, lb, &beta, dD, ld, dD, ld, &heur[0].algo, ws, wsBytes, 0);

주의: maxerr 임계를 SGEMM(1e-3)보다 느슨한 1e-2로 두는 이유 — cuBLASLt의 GELU와 CPU의 tanh 근사 GELU가 정확히 같은 공식이 아닐 수 있다(라이브러리가 erf 기반 정확 GELU를 쓸 수도, 버전 의존). 근사 일치로 검증하라(정직성 포인트). 같은 op·layout·algo로 형상이 같은 GEMM을 반복 호출하면 디스크립터·휴리스틱 비용이 한 번으로 분산 — 학습 루프 권장 패턴.

GEMM 타일 + 에필로그 융합(2.7의 확장). legacy라면 MMA 누적 직후 HBM에 쓰고 별도 커널이 다시 읽는다.

3.1.7 중간 점검 (체크포인트)

• [ ] 레벨 1·2 vs 3을 연산-바이트 비율로 구분하고 "왜 GEMM만 연산 바운드가 될 수 있나"를 제시한다.
• [ ] 행우선 버퍼를 cuBLAS에 넘길 때 인자를 B·A 순서로, m·n에 N·M($(AB)^T=B^TA^T$)으로 두는 이유를 제시한다.
• [ ] cuBLASLt 디스크립터 4종의 역할과 에필로그가 operation 디스크립터의 속성임을 안다.
• [ ] 융합이 FLOP가 아니라 HBM 왕복을 줄인다는 점, 텐서코어를 위해 차원을 8·16 배수로 맞추는 관행을 안다.

증상: MISMATCH(정사각인데 틀림) → 열우선 트릭 오적용(직사각+CPU 레퍼런스로 격리). cuBLASLt got=0 → 휴리스틱이 후보를 못 찾음(워크스페이스 상한↑·에필로그/타입/차원 조합 확인). 첫 호출만 느림 → 정상(워밍업 후 측정).

3.1.8 트레이드오프

견해: 대부분 cuBLAS로 시작해 융합·정밀도가 필요해지면 cuBLASLt, 라이브러리가 못 주는 융합이 있을 때만 CUTLASS. 비NVIDIA: CPU OpenBLAS·MKL(비교 기준선), AMD rocBLAS·hipBLASLt(API 의도적 유사). cuBLAS의 락인은 API보다 텐서코어·세대별 튜닝 커널의 성숙도 — NVIDIA 스택의 실질 해자(0.1).

3.1.9 흔한 함정

• 행/열우선 혼동(최대). 정사각에서 트릭을 빼먹으면 빌드·실행은 되고 값만 틀린다. 작은 직사각 + CPU 레퍼런스로 첫 호출 검증.
• 핸들을 호출마다 생성·파괴. cublasCreate는 부하가 크다. 한 번 만들어 재사용.
• 워밍업 없이 측정. 첫 호출은 라이브러리 초기화·커널 JIT·휴리스틱 비용 포함. 같은 형상으로 워밍업하고 cudaDeviceSynchronize 후 측정(2.13).
• 정밀도 혼동. computeType(누적)과 데이터 타입(저장)을 헷갈림. TF32 math mode를 켜면 FP32 입력이 TF32로 잘려 속도↑ 정밀도↓(1.13) — 의도한 것이 아니면 결과가 미묘하게 바뀐다.
• 워크스페이스 부족. 너무 작게 주면 빠른 split-K 알고리즘이 후보에서 빠져 느려지거나 got=0. 보통 수십 MB.
• 스트림 미바인딩. cublasSetStream 없으면 default stream에서 직렬화(2.10).

3.1.10 연습문제

연습 1 (레벨·비율). BLAS 레벨 1·2·3 각각의 연산-바이트 비율을 차원으로 표현하고 "왜 GPU에서 레벨 3만 연산 바운드가 될 수 있나"를 한 문장으로.

연습 2 (틀린 호출). sgemm_demo.cu를 3.1.4의 "틀린 호출"(A,B 순서·M,N,K)로 바꿔 빌드·실행하라. 직사각과 정사각(M=K=N=256)에서 차이는?

연습 3 (TF32 정밀도). cublasSetMathMode(h, CUBLAS_TF32_TENSOR_OP_MATH)를 켜고 maxerr를 비교하라(Ampere+). 왜 커지나?

연습 4 (융합 측정). ltgelu_demo의 융합 버전과, 같은 결과를 cuBLAS SGEMM + 별도 bias 커널 + 별도 GELU 커널로 분리한 버전을 Nsight Systems(2.13)로 프로파일하라. 커널 개수·HBM 트래픽 차이는? 어떤 형상에서 융합 이득이 큰가?

3.1.11 요약 · 다음 장

• BLAS는 레벨 1·2(메모리 바운드)·3(GEMM, 연산 바운드 가능)으로 나뉘고, 딥러닝 FLOP 대부분이 GEMM이라 최적화가 전부 여기 쏠린다.
• cuBLAS 입문 최대 함정은 열우선 레이아웃 — 행우선은 $(AB)^T=B^TA^T$ 트릭(B·A 순서)으로, 직사각+CPU 레퍼런스로 첫 호출 검증.
• cuBLASLt는 디스크립터 기반으로 에필로그 융합·정밀도 혼합·autotuning을 준다. 융합은 FLOP가 아니라 HBM 왕복을 줄인다.
• 텐서코어 활용은 휴리스틱이 타입·차원에 따라 자동 선택 — 차원을 8·16 배수로. 함정은 행/열우선·정밀도·핸들 재사용·워밍업·워크스페이스에 몰려 있다.

다음 장: 3.2 cuDNN I — 프리미티브. GEMM 너머로, 컨볼루션·풀링·정규화·활성을 다루는 cuDNN과, 컨볼루션을 계산하는 여러 알고리즘(implicit GEMM·Winograd·FFT)을 본다.