cuBLAS MatMul Tutorial

cuBLAS 소개

cuBLAS[1]는 NVIDIA CUDA runtime에서 돌아가는 BLAS를 구현한 library이다. BLAS는 Basic Linear Algebra Subprograms의 약자로 일반적으로 많이 사용되는 선형 대수의 식을 정리한 스펙(speification)[2]이다. Vector-vector 연산, Matrix-vector 연산 Matrix-matrix 연산 등의 많이 사용되는 선형 대수를 표준화한 것으로 입출력과 루틴을 정의해놓았다. BLAS 스펙을 구현한 예로는 cuBLAS(NVIDIA),  Intel MKL(Intel), clBLAS(Open source), OpenBLAS(Open source) 등이 있는데, 선형 대수는 많은 애플리케이션에서 사용되기에 많은 accelerator 제조사들은 자기들 제품을 위한 BLAS library를 제공해준다. 각 제조사는 각자의 하드웨어에 최적화된 library를 제공하기에 계산 성능이 매우 좋다. 특히 같은 제조사 내의 각 세대별 제품에도 따로 최적화를 진행하기에 성능을 항상 믿을 수 있다. 따라서 선형대수 관련 알고리즘을 직접 병렬 구현해야 할 일이 있다면, 직접 구현하는 것보다는 하드웨어 제조사에서 제공하는 BLAS library를 쓰는 것이 먼저 시도돼야 할 병렬화 전략이라고 생각한다.

오늘은 NVIDIA의 CUDA runtime에서 사용할 수 있는 BLAS implementation인 cuBLAS에 대하여 실제 사용가능한 예제를 통해 정리하고자 한다.

Fig 1. cuBLAS logo

예제는 FP32 type matrix multiplication이며 행렬 A[M,K], B[K,N], C[M,N]에 대하여 C = A*B를 계산한다. cuBLAS의 single precision floating point GEMM (general matrix multiplication) implementation인 cublasSgemm를 이용하여 행렬곱을 계산할 것이다.

 

 

 

GEneral Matrix Multiplication

BLAS가 제공하는 선형 대수 Funtionality는 3 단계(각각 level 1, level 2, level 3라고 함)로 나뉠 수 있다. Level 1은 vector-vector 연산을, level 2는 matrix-vector 연산을, level 3는 matrix-matrix 연산을 의미한다. cuBLAS는 level 3에 속한 연산 중 하나인 GEMM(general matrix multiplication)구현을 제공한다.

행렬 A, B, C와 스칼라 α, β에 대하여 GEMM의 specification은 다음과 같다.

Eq 1. GEMM의 일반적인 형식

α = 1, β = 0으로 설정할 경우, GEMM을 통해 A, B의 행렬곱을 구할 수 있다.

 

 

 

준비

cuBLAS는 CUDA toolkit 설치 시 설치되기에, cuDNN과 같이 따로 설치할 필요가 없다. cuBLAS는 library이기 때문에 header file과 library file(linux에서는 .so파일 또는 .a 파일, Windows에서는 .lib 파일 또는 dll 파일)로 존재한다. 설치 위치는 CUDA toolkit이 설치된 디렉토리를 찾아보면 된다(나의 경우, "/usr/local/cuda/include/", "/usr/local/cuda/lib64"에 각각 존재).

Library 이기에 코드 작성 시 header file만 잘 include 시키고, 컴파일과 동적 링킹 시 library file만 잘 링킹해주면 된다. 이후 cuBLAS documentation[1]를 참고하여 잘 구현된 선형대수 함수를 사용하면 된다.

 

 

 

FP32 Matrix Multiplication Example

다음 Github repository에 전체 코드를 추가하였다.

GitHub - jhson989/matmul_cublas: cuBLAS GEMM Example for FP32 MatMul

cuBLAS를 사용하기 위해서는 cublasHandle을 정의해주어야 한다.

// cublas handler 생성
cublasHandle_t handle;
cublasCreate (&handle);

// cuBLAS routine 실행
...

// cublas handler 해제
cublasDestroy(handle);

FP32 GEMM 계산을 위하여 cublasSgemm 함수를 실행한다.

// cublasStatus_t cublasSgemm(cublasHandle_t handle,
//                           cublasOperation_t transa, cublasOperation_t transb,
//                           int m, int n, int k,
//                           const float           *alpha,
//                           const float           *A, int lda,
//                           const float           *B, int ldb,
//                           const float           *beta,
//                           float           *C, int ldc)

// A, B, C or B, A, C
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N)

 

cublasSgemm은 C = alpha*A*B + beta*C 계산을 하는데, alpha를 1, beta를 0으로 설정하면 C = A*B 행렬곱 연산을 수행한다. 이때 cublasSgemm definition에서는 A, B, C 순서로 파라미터를 받는데 실제 cublasSgemm을 쓴 부분에는 B, A, C 순서로 파라미터를 입력한 것을 알 수 있다. 이렇게 A와 B의 순서를 바꾼 이유는 A, B, C를 row-major로 저장하였는데, cublasSgemm는 column-major 행렬을 입력으로 받기 때문이다.

 

 

 

Row-major VS. Column-major

행렬은 2차원 배열이다. 하지만 실제 메모리에 저장될 때는 1차원 주소 공간에 저장된다. 따라서 저장 방식은 row-major, column-mojor 2가지 방식이 있을 수 있다. Row-major는 row 방향으로 이웃한 데이터가 연속되게 저장되는 것이고, column-major는 column 방향으로 이웃한 데이터가 연속되게 저장되는 방식이다. Fig 2는 두 차이를 잘 보여주는데, 선이 연결된 순서로 데이터가 저장된다. Row-major의 경우 a11, a12, a13, a21, a22, a23, a31, a32, a33 순서로 실제 메모리에 저장된다.

Fig 2. row-major와 column-major 비교. [3]

cuBLAS의 함수는 column-major 방식의 matrix를 가정한다. 일반적으로 우리는 row-major matrix를 많이 사용하는데, 그렇기에 이것을 column-major로 변환시키거나 행렬곱의 순서(A*B -> B*A)를 바꾸는 트릭을 사용하면 된다.

Fig 3. row-major와 column-major의 관계

Fig 3.은 순서를 바꾸는 트릭이 어떻게 작동하는 지를 보여준다. Row-major에서의 A(그림의 왼쪽)과 column-major에서의 A transpose(그림의 오른쪽)는 사실 메모리 상에 같은 형태로 저장된다. 따라서 column-major에서의 C transpose를 구하는 것은 row-major상의 C를 구하는 것과 같다.

Eq 2.

따라서 Eq 2.와 같이 C transpose는 B transpose와 A transpose를 곱하면 구할 수 있고, column-major에서 C transpose를 구한다는 것은 row-major에서 C를 구하는 것과 동일하기에 우리가 원하는 row-major C를 얻을 수 있다. Column-major에서의 B transpose와 A transpose는 모두 row-major에서의 B와 A와 완전히 일치하게 메모리에 저장되기에 단순히 B*A만 해도 C를 구할 수 있다.

 

 

 

Reference

[1] https://docs.nvidia.com/cuda/cublas/index.html

[2] https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms

[3] https://en.wikipedia.org/wiki/Row-_and_column-major_order