[Multi-Core & GPU Programming] CUDA - Matrix Multiplication

"CUDA - Matrix Multiplication"

[CUDA Performance Considerations]

> Setting # of Threads

· 성능 최적화를 위해 Thread 수를 적절히 설정 → SM당 할당할 수 있는 최대 Thread 수를 파악하고 확보

> Branch Divergence (△)

· Thread 간의 분기 실행이 달라지는 Branch Divergence를 고려 → HW Resource의 활용도 저하를 유발하는 원인으로 피해야 함

> Shared Memory (O)

· 수동으로 관리할 수 있는 Cache 형태의 메모리를 사용하여 데이터의 Locality를 포착

· 병렬 접근 시 속도를 저하시키는 Bank Conflict를 피해야 함

> Coalescing(메모리 병합) (O)

· Thread들의 메모리 접근을 정렬하여 메모리 대역폭의 활용을 극대화하는 최적화 방법

> Overlap I/O with Computation

· I/O 작업과 연산 처리를 동시에 진행하여 대기 시간을 감소 → 이를 구현하기 위해 CUDA의 "Stream" 기능 활용

** GPU Programming에서의 I/O는 CPU와 GPU Memory 간의 데이터 전송을 의미

▣ CUDA Memory Model

· Register : Thread 당 할당되며 가장 빠름 (~ 1 Cycle)

· Shared Memory : Block 내 Thread들이 공유하며 매우 빠름 (~ 5 Cycle)

· Global Memory : Grid 내 모든 Thread가 접근 가능하지만 Chip 외부에 있어 매우 느림 (~ 500 Cycle)

· Constant Memory : Grid 내 Read 전용이며 Caching 될 경우 빠름 (~ 5 Cycle)

∴ CUDA 최적화 : Global Memory 접근을 줄이고, Register 또는 Shared Memory처럼 빠른 메모리를 최대한 활용

▣ CUDA Variable Type Qualifiers

· __shared__ : Shared Memory에 저장되어 Block 안의 Thread들이 공유

· __device__ : Global Memory에 위치하고 Grid 전체에서 접근 가능

· __constant__ : Constant Memory에 저장되어 Read 전용으로 사용

** 일반 지역 변수는 Register에 저장 / Thread 별 배열은 Global Memory에 저장

[Matrix Multiplication]

1 2 3 4

for i = 1 to n     for j = 1 to n         for k = 1 to n             P(i, j) = P(i, j) + M(i, k) * N(k, j)

cs

· P(i, j)가 독립적이므로 GPU 병렬화에 적합

1 2

for k = 1 to n     P(i, j) = P(i, j) + M(i, k) * N(k, j)

cs

· 각 Thread가 결과 행렬 P의 원소 하나(P_{i, j})를 계산하도록 배정 & blockIdx와 threadIdx를 사용하여 행렬의 i, j Index를 결정

→ But, 이 단순한 구현은 매 계산마다 Global Memory에 자주 접근

    ∵ 각 Thread가 M과 N의 원소를 반복적으로 Global Memory에서 Read → Bandwidth & Latency 낭비

· GTX1080의 L1 Cache는 48KB, SM 하나는 최대 2048개의 Thread를 실행 가능

→ Thread당 할당되는 Cache는 약 24B로 매우 작음

→ Spatial Locality : 활용 가능(∵ 같은 Cache Line 안의 연속 주소 접근)

→ Temporal Locality : 활용 어려움(∵ Cache가 너무 작고 Thread가 너무 많아 나중에 다시 쓸 데이터 보존 불가)

∴ 반복 재사용을 Cache에만 의존하기 어려움

1 2 3

for k = 1 to n     sum += M(i, k) * N(k, j) P(i, j) = sum

cs

· Global Memory인 P(i, j)에 직접 계속 더하는 대신, Local 변수(Register)인 "sum"을 사용

· Loop마다 "sum"을 통해 합계를 누적한 후, 마지막에 한 번 Global Memory에 기록

→ 더 많은 Register를 사용함으로써 Global Memory Bandwidth 사용량 감소 가능

▣ 2D Thread Block

· 4x4 결과 행렬 P → Block(2x2 행렬) x 4

· 각 Block은 4(= 2x2)개의 Thread를 가지고, 각 Thread는 결과 행렬의 원소 하나를 계산

** Thread Block 크기가 Warp 크기인 32보다 충분히 커야 효율적

▣ Kernel 호출

· dimBlock : 한 Block 안의 Thread 배열 크기

· dimGrid : 전체 Block 배열 크기

· (Width) x (Width) 행렬 & (TILE_WIDTH) x (TILE_WIDTH) → Grid는 각 방향으로 (Width) / (TILE_WIDTH)개의 Block 보유

    (Ex. 4x4 행렬 & 2x2 Tile → Row & Column 방향으로 2개의 Block 보유)

▣ 가장 단순한 CUDA 행렬 곱셈 Kernel

① 계산할 위치를 정하기 위해 row, col 정의

② k = 0 ~ (width - 1)까지 반복하면서 M[row][k] x N[k][col] 값을 "value"에 누적

③ 마지막에 P[row][col] = value로 결과를 Global Memory에 저장

→ 단순하지만, M과 N을 매번 Global Memory에서 Read 하므로 효율↓

▣ Work Flow

col = blockIdx.x * blockDim.x + threadIdx.x

row = blockIdx.y * blockDim.y + threadIdx.y

> Block(0, 0) - P(0, 0), P(0, 1), P(1, 0), P(1, 1)

· blockIdx.x = 0 → col = 0 * 2 + threadIdx.x

· blockIdx.y = 0 → row = 0 * 2 + threadIdx.y

→ 각 Thread는 자기 위치의 P 원소 하나를 계산

> Block(0, 1) - P(0, 2), P(0, 3), P(1, 2), P(1, 3)

· blockIdx.x = 1 → col = 1 * 2 + threadIdx.x

· blockIdx.y = 0 → row = 0 * 2 + threadIdx.y

→ 각 Thread는 자기 위치의 P 원소 하나를 계산

[Shared Memory]

· N x N 행렬 곱셈에서 각 데이터 원소는 총 N번씩 반복 사용

· But, L1 Cache가 작아서 Temporal Locality X → 같은 값을 다시 사용할 때마다 Thread가 Global Memory에서 다시 Read

▣ Shared Memory

· Shared Memory : 같은 Thread Block 안의 Thread들이 함께 사용하는 빠른 On-Chip Memory

→ Thread Block 내부에서 공유 : 같은 Thread Block에 속한 Thread들끼리 데이터 공유

→ 고성능 : Global Memory보다 훨씬 낮은 Latency와 높은 Memory Bandwidth

** L1 Cache와 달리 Programmer가 직접 제어하는 Scratchpad Memory

> Programming 방법

① 데이터를 Shared Memory에 들어갈 수 있는 작은 Subset으로 나눔

② 해당 Subset을 하나의 Thread Block이 처리

③ 필요하면 Shared Memory에서 계산한 결과를 Global Memory로 복사

· M과 N의 일부 Tile을 Global Memory에서 바로 Read X

· Block 안의 Thread들이 함께 쓸 수 있는 Shared Memory 배열 subTileM, subTileN에 복사

→ 같은 Tile 데이터를 여러 Thread가 빠르게 재사용

> Shared Memory Idea

· 여러 Thread가 Global Memory에서 Read → Global Memory 접근이 많아 느림

① Global Memory의 일부 데이터를 On-Chip Memory(Shared Memory)로 복사

② 여러 Thread가 Shared Memory에서 데이터를 Read

△ Global Memory 접근 횟수 감소

△ 같은 데이터를 Thread Block 내부에서 여러 번 재사용

▣ Work Flow - Compute Block(0, 0)

> Tile 1

① Global Memory에 있는 큰 행렬 M, N에서 1번째 Tile의 데이터를 Shared Memory로 Load

② "__syncthreads( )"로 동기화 수행

    ∵ 모든 Thread가 각자의 담당 데이터를 Shared Memory에 다 복사기 전에 계산을 하면 잘못된 값 Read 가능
③ Shared Memory에 Load된 Tile1을 계산

④ "__syncthreads( )"로 동기화 수행

    ∵ 어떤 Thread가 아직 Shared Memory 값을 사용 중인데 다른 Thread가 다음 Tile을 덮어쓰면 안됨

· P(0, 0) = M(0, 0) x N(0, 0) + M(0, 1) x N(1, 0)

· P(0, 1) = M(0, 0) x N(0, 1) + M(0, 1) x N(1, 1)

· P(1, 0) = M(1, 0) x N(0, 0) + M(1, 1) x N(1, 0)

· P(1, 1) = M(1, 0) x N(0, 1) + M(1, 1) x N(1, 1)

> Tile 2

① Global Memory에 있는 큰 행렬 M, N에서 2번째 Tile의 데이터를 Shared Memory로 Load

② "__syncthreads( )"로 동기화 수행

③ Shared Memory에 Load된 Tile2를 계산

④ "__syncthreads( )"로 동기화 수행

· P(0, 0) = M(0, 2) x N(2, 0) + M(0, 3) x N(3, 0)

· P(0, 1) = M(0, 2) x N(2, 1) + M(0, 3) x N(3, 1)

· P(1, 0) = M(1, 2) x N(2, 0) + M(1, 3) x N(3, 0)

· P(1, 1) = M(1, 2) x N(2, 1) + M(1, 3) x N(3, 1)

> Block(0, 0)

· P(0, 0) = {M(0, 0) x N(0, 0) + M(0, 1) x N(1, 0)} + {M(0, 2) x N(2, 0) + M(0, 3) x N(3, 0)}

· P(0, 1) = {M(0, 0) x N(0, 1) + M(0, 1) x N(1, 1)} + {M(0, 2) x N(2, 1) + M(0, 3) x N(3, 1)}

· P(1, 0) = {M(1, 0) x N(0, 0) + M(1, 1) x N(1, 0)} + {M(1, 2) x N(2, 0) + M(1, 3) x N(3, 0)}

· P(1, 1) = {M(1, 0) x N(0, 1) + M(1, 1) x N(1, 1)} + {M(1, 2) x N(2, 1) + M(1, 3) x N(3, 1)}

▣ Tiled Matrix Multiplication Kernel

1 ~ 2 : Shared Memory Tile 2개를 선언

5 ~ 6 : 현재 Thread가 계산할 Row, Column 정의

7 : 누적에 사용할 Local 변수 "Pvalue"를 0으로 초기화

8 ~ 10 : m Loop를 돌면서 필요한 Tile들을 하나씩 Shared Memory에 Load (동기화를 위해 Load 완료 대기)

                → 각 Thread는 M에서 원소 하나, N에서 원소 하나를 가져와서 Shared Memory에 저장

12 ~ 14 : k Loop를 돌면서 Shared Memory의 Tile을 이용해 곱셈 결과를 Pvalue에 누적(동기화를 위해 계산 완료 대기)

16 : 마지막 누적 값을 Global Memory에 저장

> __syncthreads( )

· 같은 Block 안의 모든 Thread가 "__syncthreads( )" 지점에 도달해야만 그 다음 코드로 넘어갈 수 있음

① Tile의 모든 원소가 Shared Memory에 모두 Load되었는지 보장

② 다음 Tile로 Shared Memory를 덮어쓰기 전에 기존 Tile 계산이 모두 끝났는지 보장

** "__syncthreads( )"는 같은 Thread Block 안에서만 작동 = 다른 Block과 동기화 X

    → Block(0, 0) 안의 Thread들은 서로 기다림, But, Block(0, 0)과 Block(1, 0)이 서로 기다림 X

▣ Memory 분석

· (b x b) 크기의 Shared Memory Tile 사용 시 Tile 안의 각 원소가 b번 재사용

→ Tile 크기가 커질수록 Global Memory에서 한 번 읽은 값을 더 많이 활용 가능

· 하나의 Tile을 계산하려면 (b x b) 크기의 Tile을 N/b번 Load

· 전체 Thread Block 수 : (열 방향 Tile 개수) x (행 방향 Tile 개수) = (N/b)²

· 전체 Global Memory Load 수 : N/b x b² x (N/b)² = N³/b

→ (1 Thread Block 반복 횟수) x (각 반복마다 Load하는 Element 수) x (전체 Thread Block 수)

→ 단순 구현 : N³ → 한 번 가져온 Tile Element를 b번 재사용 → b배만큼 Memory Load 감소

▣ Tile Size 결정

· SM 당 Shared Memory : 96KB

· 행렬 당 Shared Memory : 48KB (∵ M, N Tile 2개)

· 원소 개수 : 12,288개 (∵ 48KB / 4B(float))

→ Block : 110x110 (∵ √12,288 ≒ 110) = 오류

∵ Shared Memory는 한 Block만 쓰는 것이 아니라 같은 SM 안에서 동시에 실행되는 여러 Thread Block들이 나눠 써야 함

∵ 110x100 → Thread 수 = 12100개 → CUDA의 최대 Thread Block 크기인 1024개를 훨씬 초과

> 8x8 Tile Size 결정

· 8x8 Tile → Block 하나에 Thread가 64개 → 2 Warp(∵ Warp 하나당 Max Thread Block = 32)

· Block 하나당 8 x 8 x 2(M, N Tile) x 4B(float) = 512B의 Shared Memory 사용

· SM의 최대 Thread : 2048 → 64(Thread) x 32(Block) = 2048 Threads → 활용도 100% 가능

· Shared Memory 사용량 : 512B x 32 = 16KB (< 96KB)

> 16x16 Tile Size 결정

· 16x16 Tile → Block 하나에 Thread가 256개 → 8 Warp(∵ Warp 하나당 Max Thread Block = 32)

· Block 하나당 16 x 16 x 2(M, N Tile) x 4B(float) = 2048B(= 2KB)의 Shared Memory 사용

· SM의 최대 Thread : 2048 → 256(Thread) x 8(Block) = 2048 Threads → 활용도 100% 가능

· Shared Memory 사용량 : 2KB x 8 = 16KB (< 96KB)

> 32x32 Tile Size 결정

· 32x32 Tile → Block 하나에 Thread가 1024개 → 32 Warp(∵ Warp 하나당 Max Thread Block = 32)

· Block 하나당 32 x 32 x 2(M, N Tile) x 4B(float) = 8192B(= 8KB)의 Shared Memory 사용

· SM의 최대 Thread : 2048 → 1024(Thread) x 2(Block) = 2048 Threads → 활용도 100% 가능

· Shared Memory 사용량 : 8KB x 2 = 16KB (< 96KB)

** 8x8, 16x16, 32x32 모두 SM 전체 Shared Memory 사용량이 16KB로 동일

    → Tile을 키우면 Block당 Shared Memory는 증가하지만, 동시에 Block 수는 감소 & Thread 수 제한 존재

∴ 남은 Shared Memory를 더 활용하려면 Thread 하나가 결과 원소 하나만 계산하지 않고, 여러 원소를 계산하도록 Kernel 수정

    → Thread 당 Work를 늘리면 더 큰 Tile을 처리하거나 Shared Memory를 더 적극적으로 활용 가능

[Coalescing]

▣ Coalescing

M[row * width + k] * N[k * width + col]

> Single Thread : 1개의 Thread만 보면 Row-Major 방식이므로 연속 접근 → 행렬 M이 빠름

· t0 : M[0][0] → M[0][1] → M[0][2] → M[0][3]

· t1 : M[1][0] → M[1][1] → M[1][2] → M[1][3]

· t2 : M[2][0] → M[2][1] → M[2][2] → M[2][3]

· t3 : M[3][0] → M[3][1] → M[3][2] → M[3][3]

<M : No Coalescing>

> Multi Thread : Warp 전체를 보면 M은 모두 같은 값을 Read, N은 연속된 주소에 동시 접근 → 행렬 N이 빠름

· k = 0 : t0 = M[0][0], t1 = M[0][0], t2 = M[0][0], t3 = M[0][0]

· k = 1 : t0 = M[0][1], t1 = M[0][1], t2 = M[0][1], t3 = M[0][1]

· k = 2 : t0 = M[0][2], t1 = M[0][2], t2 = M[0][2], t3 = M[0][2]

· k = 3 : t0 = M[0][3], t1 = M[0][3], t2 = M[0][3], t3 = M[0][3]

<N : Coalescing>

> Multi Thread : Warp 전체를 보면 M은 모두 같은 값을 Read, N은 연속된 주소에 동시 접근 → 행렬 N이 빠름

· k = 0 : t0 = N[0][0], t1 = N[0][1], t2 = N[0][2], t3 = N[0][3]

· k = 1 : t0 = N[1][0], t1 = N[1][1], t2 = N[1][2], t3 = N[1][3]

· k = 2 : t0 = N[2][0], t1 = N[2][1], t2 = N[2][2], t3 = N[2][3]

· k = 3 : t0 = N[3][0], t1 = N[3][1], t2 = N[3][2], t3 = N[3][3]

** 만약 같은 Warp 안의 Thread들이 같은 Cache Line 안의 연속된 주소를 Read

    → GPU는 이를 하나의 Memory Transaction으로 묶음 가능

▣ Corner Turning

· Tile을 Global Memory에서 Shared Memory로 Load할 때 Thread들이 연속 주소를 Read하도록 방식 조정

    → Coalesced Pattern으로 데이터를 빠르게 가져옴

    → Shared Memory 안에서는 필요한 방향으로 Read

▣ Tiled Shared Memory Kernel

· 9 : "M" Load에서는 같은 "ty"를 가진 Thread들이 "tx"만 달라지며 연속된 "M" 원소 Read → Coalescing

· 10 : "N" Load에서도 같은 "ty"에서 "tx"가 달라지면 연속된 "N" 원소 Read → Coalescing

△ Global Memory 접근 횟수 감소

△ Global Memory 접근을 Coalesced Pattern으로 변경 (단순 Kernel보다 훨씬 빠르게 동작 가능)

[Reference]

· 5_cuda_matmul_mgp_2026 (MGP) - Yongjun Park