CUDA 기초: 스레드 구조, 메모리 복사, 벡터 덧셈과 SGEMM

CUDA를 처음 공부하면서 확인한 Host/Device의 관계, 커널 실행 구조, 메모리 복사, 벡터 덧셈과 기본적인 행렬 곱셈을 정리한다.

Host와 Device

CUDA 프로그램에서 CPU와 CPU 메모리는 Host, GPU와 GPU 메모리는 Device라고 부른다.

일반적인 CUDA 프로그램은 다음 순서로 동작한다.

Host 메모리를 할당하고 입력 데이터를 준비한다.
Device 메모리를 할당한다.
입력 데이터를 Host에서 Device로 복사한다.
GPU 커널을 실행한다.
결과를 Device에서 Host로 복사한다.
Host와 Device 메모리를 각각 해제한다.

int *a = static_cast<int *>(malloc(size));
int *d_a;

cudaMalloc(reinterpret_cast<void **>(&d_a), size);
cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice);

// GPU kernel launch

cudaMemcpy(a, d_a, size, cudaMemcpyDeviceToHost);

cudaFree(d_a);
free(a);

a와 d_a는 연결된 하나의 배열이 아니다. a는 시스템 RAM을, d_a는 GPU 메모리를 가리키며 두 공간은 서로 독립적이다. cudaMemcpy()를 호출할 때만 데이터가 복사된다.

`cudaMalloc()`과 이중 포인터

전통적인 cudaMalloc() 선언은 다음과 같다.

cudaError_t cudaMalloc(void **devPtr, size_t size);

int *d_a;
cudaMalloc(reinterpret_cast<void **>(&d_a), size);

여기에서 각 표현의 타입은 다음과 같다.

d_a                                  int *
&d_a                                 int **
reinterpret_cast<void **>(&d_a)      void **

cudaMalloc()은 할당한 Device 주소를 포인터 변수 d_a에 기록해야 한다. 따라서 d_a의 값이 아니라 d_a 자체의 주소인 &d_a를 전달한다. reinterpret_cast는 주소의 비트 값을 바꾸는 것이 아니라 컴파일러가 그 주소를 다른 포인터 타입으로 해석하도록 한다.

CUDA C++에서는 환경에 따라 다음과 같이 간단히 쓸 수도 있다.

cudaMalloc(&d_a, size);

`cudaMemcpy()`와 DMA

cudaMemcpy(d_a, a, size, cudaMemcpyHostToDevice);

인자의 의미는 다음과 같다.

d_a                       목적지 Device 주소
a                         원본 Host 주소
size                      복사할 바이트 수
cudaMemcpyHostToDevice    복사 방향

cudaMemcpyHostToDevice는 임의의 주소를 GPU 주소로 변환하지 않는다. 목적지 포인터는 cudaMalloc() 등으로 얻은 유효한 Device 주소여야 하고, 원본 포인터 역시 유효한 Host 주소여야 한다.

실제 대량 전송은 일반적으로 GPU의 copy engine과 DMA를 이용하며 PCIe 또는 NVLink를 통해 이루어진다. malloc()으로 할당한 Host 메모리는 pageable memory이므로 드라이버 내부의 pinned buffer를 한 번 거칠 수 있다.

pageable Host memory
        ↓ CPU 복사
temporary pinned memory
        ↓ DMA / PCIe
Device memory

Host 메모리를 cudaMallocHost()로 할당하면 page-locked(pinned) memory가 되어 비동기 전송과 더 효율적인 DMA 전송에 이용할 수 있다.

float *a;
cudaMallocHost(&a, size);

// 사용 후
cudaFreeHost(a);

Pinned memory는 시스템의 메모리 관리에 부담을 주므로 필요한 만큼만 사용하는 것이 좋다.

CUDA 커널과 실행 구성

__global__로 선언한 함수는 Host에서 호출하고 Device에서 실행하는 CUDA 커널이다.

__global__ void kernel()
{
    // GPU 스레드 하나가 실행할 코드
}

커널 호출의 기본 형태는 다음과 같다.

kernel<<<grid의 block 수, block당 thread 수>>>();

예를 들어 다음 호출은 블록 하나에 N개의 스레드를 만든다.

device_add<<<1, N>>>(d_a, d_b, d_c);

반대로 다음 호출은 N개의 블록마다 스레드를 하나씩 만든다.

device_add<<<N, 1>>>(d_a, d_b, d_c);

두 경우 모두 전체 스레드 수는 N이지만, GPU는 스레드를 보통 32개 단위인 warp로 실행하므로 블록당 스레드가 하나인 구성은 비효율적이다.

CUDA의 3차원 인덱스

CUDA는 Grid와 Block을 최대 3차원으로 구성할 수 있다.

threadIdx.x;
threadIdx.y;
threadIdx.z;

blockIdx.x;
blockIdx.y;
blockIdx.z;

threadIdx.x의 x는 블록 내부에서 현재 스레드의 x축 좌표다. 1차원 블록에서는 사실상 블록 내부의 스레드 번호다.

여러 블록을 사용하는 1차원 작업에서는 전체 데이터에 대응하는 global index를 다음과 같이 계산한다.

int idx = blockIdx.x * blockDim.x + threadIdx.x;

일반적인 벡터 덧셈 커널은 다음과 같다.

__global__ void vector_add(const int *a, const int *b, int *c, int n)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;

    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}

호출할 때는 필요한 블록 수를 올림 나눗셈으로 구한다.

int threads_per_block = 256;
int blocks = (N + threads_per_block - 1) / threads_per_block;

vector_add<<<blocks, threads_per_block>>>(d_a, d_b, d_c, N);

스레드 수가 N보다 많아질 수 있기 때문에 커널의 if (idx < n) 경계 검사가 필요하다.

커널에는 무엇을 작성하는가

커널에는 기본적으로 GPU 스레드 하나가 수행할 작업을 작성한다. 모든 스레드가 같은 커널 코드를 실행하지만 threadIdx와 blockIdx가 다르기 때문에 서로 다른 데이터를 처리한다.

스레드 하나가 반드시 원소 하나만 처리해야 하는 것은 아니다. 하나의 스레드가 여러 원소를 처리할 수 있고, 같은 블록의 스레드들은 shared memory와 __syncthreads()를 이용해 협력할 수도 있다.

기본 SGEMM

SGEMM은 single-precision general matrix multiplication이며 다음 연산을 뜻한다.

\[C = \alpha AB + \beta C\]

행렬 크기는 다음과 같다.

A: N × K
B: K × M
C: N × M

가장 단순한 GPU 구현에서는 스레드 하나가 결과 행렬 C의 원소 하나를 계산한다.

__global__ void sgemm_gpu_kernel(
    const float *A,
    const float *B,
    float *C,
    int N,
    int M,
    int K,
    float alpha,
    float beta)
{
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;

    if (row >= N || col >= M) {
        return;
    }

    float sum = 0.0f;
    for (int i = 0; i < K; ++i) {
        sum += A[row * K + i] * B[i * M + col];
    }

    C[row * M + col] = alpha * sum + beta * C[row * M + col];
}

row와 col은 이 스레드가 담당하는 C[row][col]을 나타낸다. 반복문에서는 A의 row번째 행과 B의 col번째 열을 내적한다.

행렬은 선형 메모리에 row-major 순서로 저장되므로, 열 개수가 width인 행렬의 2차원 좌표는 다음과 같이 변환한다.

matrix[row * width + col]

따라서 각 행렬의 올바른 접근 방식은 다음과 같다.

A[row * K + i]   // A의 열 개수: K
B[i * M + col]   // B의 열 개수: M
C[row * M + col] // C의 열 개수: M

2차원 실행 구성은 다음과 같이 만들 수 있다.

dim3 dimBlock(16, 16);
dim3 dimGrid(
    (M + dimBlock.x - 1) / dimBlock.x,
    (N + dimBlock.y - 1) / dimBlock.y
);

sgemm_gpu_kernel<<<dimGrid, dimBlock>>>(
    A, B, C, N, M, K, alpha, beta
);

이 방식은 구조를 이해하기 쉽지만 A와 B를 global memory에서 반복해서 읽기 때문에 고성능 구현은 아니다. 실제 최적화에서는 행렬 일부를 shared memory 타일에 적재하여 블록의 스레드들이 데이터를 재사용한다.

CUDA Event로 GPU 시간 측정

GPU 커널 호출은 비동기이므로 일반 CPU 타이머로 호출 전후만 측정하면 GPU의 실제 실행 시간을 제대로 얻지 못할 수 있다. CUDA Event를 같은 CUDA stream에 기록하면 GPU 작업 시간을 측정할 수 있다.

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);

cudaEventRecord(start);

for (int i = 0; i < test_iterations; ++i) {
    sgemm(A, B, C, N, M, K, alpha, beta);
}

cudaEventRecord(stop);
cudaEventSynchronize(stop);

float elapsed_ms = 0.0f;
cudaEventElapsedTime(&elapsed_ms, start, stop);

printf("Average Time = %.4f ms\n", elapsed_ms / test_iterations);

cudaEventDestroy(start);
cudaEventDestroy(stop);

cudaEventSynchronize(stop)은 stop event 이전에 제출된 GPU 작업이 끝날 때까지 Host를 기다리게 한다.

정리

CUDA 커널에는 각 GPU 스레드가 수행할 작업을 작성한다.
<<<grid, block>>>으로 블록 수와 블록당 스레드 수를 결정한다.
전체 데이터 인덱스는 일반적으로 blockIdx * blockDim + threadIdx로 계산한다.
Host 메모리와 Device 메모리는 독립적이며 cudaMemcpy()로 데이터를 이동한다.
Host와 Device 사이의 대량 전송에는 일반적으로 DMA와 GPU copy engine이 사용된다.
기본 SGEMM에서는 스레드 하나가 결과 행렬의 원소 하나를 계산한다.
GPU 실행 시간은 비동기 실행을 고려해 CUDA Event로 측정할 수 있다.