[Multi-Core & GPU Programming] Prefix Sum

"Prefix Sum"

[Sum Reduction]

▣ Sum Reduction Problem
· "worker" 함수는 배열의 일부를 받아서 "output"에 계속 더함
· 모든 Thread가 같은 "output"(= count)에 동시에 접근하는 문제 존재
→ Thread A, B가 동시에 "count += input[i]" 실행 시, 둘 중 하나의 결과가 사라질 수 있음 = Race Condition

> Try 1 - Lock
· "global_mutex.lock( )"과 "unlock( )"으로 공유 변수 "output"을 보호
· 한 번에 하나의 Thread만 "output"을 수정하므로 올바른 결과 출력
→ 모든 덧셈마다 Lock을 잡고 풀어야 하기 때문에 핵심 연산은 거의 직렬처럼 동작하여 성능↓

> Try 2 - Atomic
· "std::atomic<int>"를 사용하여 짧은 연산에서는 Lock 방식보다 더 효율적
→ 읽기 작업에만 약간의 병렬화가 적용되지만, 핵심 연산은 여전히 직렬처럼 동작하여 성능↓

▣ Serial Reduction
· 앞에서부터 차례로 더하는 방식
· N개의 숫자를 1개의 값으로 줄이려면 덧셈 연산 (N - 1)번 필요
→ 연산량 : O(N) / 실행 시간 : O(N)

▣ Parallel Reduction
· 데이터를 2개씩 짝지어 Tree 구조로 더하는 방식
· N개의 숫자를 1개의 값으로 줄이려면 덧셈 연산 (N - 1)번 필요
→ 연산량 : O(N) / 실행 시간 : O(logN)

** Serial Reduction 방식과 연산량은 같지만, 병렬 실행 덕분에 실행 시간↓

▣ Prefix Sum(Scan)
· 각 위치까지의 누적합을 구하는 연산
· Quick Sort(Partition), Histogram, Evaluate Polynomials, Radix Sort, Search for Regular Expressions 등에 활용
· N개의 숫자를 1개의 값으로 줄이려면 덧셈 연산 (N - 1)번 필요
→ 연산량 : O(N) / 실행 시간 : O(logN)

[Parallel Prefix Sum]

▣ Naive Prefix Sum
· 각 Thread가 자기 위치까지의 합을 직접 계산
  - Thread 0 : 1
  - Thread 1 : 1 + 2
  - Thread 2 : 1 + 2 + 3
  - Thread k : input[0] + … + input[k]
→ 전체 연산량 : O(N2)
→ 병렬화 자체는 쉽지만, Thread마다 중복 계산이 너무 많기 때문에 성능↓

▣ Kogge-Stone Algorithm
· 이웃한 값들을 더하고, 각 단계마다 Barrier 배치

> sum[4] : Index 4의 Prefix Sum
· Stride = 1 : 1칸 앞의 값과 덧셈
· Stride = 2 : 이전 단계의 결과와 2칸 앞에 덧셈 연산된 값과 덧셈
· Stride = 4 : 이전 단계의 결과와 4칸 앞에 덧셈 연산된 값과 덧셈
→ Index 4는 자기 앞의 모든 값을 포함

> sum[5] : Index 5의 Prefix Sum

> sum[6] : Index 6의 Prefix Sum

> sum[7] : Index 7의 Prefix Sum

> 전체 배열에 대한 Kogge-Stone Algorithm

· k번째 반복에서의 계산 식 : Ak+1[i] = Ak[i - 2k] + Ak[i] (Stride = 2k)
→ 전체 실행 시간 : O(logN)

> Kogge-Stone Algorithm의 한계
· 현실에는 무한한 Core X
· Kogge-Stone Algorithm은 단계 수가 적지만, 각 단계마다 많은 연산 진행
    → (N - 1) + (N - 2) + (N - 4) + … + (N - N/2) = O(N logN)
· P개의 Core : 실행 시간 : O((N/P) logN)
· 각 단계의 값을 저장하려면 Memory도 많이 필요해서, 단순 구현에서는 logN배 추가 Memory가 필요할 수 있음
 
△ 빠른 병렬 시간
▼ 연산량과 Memory 사용량↑

** Double Buffering을 사용하여, 모든 단계의 결과를 따로 저장하지 않고, 2개의 Buffer만 번갈아 사용
→ Memory Overhead : logN배 → 2배(현재 단계 값을 읽는 Buffer / 다음 단계 결과를 쓰는 Buffer)
→ But, 연산량 자체가 O(N logN)인 문제는 여전히 존재

▣ Brent-Kung Algorithm
· Kogge-Stone Algorithm보다 연산량을 줄이기 위한 Prefix Sum Algorithm
· Balanced Binary Tree 형태의 계산 패턴 사용

> Phase 1
· Parallel Reduction와 동일하게 전체 합을 구하는 Reduction 과정(Bottom-Up 방식)

> Phase 2
· 전체 배열 앞에는 아무 값도 없으니깐 시작 Prefix는 0부터 시작(Top-Down 방식)

· Prefix[Left(v)] = Prefix[v] / Prefix[Right(v)] = Prefix[v] + Sum[Left(v)]

· Prefix[Left(v)] = Prefix[v] / Prefix[Right(v)] = Prefix[v] + Sum[Left(v)]

· Prefix[Left(v)] = Prefix[v] / Prefix[Right(v)] = Prefix[v] + Sum[Left(v)]

· Exclusive Scan = [0, 1, 3, 6, 10, 15, 21, 28] : 현재 원소를 포함하지 않고, 자기 앞 원소들의 합만 저장
· Inclusive Scan = [1, 3, 6, 10, 15, 21, 28, 36] : 각 위치에 자기 자신의 값을 더함

▶ Brent-Kung Algorithm - Array View
> Phase 1
· Parallel Reduction처럼 부분합 생성

△ Extra Memory 필요 X (별도의 Tree 자료 구조를 만들지 않고 입력 배열 자체를 갱신하면서 계산 가능)

> Phase 2
· Prefix 값을 아래로 전파 (왼쪽 : Prefix[v] / 오른쪽 : Sum[Left(v)])

> Complexity
· Phase 1 : 아래에서 위로 Reduction (실행 시간 : logN / 연산량 : N - 1)
· Phase 2 : 위에서 아래로 Prefix 전파 (실행 시간 : logN / 연산량 : N - 1)
→ 전체 실행 시간 : 2 logN = O(logN) / 전체 연산량 : 2(N - 1) = O(N)

· Kogge-Stone보다 단계수는 조금 많지만, 전체 연산량은 훨씬 적음 & P개의 Processor : O(N/P) 수준으로 효율적

[GPU Parallel Prefix Sum]

▣ Naive Prefix Sum
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
__global__ void scan(float *g_odata, float *g_idata, int n) {
    extern __shared__ float temp[];        // Allocated On Invocation
    int tid = threadIdx.x;
    int pout = 0, pin = 1;
 
    // Load Input Into Shared Memory
    // This is exclusive scan, so shift right by one and set first elt to 0
    temp[pout * n + tid] = (tid > 0) ? g_idata[tid - 1] : 0;
    __syncthreads();
    for (int offset = 1; offset < n; offset *= 2) {
        pout = 1 - pout;                // Swap Double Buffer Indices
        pin = 1 - pout;
        if (tid >= offset) {
            temp[pout * n + tid] += temp[pin * n + tid - offset];
        } else {
            temp[pout * n + tid] = temp[pin * n + tid];
        }
        __syncthreads();
    }
    g_odata[tid] = temp[pout * n + tid]; // Write Output
}
cs

· 입력을 Shared Memory인 "temp[]"에 복사한 후, "offset = 1, 2, 4, …"로 증가시키면서 Prefix Sum 계산
· pin, pout : 2개의 Buffer를 번갈아 사용해서 이전 단계 값을 읽고, 다음 단계 결과를 다른 위치에 저장하는 Double Buffering 사용
    (∵ 같은 배열을 바로 덮어쓰면 다른 Thread가 아직 읽어야 할 값을 변경할 수 있음)
· 현재 단계에서 Offset만큼 왼쪽 값을 읽을 수 있는 Thread는 값을 합하고, 읽을 수 없는 Thread는 현재 값 유지

8 : Exclusive Scan Code로 입력을 한 칸 오른쪽으로 밀고 1번째 값을 0으로 둠
9 : 모든 Thread가 Shared Memory에 초기값을 다 넣을 때까지 대기

** "__syncthreads( )"에 의해 Code가 Architecture의 한계를 초과하지 않는 원소 수에 대해서만 올바르게 동작
∵ "__syncthreads( )"는 동일 Thread Block 내의 Thread들만 동기화 가능
    → 여러 Block에 분산된 데이터는 "__syncthreads( )"로 동기화가 불가능하여 올바른 Prefix Sum 의존성을 보장 불가
∴ Naive Scan 구현은 입력 배열 전체가 하나의 Thread Block 안에 포함될 때만 올바르게 동작

▣ Brent-Kung Algorithm
· 배열 값을 Shared Memory로 불러온 후, Up-Sweep과 Down-Sweep 단계를 거쳐 Scan 수행
· 1번째 For 구문 : Up-Sweep 단계로, 부분합을 만들며 Reduction 하는 단계
· if (thid == 0) {temp[n - 1] = 0;} : 마지막 원소를 0으로 만들어 Exclusive Scan 적용
· 2번째 For 구문 : Down-Sweep 단계로, Prefix 전파 단계
· Shared Memory의 결과를 Global Memory로 다시 저장
→ Shared Memory에 접근하는 방식(temp[bi] += temp[ai];)이 심각한 Bank Conflict 유발

** GPU Shared Memory는 여러 개의 Bank로 나뉘어 있어서 여러 Thread가 동시에 Shared Memory에 접근 가능
→ 여러 Thread가 서로 다른 Bank에 접근하면 병렬로 빠르게 처리
→ 여러 Thread가 서로 같은 Bank에 접근하면 충돌 발생 ▶ 접근이 직렬로 전환

> Without Padding - Bank Conflict O
· Offset = 1 : Memory 주소(ai) 접근 간격이 2가 되면서 2-Way Bank Conflict 발생
· Offset = 2 : Memory 주소(ai) 접근 간격이 4가 되면서 4-Way Bank Conflict 발생

> With Padding - Bank Conflict X
· Shared Memory Index에 추가 Offset을 더하여 Memory Bank 접근을 의도적으로 엇갈리게 만드는 방식 사용
→ Bank Conflict 제거

▣ Large Array
· 하나의 Thread Block 크기를 넘어가서 1 Block으로 처리할 수 없는 큰 배열은 여러 Block으로 나눠서 단계별로 처리
① 전체 배열을 Block 단위로 나눔
② 각 Block 내부에서 Scan 수행
③ 각 Block의 총합을 Auxiliary Array에 저장
④ Auxiliary Array에 대해 다시 Scan 수행
⑤ 각 Block의 결과에 이전 Block들의 누적합을 더함(Ex. Block 2 = Block 0 + Block 1)

[Max Reduction]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#include <cuda_runtime.h>
#include <float.h>
 
__global__ void max_reduce_kernel(float* g_odata, const float* g_idata, int n) {
    extern __shared__ float sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int global_idx = blockIdx.x * blockDim.x + threadIdx.x;
 
    // 1. Global memory에서 Shared memory로 로드
    // 범위를 벗어난 thread는 매우 작은 값으로 채움
    if (global_idx < n) {sdata[tid] = g_idata[global_idx];}
    else {sdata[tid] = -FLT_MAX;}
 
    __syncthreads();
 
    // 2. Block 내부에서 parallel max reduction
    for (unsigned int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
        if (tid < stride) {sdata[tid] = fmaxf(sdata[tid], sdata[tid + stride]);}
        __syncthreads();
    }
 
    // 3. 각 block의 최댓값 저장
    if (tid == 0) {g_odata[blockIdx.x] = sdata[0];}
}
 
float run_max_reduce(float* d_input, int n) {
    int threads = 256;
    int blocks = (n + threads - 1/ threads;
 
    float* d_buf1;
    float* d_buf2;
 
    cudaMalloc(&d_buf1, blocks * sizeof(float));
    cudaMalloc(&d_buf2, blocks * sizeof(float));
 
    float* d_in = d_input;
    float* d_out = d_buf1;
 
    int current_n = n;
 
    while (current_n > 1) {
        blocks = (current_n + threads - 1/ threads;
 
        max_reduce_kernel<<<blocks, threads, threads * sizeof(float)>>>(
            d_out,
            d_in,
            current_n
        );
 
        cudaDeviceSynchronize();
 
        current_n = blocks;
        std::swap(d_in, d_out);
    }
 
    float result;
    cudaMemcpy(&result, d_in, sizeof(float), cudaMemcpyDeviceToHost);
 
    cudaFree(d_buf1);
    cudaFree(d_buf2);
 
    return result;
}
cs

① 각 Thread가 원소 하나를 Shared Memory에 Load
② stride = blockDim.x / 2부터 시작
③ sdata[tid]와 sdata[tid + stride] 중 큰 값 저장
④ stride를 절반으로 줄이며 반복
⑤ sdata[0]에 Block의 Max가 남음

** Sum Reduction처럼 초기 항등원 값을 0으로 설정 시 잘못된 결과값 출력
∴ "- FLT_MAX(- ∞)"로 초기화

** "temp[bi] += temp[ai]" 대신 "temp[bi] = max(temp[bi], temp[ai])" 형태로 비교 연산 수행(float인 경우 "fmaxf" 사용)

[Reference]

· 10_prefix_sum_mgp_2026 (MGP) - Yongjun Park

댓글

이 블로그의 인기 게시물

[Mini-NPU RTL] NN Reference Model

[Mini-NPU RTL] TPU (Study Paper)