"CUDA Sum Reduction"

[Objective]

CUDA Reduction Kernel의 단계별 최적화를 적용하여 대용량 데이터의 합산 연산 성능을 향상

[Implementation]

▣ reduce1

> Kernel Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
__global__ void reduce1(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (i < n) {
        sdata[tid] = g_idata[i];
    } else {
        sdata[tid] = 0;
    }
    __syncthreads();
 
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

2 : Kernel 실행 시 동적으로 Shared Memory 할당(각 Block마다 따로 존재 & Block 내부 Thread들이 함께 사용)

4 : Block 내부 Thread 번호

5 : 전체 입력 배열에서 해당 Thread가 담당할 Global Index

7~11 : 입력 범위를 벗어나지 않는 경우 = Global Memory 값을 Shared Memory로 복사 / 범위를 벗어나는 경우 = 0 할당

→ Reduction을 빠른 Shared Memory에서 수행 & 마지막 Block에서 원소 개수가 부족해도 안전하게 동작

12 : 모든 Thread가 sdata[tid]에 값을 할당할 때까지 대기 (∵ 아직 값을 안 넣었는데 다른 Thread가 해당 값 Read하는 것을 방지)

14~19 : Partial Sum 계산

> s = 1

tid 0 : sdata[0] += sdata[1]

tid 2 : sdata[2] += sdata[3]

tid 4 : sdata[4] += sdata[5]

tid 6 : sdata[6] += sdata[7]

→ sdata = [a₀ + a₁, a₁, a₂ + a₃, a₃, a₄ + a₅, a₅, a₆ + a₇, a₇]

> s = 2

tid 0 : sdata[0] += sdata[2]

tid 4 : sdata[4] += sdata[6]

→ sdata = [a₀ + a₁ + a₂ + a₃, a₁, a₂ + a₃, a₃, a₄ + a₅ + a₆ + a₇, a₅, a₆ + a₇, a₇]

> s = 4

tid 0 : sdata[0] += sdata[4]

→ sdata = [a₀ + a₁ + a₂ + a₃ + a₄ + a₅ + a₆ + a₇, a₁, a₂ + a₃, a₃, a₄ + a₅ + a₆ + a₇, a₅, a₆ + a₇, a₇]

18 : 각 Step의 덧셈이 끝난 후 다음 "s" 단계로 넘어가기 전에 모든 Thread의 계산이 끝날 때까지 대기

20~22 : 각 Block의 Thread 0만 결과를 Global Memory에 저장

> Host Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
void reduce_optimize(const int* const g_idata, int* const g_odata, const int* const d_idata, int* const d_odata, const int n) {
    int threads = 512;
    int* d_in = (int*)d_idata;
    int* d_out = d_odata;
    int* d_partial;
 
    int numElements = n;
    while (numElements > 1) {
        int blocks = (numElements + threads - 1) / threads;
        reduce1<<<blocks, threads, threads * sizeof(int)>>>(d_in, d_out, numElements);
        numElements = blocks;
        d_partial = d_in;
        d_in = d_out;
        d_out = d_partial;
    }
    if (d_in != d_odata) {
        cudaMemcpy(d_odata, d_in, sizeof(int), cudaMemcpyDeviceToDevice);
    }
Colored by Color Scripter
cs

2 : Block 하나당 Thread 수를 512개로 설정

3 : 현재 입력 할당

4 : 현재 출력 할당

5 : "d_in"과 "d_out"을 바꾸기 위한 임시 Pointer 변수

8 : 현재 남아있는 원소 개수가 1개가 될 때까지 반복

9 : 현재 원소 개수 기준으로 필요한 Block 개수 계산

10 : 현재 입력 "d_in"을 "reduce1" Kernel로 Reduction하여 Block별 Partial Sum을 "d_out"에 저장

11 : Kernel 실행 후 결과 원소 개수는 Block 개수와 동일하므로 원소 개수 Update

12~14 : "d_in"과 "d_out" Swap

16~18 : 반복이 끝났을 때 최종 결과가 "d_odata[0]"에 없는 경우 최종 결과가 있는 "d_in[0]"을 "d_odata[0]"로 복사

** 각 Block 안에서만 합을 구해서 g_odata[blockIdx.x]에 Partial Sum을 저장 = 한 단계 Sum Reduction을 수행하는 Kernel

→ Host 함수에서 "reduce1" Kernel을 반복 호출하여 전체 합 계산

▼ "if (tid % (2 * s) == 0)" 때문에 일부 Thread만 실행(Warp Divergence↑) & "%" 연산은 GPU에서 상대적으로 느림

▣ reduce2

> Kernel Code

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
__global__ void reduce2(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (i < n) {
        sdata[tid] = g_idata[i];
    } else {
        sdata[tid] = 0;
    }
    __syncthreads();
 
    for (unsigned int s = 1; s < blockDim.x; s *= 2) {
        unsigned int index = 2 * s * tid;
 
        if (index < blockDim.x) {
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

· "reduce1" Kernel의 단점만 보완하여 "reduce1" Kernel과 대체로 동일

14~21 : 계산에 참여할 Shared Memory 위치를 직접 만들어 GPU에서 느린 "%" 연산 제거

> s = 1

tid 0 = index 0 : sdata[0] += sdata[1]

tid 1 = index 2 : sdata[2] += sdata[3]

tid 2 = index 4 : sdata[4] += sdata[5]

tid 3 = index 6 : sdata[6] += sdata[7]

tid 4 = index 8 : 실행 X

→ sdata = [a₀ + a₁, a₁, a₂ + a₃, a₃, a₄ + a₅, a₅, a₆ + a₇, a₇]

> s = 2

tid 0 = index 0 : sdata[0] += sdata[2]

tid 1 = index 4 : sdata[4] += sdata[6]

tid 2 = index 8 : 실행 X

→ sdata = [a₀ + a₁ + a₂ + a₃, a₁, a₂ + a₃, a₃, a₄ + a₅ + a₆ + a₇, a₅, a₆ + a₇, a₇]

> s = 4

tid 0 = index 0 : sdata[0] += sdata[4]

tid 1 = index 8 : 실행 X

→ sdata = [a₀ + a₁ + a₂ + a₃ + a₄ + a₅ + a₆ + a₇, a₁, a₂ + a₃, a₃, a₄ + a₅ + a₆ + a₇, a₅, a₆ + a₇, a₇]

> Host Code

· "reduce1" Kernel을 실행하는 Host Code와 대부분 동일

· 현재 입력을 "reduce2" Kernel로 Reduction 하도록 변경 (reduce1<<<…>>>(…); ▶ reduce2<<<…>>>(…);)

** "reduce1" Kernel처럼 각 Block 안에서만 합을 구해서 Partial Sum을 저장 = 한 단계 Sum Reduction을 수행하는 Kernel

→ Host 함수에서 "reduce2" Kernel을 반복 호출하여 전체 합 계산

▼ Shared Memory 위치가 (Stride 2 = 0, 2, 4, 6, …), (Stride 3 = 0, 4, 8, …)처럼 간격 증가

→ 여러 Thread가 같은 Memory Bank에 동시에 접근하면서 Shared Memory 접근 속도↓ = Bank Conflict

** reduce 1 : 띄엄띄엄 있는 Thread가 띄엄띄엄 있는 주소에 접근 → Bank Conflict 영향↓

** reduce 2 : 연속된 Thread가 띄엄띄엄 있는 주소에 접근 = 같은 Warp 안에서 Bank 겹치기 쉬움 → Bank Conflict 영향↑

▣ reduce3

> Kernel Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
__global__ void reduce3(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
 
    if (i < n) {
        sdata[tid] = g_idata[i];
    } else {
        sdata[tid] = 0;
    }
    __syncthreads();
 
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

· "reduce2" Kernel의 단점만 보완하여 "reduce2" Kernel과 대체로 동일

· 연속 Thread가 띄엄띄엄 떨어진 주소에 접근한 "reduce2" Kernel과 달리 연속된 Shared Memory 주소 사용

→ Bank Conflict↓

14~19 : "index" 방식을 제거하고, "blockDim.x / 2"부터 절반씩 줄여나가는 Reverse Loop 사용

> s = 4

tid 0 : sdata[0] += sdata[4]

tid 1 : sdata[1] += sdata[5]

tid 2 : sdata[2] += sdata[6]

tid 3 : sdata[3] += sdata[7]

→ sdata = [a₀ + a₄, a₁ + a₅, a₂ + a₆, a₃ + a₇, a₄, a₅, a₆, a₇]

> s = 2

tid 0 : sdata[0] += sdata[2]

tid 1 : sdata[1] += sdata[3]

→ sdata = [a₀ + a₄ + a₂ + a₆, a₁ + a₅ + a₃ + a₇, a₂ + a₆, a₃ + a₇, a₄, a₅, a₆, a₇]

> s = 1

tid 0 : sdata[0] += sdata[1]

→ sdata = [a₀ + a₄ + a₂ + a₆ + a₁ + a₅ + a₃ + a₇, a₁ + a₅ + a₃ + a₇, a₂ + a₆, a₃ + a₇, a₄, a₅, a₆, a₇]

> Host Code

· "reduce2" Kernel을 실행하는 Host Code와 대부분 동일

· 현재 입력을 "reduce3" Kernel로 Reduction 하도록 변경 (reduce2<<<…>>>(…); ▶ reduce3<<<…>>>(…);)

** "reduce2" Kernel처럼 각 Block 안에서만 합을 구해서 Partial Sum을 저장 = 한 단계 Sum Reduction을 수행하는 Kernel

→ Host 함수에서 "reduce3" Kernel을 반복 호출하여 전체 합 계산

▼ 첫 반복이 "s = blockDim.x / 2"이므로 Thread의 절반만 동작, 절반은 IDLE 상태 → GPU 자원 낭비

▣ reduce4

> Kernel Code

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
__global__ void reduce4(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
 
    int sum = 0;
    if (i < n) {
        sum += g_idata[i];
    }
    if (i + blockDim.x < n) {
        sum += g_idata[i + blockDim.x];
    }
    sdata[tid] = sum;
    __syncthreads();
 
    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

· "reduce3" Kernel의 단점만 보완하여 "reduce3" Kernel과 대체로 동일

5 : 한 Block이 "blockDim.x"개의 원소를 처리한 "reduce1, 2, 3"과 달리, 한 Block이 "blockDim.x * 2"개 원소 처리

7 : Thread가 담당할 부분합을 0으로 초기화

8~10 : "i"가 전체 배열 크기 "n"보다 작은 경우 Thread 번호에 해당하는 Index에 있는 원소를 Read

11~13 : "blockDim.x"만큼 떨어진 위치의 원소도 Read하여 "i" Index에 있는 원소와 더함

14 : 각 Thread는 2개의 값을 Read하여 더한 값을 Shared Memory에 할당

→ sdata = [a₀ + a₈, a₁ + a₉, a₂ + a₁₀, a₃ + a₁₁, a₄ + a₁₂, a₅ + a₁₃, a₆ + a₁₄, a₇ + a₁₅]

> s = 4

tid 0 : sdata[0] += sdata[4]

tid 1 : sdata[1] += sdata[5]

tid 2 : sdata[2] += sdata[6]

tid 3 : sdata[3] += sdata[7]

→ sdata = [a₀ + a₈ + a₄ + a₁₂, a₁ + a₉ + a₅ + a₁₃, a₂ + a₁₀ + a₆ + a₁₄, a₃ + a₁₁ + a₇ + a₁₅, a₄ + a₁₂, a₅ + a₁₃, a₆ + a₁₄, a₇ + a₁₅]

> s = 2

tid 0 : sdata[0] += sdata[2]

tid 1 : sdata[1] += sdata[3]

→ sdata = [a₀ + a₈ + a₄ + a₁₂ + a₂ + a₁₀ + a₆ + a₁₄, a₁ + a₉ + a₅ + a₁₃ + a₃ + a₁₁ + a₇ + a₁₅, 생략, 생략, …, 생략]

> s = 1

tid 0 : sdata[0] += sdata[1]

→ sdata = [a₀ + a₈ + a₄ + a₁₂ + a₂ + a₁₀ + a₆ + a₁₄ + a₁ + a₉ + a₅ + a₁₃ + a₃ + a₁₁ + a₇ + a₁₅, 생략, 생략, …, 생략]

> Host Code

· "reduce3" Kernel을 실행하는 Host Code와 대부분 동일

· 한 Block이 처리하는 원소 개수 증가 ((num_E + threads - 1) / threads ▶ (num_E + (threads * 2) - 1) / (threads * 2))

· 현재 입력을 "reduce4" Kernel로 Reduction 하도록 변경 (reduce3<<<…>>>(…); ▶ reduce4<<<…>>>(…);)

** "reduce1, 2, 3" Kernel과 달리 각 Block 안에서 "2 * blockDim.x"개 원소의 합을 구해서 Partial Sum을 저장

= 한 단계 Sum Reduction을 수행하는 Kernel인 것은 동일

→ Host 함수에서 "reduce4" Kernel을 반복 호출하여 전체 합 계산

▼ Sum Reduction은 연산량이 적고 Memory 접근이 중요한 작업

→ Memory 자체가 병목이라기 보다는 명령어 Overhead가 병목일 가능성 존재

** 명령어 Overhead : 실제 연산에 쓰이는 명령어가 아닌, 주소 계산, 반복문 조건 검사, Iteration 증가 같은 부가 명령어

▣ reduce5

> Kernel Code

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
__global__ void reduce5(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
 
    int sum = 0;
    if (i < n) {
        sum += g_idata[i];
    }
    if (i + blockDim.x < n) {
        sum += g_idata[i + blockDim.x];
    }
    sdata[tid] = sum;
    __syncthreads();
 
    for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid < 32) {
        volatile int* smem = sdata;
        smem[tid] += smem[tid + 32];
        smem[tid] += smem[tid + 16];
        smem[tid] += smem[tid + 8];
        smem[tid] += smem[tid + 4];
        smem[tid] += smem[tid + 2];
        smem[tid] += smem[tid + 1];
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

· "reduce4" Kernel과 대체로 동일하지만, 마지막 32개 Thread가 남은 구간의 반복문을 Unrolling

17 : "s > 0"이 아닌 "s > 32" 까지만 반복

23~31 : "s ≤ 32"부터는 하나의 Warp 안에서 실행되므로 동기화 없이 Sum Reduction 단계를 Unrolling

24 : 동기화를 하지 않으므로 "volatile"을 사용하여 항상 Shared Memory의 최신 값을 반드시 읽도록 함

(※ volatile : Compiler가 Shared Memory 값을 Register에만 저장하여 사용하는 최적화를 방지)

** 1 Warp = 32 Threads → 반복문 Overhead(Iteration 증가, 조건 확인 등) & 불필요한 동기화 제거

(동기화 없어도 되는 이유 : 하나의 Warp 안에서는 모든 Thread가 동일한 명령어를 동시에 실행)

> Host Code

· "reduce4" Kernel을 실행하는 Host Code와 대부분 동일

· 현재 입력을 "reduce5" Kernel로 Reduction 하도록 변경 (reduce4<<<…>>>(…); ▶ reduce5<<<…>>>(…);)

** "reduce4" Kernel처럼 각 Block 안에서 "2 * blockDim.x"개 원소의 합을 구해서 Partial Sum을 저장

= 한 단계 Sum Reduction을 수행하는 Kernel

→ Host 함수에서 "reduce5" Kernel을 반복 호출하여 전체 합 계산

▣ reduce6

> Kernel Code

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
template <unsigned int blockSize>
__global__ void reduce6(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
 
    int sum = 0;
    if (i < n) {
        sum += g_idata[i];
    }
    if (i + blockDim.x < n) {
        sum += g_idata[i + blockDim.x];
    }
    sdata[tid] = sum;
    __syncthreads();
 
    if (blockSize >= 512) {
        if (tid < 256) {
            sdata[tid] += sdata[tid + 256];
        }   __syncthreads();
    }
    if (blockSize >= 256) {
        if (tid < 128) {
            sdata[tid] += sdata[tid + 128];
        }   __syncthreads();
    }
    if (blockSize >= 128) {
        if (tid < 64) {
            sdata[tid] += sdata[tid + 64];
        }   __syncthreads();
    }
    if (tid < 32) {
        volatile int* smem = sdata;
        if (blockSize >= 64) {
            smem[tid] += smem[tid + 32];
        }
        if (blockSize >= 32) {
            smem[tid] += smem[tid + 16];
        }
        if (blockSize >= 16) {
            smem[tid] += smem[tid + 8];
        }
        if (blockSize >= 8) {
            smem[tid] += smem[tid + 4];
        }
        if (blockSize >= 4) {
            smem[tid] += smem[tid + 2];
        }
        if (blockSize >= 2) {
            smem[tid] += smem[tid + 1];
        }
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

· "reduce5" Kernel과 대체로 동일하지만, 마지막 Warp뿐만 아니라 Block 내부 Reduction 전체를 Unrolling

1 : "blockSize"를 Template Parameter로 Load → Compile Time에 결정

18~32 : 앞쪽 Thread가 뒤쪽 Thread를 더하는 반복문을 Unrolling(Runtime이 아닌 Compile Time에 조건 검사)

> Host Code

· "reduce5" Kernel을 실행하는 Host Code와 대부분 동일

· 현재 입력을 "reduce6" Kernel로 Reduction 하도록 변경 (reduce5<<<…>>>(…); ▶ reduce6<512><<<…>>>(…);)

→ <512> : blockSize = 512 버전으로 실행 ("threads" 값과 일치해야 올바르게 동작)

** "reduce5" Kernel처럼 각 Block 안에서 "2 * blockDim.x"개 원소의 합을 구해서 Partial Sum을 저장

= 한 단계 Sum Reduction을 수행하는 Kernel

→ Host 함수에서 "reduce6" Kernel을 반복 호출하여 전체 합 계산

▣ reduce7

> Kernel Code

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
template <unsigned int blockSize>
__global__ void reduce7(const int* g_idata, int* g_odata, int n) {
    extern __shared__ int sdata[];
 
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockSize * 2) + threadIdx.x;
    unsigned int gridSize = blockSize * 2 * gridDim.x;
    sdata[tid] = 0;
    while (i < n) {
        sdata[tid] += g_idata[i];
 
        if (i + blockSize < n) {
            sdata[tid] += g_idata[i + blockSize];
        }
        i += gridSize;
    }
    __syncthreads();
 
    if (blockSize >= 512) {
        if (tid < 256) {
            sdata[tid] += sdata[tid + 256];
        }   __syncthreads();
    }
    if (blockSize >= 256) {
        if (tid < 128) {
            sdata[tid] += sdata[tid + 128];
        }   __syncthreads();
    }
    if (blockSize >= 128) {
        if (tid < 64) {
            sdata[tid] += sdata[tid + 64];
        }   __syncthreads();
    }
    if (tid < 32) {
        volatile int* smem = sdata;
        if (blockSize >= 64) {
            smem[tid] += smem[tid + 32];
        }
        if (blockSize >= 32) {
            smem[tid] += smem[tid + 16];
        }
        if (blockSize >= 16) {
            smem[tid] += smem[tid + 8];
        }
        if (blockSize >= 8) {
            smem[tid] += smem[tid + 4];
        }
        if (blockSize >= 4) {
            smem[tid] += smem[tid + 2];
        }
        if (blockSize >= 2) {
            smem[tid] += smem[tid + 1];
        }
    }
    if (tid == 0) {
        g_odata[blockIdx.x] = sdata[0];
    }
}
Colored by Color Scripter
cs

· "reduce6" Kernel과 대체로 동일하지만, 각 Thread가 처음에 여러 개의 원소를 반복해서 더하는 순차 처리 작업이 증가

→ 각 Thread마다 작업량 증가

7 : Grid 전체가 한 번에 처리하는 원소 수

9~16 : 각 Thread가 2개의 원소만 Read 하지 않고, Grid-Stride 방식으로 여러 원소를 반복해서 더함

> Host Code

1
2
3
4
5
6
7
8
9
10
11
void reduce_optimize(const int* const g_idata, int* const g_odata, const int* const d_idata, int* const d_odata, const int n) {
    int threads = 512;
    int blocks = 128;
    int* d_partial = 0;
    cudaMalloc(&d_partial, blocks * sizeof(int));
 
    reduce7<512><<<blocks, threads, threads * sizeof(int)>>>(d_idata, d_partial, n);
    reduce7<512><<<1, threads, threads * sizeof(int)>>>(d_partial, d_odata, blocks);
 
    cudaFree(d_partial);
}
Colored by Color Scripter
cs