Parallel Prefix Sum (2)

References

• Programming Massively Parallel Processors

Contents

• Brent-Kung adder Algorithm

---

A More Work-Efficient Parallel Scan

Parallel Prefix Sum (1)

이전 포스팅에서 살펴본 Kogge-Stone 커널은 단순하고, 실제 어플리케이션에서의 효율성이 상당히 낮습니다.

위에서 살펴봤듯이, 어떠한 값들의 집합의 합을 계산하는 가장 빠른 병렬 방법은 reduction tree입니다. 충분한 execution units이 있다면, reduction tree는 \(log_2 N\)의 time unit으로 N개의 값에 대한 합을 계산할 수 있습니다.

 

Tree는 출력값 계산에 사용할 수 있는 여러 개의 sub-sum을 생성할 수 있는데, 이는 Brent-Kung adder 디자인의 기초가 되는 알고리즘입니다.

아래 그림에서 단지 4단계만에 16개 원소의 합을 구합니다.

A parallel inclusive scan algorithm based on the Brent-Kung adder design.

알고리즘의 첫 번째 파트은 그림 상단에 표시된 트리를 사용합니다.

첫 번째 단계에서 XY[i]의 홀수 요소만 XY[i-1] + XY[i]로 업데이트됩니다. 그리고 두 번째 단계에서 XY[4n-1]만 업데이트되는데, 이때 인덱스 3, 7, 11, 15가 이에 해당합니다. 세 번째 단계에서는 XY[8n-1] 원소들만 업데이트되고, 위 그림에서 7과 15가 이에 해당합니다. 그리고 마지막으로 XY[15]가 업데이트됩니다.

이 과정에서 수행된 총 연산의 수는 8+4+2+1=15 입니다. 일반적으로 N개의 요소의 scan section에 대해서 N-1번의 연산을 수행합니다.

알고리즘의 두 번째 파트는 위 그림 하단에 표시된 역 트리를 사용합니다.

 

reduction phase가 끝나면 몇 가지 쓸만한 부분합이 계산됩니다.

아래 표의 첫 번째 행은 위 그림의 위쪽에 있는 reduction tree가 수행된 직후의 부분합을 보여줍니다.

첫 번째 행에서 XY[0], XY[7], XY[15]는 최종 결과를 포함하고 있는 것을 확인할 수 있습니다. 그러므로, 남은 XY 원소들은 4칸 이내에 위치하는 값들로부터 부분합을 계산할 수 있습니다.

Partial sums available in each XY element after the reduction tree phase.

예를 들어, XY[14]를 살펴보겠습니다. XY[14]가 최종값을 얻으려면 XY[7](\(x_0...x_7\)), XY[11](\(x_8...x_{11}\)), XY[13](\(x_{12}...x_{13}\))이 필요합니다. 위쪽에 위치한 reduction tree를 수행한 뒤에 아래쪽의 reduction tree는 나머지 절반의 덧셈 연산을 수행하는데, 첫 번째 reduction phase에서 XY[7]을 4칸 떨어진 위치에 있는 XY[11]에 더합니다. 그 결과 XY[11] = XY[7] + XY[11] 이므로, 이 phase에서 XY[11]은 최종값을 가지게 됩니다. 이는 위의 표의 두 번째 행에 해당합니다.

 

그 다음 단계에서는 2칸 떨어진 위치의 합계를 사용하여 덧셈 연산을 수행합니다. XY[5]는 두 칸 떨어진 XY[3]으로부터 부분합을 계산하고, XY[9]는 XY[7]로부터, XY[13]은 XY[11]로부터 부분합을 계산합니다. XY[5], XY[9], XY[13]은 이 단계에서 최종값을 가지게 됩니다.

그리고 마지막 단계에서 한 칸 떨어진 위치로부터 남은 원소들의 부분합을 계산합니다. XY[2]는 XY[1]로부터, XY[4]는 XY[3]으로부터, XY[6]은 XY[5], XY[8]은 XY[4], .. 이런식으로 남은 XY[2], XY[4], XY[6], XY[8], XY[10], XY[12], XY[14]의 부분합을 계산하게 됩니다.

 

 

알고리즘을 자세히 살펴보면, 사용하는 스레드의 최대 갯수가 전체 16개의 원소의 절반입니다. 이는 reduction tree의 첫 번째 단계를 살펴보면 확인 가능한데, 첫 번째 단계에서 홀수의 인덱스의 XY 원소들만 업데이트된다는 것을 보면 알 수 있습니다. 따라서, CUDA 구현에서 블록의 크기는 SECTION_SIZE의 반이 된다는 것에 생각하면서 커널 사이즈를 설정해야하고, 커널 함수 내에서 blockDim.x의 값을 사용할 때 유의해야 합니다.

(SECTION_SIZE가 256이라면, 블록 내 스레드의 수는 128이고, 이 블록에서 256개의 원소를 커버합니다.)

Brent-Kung Scan 알고리즘을 구현하면 다음과 같습니다.

__global__
void brentKungScan(float* X, float* Y, int n)
{
    __shared__ float XY[SECTION_SIZE];
    int i = 2*blockIdx.x*blockDim.x + threadIdx.x;
    if (i < n)
        XY[threadIdx.x] = X[i];
    if (i + blockDim.x < n)
        XY[threadIdx.x + blockDim.x] = X[i + blockDim.x];

    for (unsigned int stride = 1; stride <= blockDim.x; stride *= 2) {
        __syncthreads();
        int index = ((threadIdx.x + 1) * stride * 2) - 1;
        if (index < SECTION_SIZE) {
            XY[index] += XY[index - stride];
        }
    }

    for (unsigned int stride = SECTION_SIZE/4; stride > 0; stride /= 2) {
        __syncthreads();
        int index = ((threadIdx.x + 1) * stride * 2) - 1;
        if (index + stride < SECTION_SIZE) {
            XY[index + stride] += XY[index];
        }
    }

    __syncthreads();
    if (i < n)
        Y[i] = XY[threadIdx.x];
    if (i + blockDim.x < n)
        Y[i + blockDim.x] = XY[threadIdx.x + blockDim.x];
}

line 4-9는 XY 배열의 초기값을 각 섹션에서의 입력 배열의 값으로 설정하는 부분입니다. 블록의 사이즈가 SECTION_SIZE의 절반이기 때문에 두 번의 if문에 걸쳐서 초기화를 진행해주고 있습니다.

 

다음 line 11-17은 처음에 봤던 그림에서 위쪽의 reduction tree의 각 단계를 반복하는 부분입니다.

첫 번째 반복에서 두 칸씩 떨어진 홀수 인덱스인 1, 3, 5, 7, 9, ... 의 XY 원소들의 값에 바로 왼쪽에 위치한 원소의 값을 더해줍니다.

두 번째 반복에서는 네 칸씩 떨어진 인덱스인 3, 7, 11, 15, ...의 XY 원소들의 값을 업데이트합니다. 그 다음은 7, 15, ..., 이런식으로 진행됩니다.

 

line 19-25는 reverse reduction tree의 각 단계를 반복하는 부분입니다.

SECTION_SIZE가 16이라면, 첫 번째 반복에서는 인덱스 7만 활성화되고, stride=4만큼 떨어진 인덱스 11, XY[11]에 XY[7]의 부분합을 더합니다.

두 번째 반복에서는 인덱스 3, 7, 11 이 stride=2만큼 떨어진 XY[5]에 XY[3]을, XY[9]에 XY[7]을, XY[13]에 XY[11]의 부분합을 더합니다. 

그리고 마지막 단계에서 1,3,5,7,9,11,13 인덱스의 XY 값을 2,4,6,8,10,12,14 인덱스의 XY 값에 더해주고 반복은 끝납니다.

 

 

각 단계에서 연산의 횟수를 분석해보겠습니다.

N개의 원소를 입력으로 사용한다고 했을 때, 첫 단계의 reduction tree에서의 반복 횟수는 

\[1 + 2 + 4 + \cdots + \frac{N}{4} + \frac{N}{2} = N-1\]

입니다. 

그리고 두 번째 reverse reduction tree에서의 횟수는

\[(2-1) + (4-1) + \cdots + (\frac{N}{4}-1) + (\frac{N}{2}-1) = N-1-log_2 N\]

입니다. 

따라서 총 횟수는 \(2N-2-log_2 N\)이 되고, 연산 횟수는 N에 비례한다는 것을 볼 수 있습니다.

 

이전 포스팅에서 본 Kogge-Stone 알고리즘보다 이론적으로 더 좋은 작업 효율성을 보여줍니다.

특히 CUDA 커널 구현에서 블록 당 N/2 개의 스레드를 사용하고, 활성화되는 스레드 수가 Kogge-Stone 커널보다 적습니다. 따라서 reduction tree를 통과하는 속도가 Kogge-Stone 커널보다 빠릅니다. 다만, 비활성화되는 스레드도 여전히 실행 리소스에 포함되므로, Brent-Kung 커널이 소비하는 실제 리소스의 양은 \(\frac{N}{2}*(2log_2(N) - 1)\)에 가깝습니다. 만약 1024개의 입력 원소를 32개의 execution unit으로 실행한다면, Brent-Kung 커널은 약 \(512*(2*10 - 1) = 304\)의 time units이 필요하고, 대략 3.4배의 속도 향상을 이끌어낼 수 있습니다.

 

1,000,000개의 입력으로 커널을 수행하면 다음의 결과를 얻을 수 있습니다. (--kernel=1이 Brent-Kung 알고리즘을 선택한다는 의미입니다.)

 

마찬가지로 RTX3080으로 커널을 실행해봤습니다.

 

이전 포스팅에서 블록당 스레드를 512개로 설정했는데, 똑같이 512개를 사용하도록 설정하고 커널을 실행해보면 다음의 결과를 얻을 수 있습니다.

 

 

전체 코드는 아래 링크를 참조하시길 바랍니다.

https://github.com/junstar92/parallel_programming_study/blob/master/CUDA/prefixSum/prefixSum.cu.

GitHub - junstar92/parallel_programming_study: Study parallel programming - CUDA, OpenMP, MPI, Pthread