행렬곱 병렬화, 타일링 알고리즘

큰 행렬을 작은 블록으로 나누어 각 블록을 병렬 처리한다.

각 블록은 결과 행렬의 부분을 계산하며, 내부적으로 공유 메모리를 사용해 효율성을 높인다.

이 방식은 타일 단위로 데이터를 공유 메모리에 로드하여 전역 메모리 접근을 최소화하며, 계산 시 반복해서 사용하는 데이터를 공유 메모리에서 읽어오도록 최적화한다.

 

1. 행렬 분할

행렬 A, B 를 각 블록 크기에 맞춰 여러 개의 작은 타일로 나눈다. 각 블록은 행렬 A 와 B 의 한 타일을 가져와, 결과 행렬 C 의 부분 결과를 계산한다.

2. 공유 메모리를 사용한 타일 연산

각 타일의 데이터를 전역 메모리에서 공유 메모리로 로드한다. 이는 GPU 에서 빠르게 접근할 수 있어, 전역 메모리 접근을 최소화한다. 

각 스레드는 공유 메모리에 있는 데이터만을 사용해 타일 내 연산을 수행하여, 최종적으로 C 의 한 타일에 대한 결과를 계산한다. 

3. 동기화

타일 연산이 완료될 때마다 __syncthreads() 를 호출해 블록 내의 모든 스레드가 같은 단계로 이동

이를 통해 모든 스레드가 다음 타일을 읽기 전에 현재 타일의 연산이 끝났음을 보장한다.

#include <cuda_runtime.h>
#include <iostream>

#define TILE_SIZE 16  // 타일 크기
#define N 1024  // 행렬 크기

__global__ void matrixMulTiled(float *A, float *B, float *C, int n) {
    __shared__ float shared_A[TILE_SIZE][TILE_SIZE];
    __shared__ float shared_B[TILE_SIZE][TILE_SIZE];

    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float value = 0;

    for (int i = 0; i < n / TILE_SIZE; i++) {
        shared_A[threadIdx.y][threadIdx.x] = A[row * n + i * TILE_SIZE + threadIdx.x];
        shared_B[threadIdx.y][threadIdx.x] = B[(i * TILE_SIZE + threadIdx.y) * n + col];
        __syncthreads();  // 타일이 모두 로드될 때까지 대기

        for (int j = 0; j < TILE_SIZE; j++) {
            value += shared_A[threadIdx.y][j] * shared_B[j][threadIdx.x];
        }
        __syncthreads();  // 현재 타일의 계산이 끝나면 다음 타일로 넘어가기 위해 대기
    }

    if (row < n && col < n) {
        C[row * n + col] = value;
    }
}

int main() {
    int size = N * N * sizeof(float);
    float *h_A = (float*)malloc(size);
    float *h_B = (float*)malloc(size);
    float *h_C = (float*)malloc(size);

    for (int i = 0; i < N * N; i++) {
        h_A[i] = 1.0f;
        h_B[i] = 1.0f;
    }

    float *d_A, *d_B, *d_C;
    cudaMalloc((void**)&d_A, size);
    cudaMalloc((void**)&d_B, size);
    cudaMalloc((void**)&d_C, size);

    cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

    dim3 dimBlock(TILE_SIZE, TILE_SIZE);
    dim3 dimGrid(N / TILE_SIZE, N / TILE_SIZE);
    matrixMulTiled<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, N);

    cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

    // 일부 결과 출력
    std::cout << "C[0]: " << h_C[0] << " C[N*N-1]: " << h_C[N * N - 1] << std::endl;

    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    free(h_A);
    free(h_B);
    free(h_C);

    return 0;
}

공유 메모리 할당 : shared_A,B 는 각 타일의 데이터를 저장할 공유 메모리 배열, 16 x 16 크기로 설정, 

블록과 스레드 구성 : 행렬 크기와 타일 크기를 바탕으로 블록과 그리드 구성

타일링을 사용한 연산의 목표는

전역 메모리 접근 최소화, 공유 메모리 활용, 데이터 로컬리티 향상

 

결과 행렬, C[i][j] 를 계산하기 위해 A, B 의 여러 요소에 접근해야 하므로, 큰 행렬에서는 전역 메모리 접근이 많아 성능 저하가 발생할 수 있다.

타일링을 이용한 최적화

1. 타일 단위로 분할

2. 공유 메모리에 로드

3. 부분 결과 계산

4. 동기화 및 다음 타일로 이동

타일링 유무의 속도 비교

#include <cuda_runtime.h>
#include <iostream>
#include <chrono>

#define TILE_SIZE 16  // 타일 크기
#define N 1024        // 행렬 크기

// 타일링을 사용하지 않은 행렬 곱셈 커널
__global__ void matrixMulNaive(float *A, float *B, float *C, int n) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    float value = 0;
    if (row < n && col < n) {
        for (int k = 0; k < n; k++) {
            value += A[row * n + k] * B[k * n + col];
        }
        C[row * n + col] = value;
    }
}

// 타일링을 사용한 행렬 곱셈 커널
__global__ void matrixMulTiled(float *A, float *B, float *C, int n) {
    __shared__ float shared_A[TILE_SIZE][TILE_SIZE];
    __shared__ float shared_B[TILE_SIZE][TILE_SIZE];

    int row = blockIdx.y * TILE_SIZE + threadIdx.y;
    int col = blockIdx.x * TILE_SIZE + threadIdx.x;
    float value = 0;

    for (int i = 0; i < n / TILE_SIZE; i++) {
        shared_A[threadIdx.y][threadIdx.x] = A[row * n + i * TILE_SIZE + threadIdx.x];
        shared_B[threadIdx.y][threadIdx.x] = B[(i * TILE_SIZE + threadIdx.y) * n + col];
        __syncthreads();

        for (int j = 0; j < TILE_SIZE; j++) {
            value += shared_A[threadIdx.y][j] * shared_B[j][threadIdx.x];
        }
        __syncthreads();
    }

    if (row < n && col < n) {
        C[row * n + col] = value;
    }
}

int main() {
    int size = N * N * sizeof(float);
    float *h_A = (float*)malloc(size);
    float *h_B = (float*)malloc(size);
    float *h_C_naive = (float*)malloc(size);
    float *h_C_tiled = (float*)malloc(size);

    for (int i = 0; i < N * N; i++) {
        h_A[i] = 1.0f;
        h_B[i] = 1.0f;
    }

    float *d_A, *d_B, *d_C;
    cudaMalloc((void**)&d_A, size);
    cudaMalloc((void**)&d_B, size);
    cudaMalloc((void**)&d_C, size);

    cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

    dim3 dimBlock(TILE_SIZE, TILE_SIZE);
    dim3 dimGrid(N / TILE_SIZE, N / TILE_SIZE);

    // 타일링을 사용하지 않은 행렬 곱셈
    auto start_naive = std::chrono::high_resolution_clock::now();
    matrixMulNaive<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, N);
    cudaDeviceSynchronize();
    auto end_naive = std::chrono::high_resolution_clock::now();
    std::chrono::duration<float, std::milli> naive_duration = end_naive - start_naive;
    cudaMemcpy(h_C_naive, d_C, size, cudaMemcpyDeviceToHost);

    // 타일링을 사용한 행렬 곱셈
    auto start_tiled = std::chrono::high_resolution_clock::now();
    matrixMulTiled<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, N);
    cudaDeviceSynchronize();
    auto end_tiled = std::chrono::high_resolution_clock::now();
    std::chrono::duration<float, std::milli> tiled_duration = end_tiled - start_tiled;
    cudaMemcpy(h_C_tiled, d_C, size, cudaMemcpyDeviceToHost);

    // 성능 비교 출력
    std::cout << "Naive matrix multiplication duration: " << naive_duration.count() << " ms" << std::endl;
    std::cout << "Tiled matrix multiplication duration: " << tiled_duration.count() << " ms" << std::endl;

    // 메모리 해제
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    free(h_A);
    free(h_B);
    free(h_C_naive);
    free(h_C_tiled);

    return 0;
}

>>>
Naive matrix multiplication duration: 22.2422 ms
Tiled matrix multiplication duration: 1.6547 ms

약 15배의 속도 차이 확인 가능

 

(4,4) 타일을 이용한 (8,8) 행렬 곱셈

행렬 A, B 를 (4,4) 타일로 분할하여, (2,2) 크기의 그리드로 나눈다. 

첫 번째 타일에서 A[0:4][0:4] 와 B[0:4][0:4] 를 곱하여 C[0:4][0:4] 의 부분 결과 계산

두 번째 타일에서 A[0:4][4:8] 와 B[4:8][0:4] 을 곱하여  C[0:4][0:4] 에 추가로 더해준다.