junyeokk
Blog
Parallel Computing·2025. 10. 23

행렬 덧셈 병렬화

행렬 기초 개념

Matrix의 기본 개념

Matrix는 2D array로 표현된다. 수학적으로는 m×nm \times n 행렬이라고 쓰며, mm은 행의 개수, nn은 열의 개수를 나타낸다. 프로그래밍에서는 nrow와 ncol 변수를 사용하여 이를 표현한다.

각 원소의 위치는 (row, col) 또는 (y, x)로 접근하며, 수학에서는 보통 (i, j)를 사용한다. ii는 행 인덱스, jj는 열 인덱스다. 예를 들어 4행 3열 행렬은 nrow = 4, ncol = 3으로 표현되며, 원소는 M0,0M_{0,0}부터 M3,2M_{3,2}까지 존재한다. 인덱스는 0부터 시작하므로 row #0부터 row #3까지, col #0부터 col #2까지의 원소가 있다.

좌표계를 생각하면 x축은 가로 방향(열), y축은 세로 방향(행)이다. 하지만 주의할 점이 있다. 이미지나 행렬에서 y축은 보통 아래가 아니라 위로 증가한다. (0, 0)이 왼쪽 위 모서리고, (nrow-1, ncol-1)이 오른쪽 아래 모서리가 된다.

Matrix Addition 연산

행렬 덧셈은 같은 크기의 두 행렬을 더하는 연산이다. C=A+BC = A + B로 표현하며, 행렬 AA와 BB는 반드시 같은 크기를 가져야 한다. 둘 다 m×nm \times n이어야 계산이 가능하다. 연산 방식은 간단하다. CC의 각 원소 cijc_{ij}는 aij+bija_{ij} + b_{ij}로 계산된다. 즉, 같은 위치의 원소끼리 더하는 것이다.

수학적으로 표현하면 다음과 같다.

C[i,j]=A[i,j]+B[i,j]for all i,jC[i, j] = A[i, j] + B[i, j] \quad \text{for all } i, j

ii는 0부터 nrow-1까지, jj는 0부터 ncol-1까지 반복한다. 4×44 \times 4 행렬을 예로 들면 총 16번의 덧셈이 필요하며, 일반적으로 m×nm \times n 행렬의 경우 m×nm \times n번의 덧셈을 수행한다.

병렬화 가능성을 생각해보자. 각 원소의 계산은 완전히 독립적이다. c00c_{00}을 계산할 때 c01c_{01}이나 다른 원소가 필요 없다. 모든 원소를 동시에 계산할 수 있다. 이것이 바로 GPU에 적합한 이유다. 1차원 벡터 덧셈에서는 1차원 스레드 인덱스만 있으면 됐지만, 행렬 덧셈에서는 행과 열, 두 개의 인덱스가 필요하다. 다행히 CUDA는 2차원(또는 3차원) 스레드 블록을 지원한다.

Row-major vs Column-major 메모리 레이아웃

메모리는 1차원이다. RAM은 일렬로 늘어선 바이트 배열이다. 2차원 행렬을 1차원 메모리에 저장하려면 변환이 필요하다. 두 가지 방법이 있다.

Row-major 방식은 행을 먼저 저장한다. 첫 번째 행의 모든 원소를 연속으로 저장하고, 그 다음 두 번째 행을 저장한다. C, C++, Python(NumPy)이 이 방식을 사용한다. 4×34 \times 3 행렬을 row-major로 저장하면 M0,0M_{0,0} M0,1M_{0,1} M0,2M_{0,2} M1,0M_{1,0} M1,1M_{1,1} M1,2M_{1,2} M2,0M_{2,0} M2,1M_{2,1} M2,2M_{2,2} M3,0M_{3,0} M3,1M_{3,1} M3,2M_{3,2} 순서가 된다. 첫 번째 행, 두 번째 행, 세 번째 행, 네 번째 행 순서로 배치된다.

반면 Column-major 방식은 열을 먼저 저장한다. 첫 번째 열의 모든 원소를 연속으로 저장하고, 그 다음 두 번째 열을 저장한다. Fortran, MATLAB, R이 이 방식을 사용한다. 같은 4×34 \times 3 행렬을 column-major로 저장하면 M0,0M_{0,0} M1,0M_{1,0} M2,0M_{2,0} M3,0M_{3,0} M0,1M_{0,1} M1,1M_{1,1} M2,1M_{2,1} M3,1M_{3,1} M0,2M_{0,2} M1,2M_{1,2} M2,2M_{2,2} M3,2M_{3,2} 순서가 된다.

C/C++에서 Row-major Layout

C/C++에서 2D array는 row-major로 저장된다. 선언 방식은 정적 할당과 동적 할당 두 가지가 있다. 정적 할당은 float M[4][3] 형태로 컴파일 타임에 크기가 고정되며 스택에 할당된다. 크기가 크면 스택 오버플로우가 발생할 수 있다.

동적 할당은 float* M = new float[4 * 3] 형태로 런타임에 크기를 결정할 수 있다. 힙에 할당되므로 큰 행렬도 문제없다. CUDA에서는 항상 이 방식을 사용한다.

2D 인덱스에서 1D 인덱스로 변환하는 공식은 다음과 같다.

index=row×ncol+col\text{index} = \text{row} \times \text{ncol} + \text{col}

이것이 row-major 레이아웃의 핵심 공식이다. 왜 이 공식이 맞는가? 각 행은 ncol개 원소를 가진다. row 번째 행은 row * ncol 위치에서 시작한다. 그 행에서 col 번째 원소는 row * ncol + col 위치에 있다.

구체적인 예를 보자. 4×34 \times 3 행렬에서 M[2][1] 원소를 찾는다. row = 2, col = 1이므로 index = 2 * 3 + 1 = 7이다. 메모리 배열의 7번째 위치(0부터 시작하므로 8번째 원소)에 있다. C/C++ 2D array 표기법 M[row][col]은 사실 포인터 연산이며, 결과적으로 row * ncol + col 계산과 동일하다.

CUDA에서는 1D 배열로 할당하므로 직접 인덱스 변환을 해야 한다. cudaMalloc(&M, nrow * ncol * sizeof(float))로 할당하고, M[row * ncol + col] 형태로 접근한다. 이것이 모든 CUDA 행렬 연산의 기본 패턴이다.

CPU 구현

CPU 버전 구현

10,000 × 10,000 행렬을 CPU에서 double for loop로 처리하는 구현을 살펴본다.

c
unsigned nrow = 10000;
unsigned ncol = 10000;

int main(const int argc, const char* argv[]) {
    float* matA = new float[nrow * ncol];
    float* matB = new float[nrow * ncol];
    float* matC = new float[nrow * ncol];

    srand(0);
    setNormalizedRandomData(matA, nrow * ncol);
    setNormalizedRandomData(matB, nrow * ncol);

    ELAPSED_TIME_BEGIN(0);
    for (register unsigned r = 0; r < nrow; ++r) {
        for (register unsigned c = 0; c < ncol; ++c) {
            unsigned i = r * ncol + c;
            matC[i] = matA[i] + matB[i];
        }
    }
    ELAPSED_TIME_END(0);

    printMat("matC", matC, nrow, ncol);

    delete[] matA;
    delete[] matB;
    delete[] matC;
}

위 코드의 동작을 살펴보면 먼저 10,000×10,00010,000 \times 10,000 행렬 세 개를 할당한다. 총 100,000,000개 원소다. 각 행렬은 400MB 메모리를 차지한다. 외부 루프는 행을 순회하고 내부 루프는 열을 순회한다. index = row * ncol + col 공식으로 1D 인덱스를 계산한 뒤 matC[i] = matA[i] + matB[i] 덧셈을 수행한다.

register 키워드는 컴파일러에게 루프 변수를 가능하면 레지스터에 유지하라는 힌트를 준다. 현대 컴파일러는 알아서 최적화하므로 큰 효과는 없다. Intel Core i5-3570에서 168,857 usec가 소요된다. 약 169ms다.

printMat() 헬퍼 함수

행렬 출력을 위한 헬퍼 함수가 ./common.cpp에 구현되어 있다. void printMat(const char* name, float* mat, int nrow, int ncol) 형태로 nrow by ncol 크기의 행렬을 출력한다. 하지만 전체를 다 출력하지 않는다. 10,000×10,00010,000 \times 10,000 행렬을 모두 출력하면 화면이 넘친다.

대신 처음과 마지막 3개 행과 열만 출력한다. 디버그 용도로 충분하다. 경계 케이스에서 C[i] = A[i] + B[i] 계산이 맞는지 확인한다. 첫 번째와 마지막 원소가 제대로 계산되었는지 보는 것이 중요하다. 인덱스 계산 실수는 보통 경계에서 드러난다.

CUDA 병렬화 구현

Thread Block 설계와 Grid 구성

행렬 크기가 10,000×10,00010,000 \times 10,000이다. 2차원 데이터이므로 2차원 스레드 블록을 사용한다. 블록 크기로 32×3232 \times 32를 선택한다. 왜 32인가? Warp 크기가 32이기 때문이다. 32의 배수로 블록을 구성하면 warp 경계와 일치하여 비효율이 없다. 32×32=1,02432 \times 32 = 1,024개 스레드는 CUDA가 허용하는 블록당 최대 스레드 개수다. 정확히 최대치를 사용한다.

그렇다면 Grid 크기는 어떻게 계산하는가? 전체 행렬 크기를 블록 크기로 나눈다. 하지만 나머지가 있을 수 있으므로 올림 연산이 필요하다.

c
gridDim.x = (ncol + blockDim.x - 1) / blockDim.x;
gridDim.y = (nrow + blockDim.y - 1) / blockDim.y;

이것은 올림 나눗셈(ceiling division) 공식이다. 정수 나눗셈에서 올림을 구현한다. 구체적으로 계산하면 다음과 같다.

gridDim.x=10000+32−132=1003132=313\text{gridDim.x} = \frac{10000 + 32 - 1}{32} = \frac{10031}{32} = 313

gridDim.y도 동일하게 313이다. Grid는 313×313313 \times 313 블록을 가진다.

총 스레드 개수는 얼마나 될까? 313×32×313×32=10,016×10,016=100,320,256313 \times 32 \times 313 \times 32 = 10,016 \times 10,016 = 100,320,256개다. 실제 필요한 것은 10,000×10,000=100,000,00010,000 \times 10,000 = 100,000,000개이므로 약 32만 개 스레드가 남는다. 행과 열 방향 각각 16개씩 경계에서 버려지는 스레드들이다.

이 스레드들은 커널 함수에서 if (row < nrow && col < ncol) 조건으로 걸러진다. 아무 일도 하지 않고 종료한다. 낭비처럼 보이지만 불가피하다. 블록 크기를 32로 고정하면 10,000은 32의 배수가 아니므로 남는 스레드가 생긴다. 32만 개 스레드는 전체의 0.3%다. 무시할 만한 수준이다.

c
dim3 dimBlock(32, 32, 1);
dim3 dimGrid((ncol + dimBlock.x - 1) / dimBlock.x,
             (nrow + dimBlock.y - 1) / dimBlock.y, 1);

dim3는 CUDA의 3차원 차원 타입이다. (x, y, z) 세 값을 가지며, 2차원 그리드이므로 z = 1이다.

CUDA 커널 구현과 인덱싱

matadd-dev.cu에서 CUDA kernel function을 작성한다.

c
__global__ void kernel_matadd(float* c, const float* a, const float* b,
                               unsigned nrow, unsigned ncol) {
    unsigned col = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned row = blockIdx.y * blockDim.y + threadIdx.y;

    if (row < nrow && col < ncol) {
        unsigned i = row * ncol + col;
        c[i] = a[i] + b[i];
    }
}

위 커널의 동작을 단계별로 살펴보자.

첫째, 각 스레드가 2D 인덱스를 계산한다. col = blockIdx.x * blockDim.x + threadIdx.x는 열 인덱스를 계산한다. blockIdx.x는 블록의 x 좌표(0부터 312까지), blockDim.x는 블록의 x 크기(32), threadIdx.x는 블록 내 스레드의 x 좌표(0부터 31까지)다. 예를 들어 블록 (5, 3)의 스레드 (10, 20)은 col = 5 * 32 + 10 = 170, row = 3 * 32 + 20 = 116이 된다. 이 스레드는 행렬의 (116, 170) 위치를 담당한다.

둘째, 범위 체크를 수행한다. if (row < nrow && col < ncol) 조건이 필요한 이유는 경계 스레드 때문이다. row나 col이 10,000 이상이면 행렬 범위를 벗어나므로 아무것도 하지 않고 종료한다.

셋째, 1D 인덱스로 변환한다. i = row * ncol + col은 row-major 레이아웃 공식이다. (row, col) 2D 좌표를 메모리 배열의 1D 인덱스로 변환한다.

넷째, 덧셈을 수행한다. c[i] = a[i] + b[i]는 간단하다. 세 배열 모두 같은 인덱스를 사용한다. Global memory에서 읽기 2번, 쓰기 1번이 발생한다.

커널 호출 코드를 보자.

c
int main(const int argc, const char* argv[]) {
    unsigned nrow = 10000;
    unsigned ncol = 10000;

    float* matA = new float[nrow * ncol];
    float* matB = new float[nrow * ncol];
    float* matC = new float[nrow * ncol];
    setNormalizedRandomData(matA, nrow * ncol);
    setNormalizedRandomData(matB, nrow * ncol);

    float *dev_matA, *dev_matB, *dev_matC;
    cudaMalloc((void**)&dev_matA, nrow * ncol * sizeof(float));
    cudaMalloc((void**)&dev_matB, nrow * ncol * sizeof(float));
    cudaMalloc((void**)&dev_matC, nrow * ncol * sizeof(float));

    cudaMemcpy(dev_matA, matA, nrow * ncol * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_matB, matB, nrow * ncol * sizeof(float), cudaMemcpyHostToDevice);

    ELAPSED_TIME_BEGIN(0);
    dim3 dimBlock(32, 32, 1);
    dim3 dimGrid((ncol + dimBlock.x - 1) / dimBlock.x,
                 (nrow + dimBlock.y - 1) / dimBlock.y, 1);
    kernel_matadd <<<dimGrid, dimBlock>>> (dev_matC, dev_matA, dev_matB, nrow, ncol);
    cudaDeviceSynchronize();
    CUDA_CHECK_ERROR();
    ELAPSED_TIME_END(0);

    cudaMemcpy(matC, dev_matC, nrow * ncol * sizeof(float), cudaMemcpyDeviceToHost);

    printMat("matC", matC, nrow, ncol);

    cudaFree(dev_matA);
    cudaFree(dev_matB);
    cudaFree(dev_matC);
    delete[] matA;
    delete[] matB;
    delete[] matC;
}

위 코드는 표준적인 CUDA 패턴을 따른다. cudaMalloc으로 디바이스 메모리를 할당하고, 세 행렬 모두 nrow * ncol * sizeof(float) = 400MB가 필요하다. cudaMemcpy로 입력 행렬 A, B를 호스트에서 디바이스로 복사한다. 커널 실행 후 결과 행렬 C를 디바이스에서 호스트로 복사한다.

커널 실행은 <<<dimGrid, dimBlock>>> 구문을 사용하며 313×313313 \times 313 그리드와 32×3232 \times 32 블록으로 실행한다. cudaDeviceSynchronize()로 커널 완료를 기다린다. CUDA 커널은 비동기로 실행되므로 이 함수를 호출해야 GPU가 끝날 때까지 CPU가 기다린다.

GeForce RTX 2070에서 3,282 usec가 소요된다. 약 3.3ms다. CPU 버전이 168,857 usec(169ms)였으므로 168.857/3.282≈51168.857 / 3.282 \approx 51배 빠르다.

메모리 최적화

Memory Coalescing의 중요성

Memory coalescing은 GPU 성능의 핵심이다. 여러 스레드의 메모리 접근을 하나의 트랜잭션으로 합치는 기법이다.

Warp는 32개 스레드를 가진다. 이 32개 스레드가 메모리 접근을 할 때, 만약 연속된 주소에 접근한다면 하나의 메모리 트랜잭션으로 처리할 수 있다. 예를 들어 warp의 스레드 0이 주소 100을 읽고, 스레드 1이 주소 104를 읽고, 스레드 31이 주소 224를 읽는다고 하자. 32개 float, 128바이트를 연속으로 읽는다. GPU는 이것을 하나의 128바이트 트랜잭션으로 합친다. 32번 메모리 접근 대신 1번만 한다. 32배 효율적이다.

반대로 스레드들이 흩어진 주소에 접근하면 어떻게 될까? 각 스레드마다 별도의 메모리 트랜잭션이 필요하다. 32번 메모리 접근을 해야 한다. 32배 느리다.

우리 커널에서 메모리 접근 패턴을 보자. Warp는 연속된 32개 스레드로 구성된다. 2D 블록 (32, 32)에서 첫 번째 warp는 threadIdx.y = 0, threadIdx.x = 0~31이다. col = blockIdx.x * 32 + 0~31, row = blockIdx.y * 32 + 0이므로 같은 행의 연속된 열들이다.

인덱스 i = row * ncol + col을 계산하면 다음과 같다.

  • 스레드 0: row * ncol + col₀
  • 스레드 1: row * ncol + col₀ + 1
  • 스레드 2: row * ncol + col₀ + 2
  • ...
  • 스레드 31: row * ncol + col₀ + 31

연속된 메모리 주소다. Memory coalescing이 완벽하게 작동한다.

다음 warp는 threadIdx.y = 1, threadIdx.x = 0~31이다. 다음 행의 연속된 열들이며 여전히 연속된 메모리 주소다. 결과적으로 모든 warp가 연속된 메모리에 접근한다. Coalescing 효율이 100%다. 이것이 높은 메모리 대역폭을 달성한 이유다.

잘못된 Memory 접근 패턴

만약 row와 col을 바꾸면 어떻게 될까? 다음 커널을 보자.

c
__global__ void kernel_matadd_wrong(float* c, const float* a, const float* b,
                                     unsigned nrow, unsigned ncol) {
    unsigned row = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned col = blockIdx.y * blockDim.y + threadIdx.y;

    if (row < nrow && col < ncol) {
        unsigned i = row * ncol + col;
        c[i] = a[i] + b[i];
    }
}

row = blockIdx.x * blockDim.x + threadIdx.x로 바뀌었고, col = blockIdx.y * blockDim.y + threadIdx.y로 바뀌었다. 한 warp 내 스레드들은 threadIdx.x = 0~31, threadIdx.y = 0으로 고정이다. row = blockIdx.x * 32 + 0~31, col = blockIdx.y * 32 + 0이므로 연속된 행, 같은 열이다.

인덱스 i = row * ncol + col을 계산하면 다음과 같다.

  • 스레드 0: 0 * ncol + col
  • 스레드 1: 1 * ncol + col
  • 스레드 2: 2 * ncol + col
  • ...
  • 스레드 31: 31 * ncol + col

ncol = 10,000이므로 각 스레드의 메모리 주소는 10,000 간격이다. 40,000 바이트 간격이다. 완전히 흩어져 있다.

Memory coalescing이 불가능하다. 각 스레드마다 별도의 메모리 트랜잭션이 필요하다. 32배 더 많은 메모리 접근이 발생한다. GeForce RTX 2070에서 11,046 usec가 소요된다. 정상 버전이 3,282 usec였으므로 11.046/3.282≈3.411.046 / 3.282 \approx 3.4배 느려진다.

왜 32배가 아니라 3.4배일까? 몇 가지 완화 요인이 있다.

첫째, GPU 캐시가 일부 접근을 빠르게 처리한다. L2 캐시는 모든 SM이 공유하므로 여러 warp가 비슷한 영역에 접근하면 캐시 히트가 발생한다.

둘째, 메모리 컨트롤러가 요청을 재정렬해서 일부 효율을 회복한다.

셋째, 계산이 전혀 없는 것은 아니다. 덧셈 연산 자체에도 시간이 걸리므로 메모리 접근이 3배 느려져도 전체 실행 시간은 3배 미만으로 증가한다.

그래도 3.4배 차이는 크다. Memory coalescing의 중요성을 명확히 보여준다.

관련 문서