CUDAとOpenCLによるGPUを用いた行列-行列積(GEMM)の実装テスト備忘録。OpenMPを用いたCPUのマルチスレッド型の並列実装とも比較する
実装コード
実装したコード、環境構築用スクリプトはGitHubリポジトリ参照
評価環境
| CPU | Intel Core i7-9700K |
| GPU | NVIDIA GeForce RTX 2060 Super |
| DDR Memory | DDR 16GB 2666MHz |
| OS | Ubuntu 24.04.4 LTS |
| Compiler | CUDA: nvcc V12.0.140 OpenCL: gcc (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0 |
| Option | CUDA: -O2 -arch=sm_75 OpenCL: -O2 -lOpenCL -mcmodel=medium |
| Time measure | CUDA: cudaEventRecord OpenCL: clGetEventProfilingInfo |
環境構築
詳細はGithub上の "env"ディレクトリと、READMEファイル参照
パッケージリスト更新等
sudo apt update -y
sudo apt upgrade -y
sudo apt autoremove -y
GPU認識確認コマンド：
lspci | grep -i nvidiaCUDA環境構築
ドライバインストール コマンド
ubuntu-drivers devicessudo apt install nvidia-driver-580-open再起動：
sudo rebootsudo apt install nvidia-cuda-toolkit＜環境確認コマンド＞
ドライバ確認：
ubuntu-drivers devicesnvidia-sminvcc --versionOpenCL
ドライバーインストール：
sudo apt install nvidia-driver-580-open nvidia-opencl-dev -ysudo apt install clinfo -y＜環境確認コマンド＞
プラットフォーム確認：
clinfoCPU-GPU間のデータ受け渡し
CUDA
GPUメモリの動的確保
cudaMalloc(void** devPtr, size_t size)cudaMemcpy(void* dst , const void* src , size_t size , cudaMemcpyKind cudaMemcpyHostToDevice)cudaMemcpy(void* dst , const void* src , size_t size , cudaMemcpyKind cudaMemcpyDeviceToHost)OpenCL
GPUメモリ == Bufferオブジェクト
として表現されているため、下記キュー型の実行を行う
clCreateBuffer
↓
clEnqueueWriteBuffer
↓
kernel実行
↓
clEnqueueReadBuffer
GPUメモリの動的確保
cl_mem clCreateBuffer( cl_context context, cl_mem_flags flags, size_t size, void* host_ptr, cl_int* errcode_ret )cl_int clEnqueueWriteBuffer( cl_command_queue command_queue, cl_mem buffer, cl_bool blocking_write, size_t offset, size_t size, const void* ptr, cl_uint num_events_in_wait_list, const cl_event* event_wait_list, cl_event* event )cl_int clEnqueueReadBuffer( cl_command_queue command_queue, cl_mem buffer, cl_bool blocking_read, size_t offset, size_t size, void* ptr, cl_uint num_events_in_wait_list, const cl_event* event_wait_list, cl_event* event )GPUアーキテクチャの勉強
評価に使用したRTX2060 Superは、2018年後半のNvidia Turing（チューリング）アーキテクチャを採用しています。ですので、今回はこのアーキテクチャをベースとして、CUDA実装をベースにGPUの構成を見ていきます。
参照：
- https://developer.nvidia.com/blog/nvidia-turing-architecture-in-depth/
- https://www.nvidia.com/content/dam/en-zz/Solutions/design-visualization/technologies/turing-architecture/NVIDIA-Turing-Architecture-Whitepaper.pdf
GPUアーキテクチャの構成要素
CUDAで高性能な行列演算を実装するために、次の要素が重要
- CUDA Core
- Warp
- Shared Memory
- Tensor Core
- WMMA API
GPU アーキテクチャの基本構造
GPUは以下の階層構造で構成される。
GPU
├ GPC (Graphics Processing Cluster)
│ ├ TPC (Texture Processing Cluster)
│ │ ├ SM (Streaming Multiprocessor)
│ │ │ ├ CUDA Core
│ │ │ ├ Tensor Core
│ │ │ ├ Shared Memory
│ │ │ └ Warp Scheduler
CUDAプログラムは最終的に SM (Streaming Multiprocessor) 上で実行される。
Turing世代のSMは以下のような構成
- 64 FP32 CUDA cores
- 64 INT32 cores
- 8 Tensor Cores
- Warp Scheduler
Turingでは FP32 と INT32 を同時実行可能で、浮動小数演算と整数演算を並列に処理可能
Warp（GPUスレッドの実行単位）
CUDAではスレッドは1つずつ実行されるのではなく、
Warp単位で実行
1 Warp = 32 threads
今回のコードも Warp を前提とした実装
int tid = threadIdx.x;
int warpId = tid / 32;
int laneId = tid % 32;
| 変数 | 意味 |
|---|---|
| tid | スレッドID |
| warpId | Warp番号 |
| laneId | Warp内スレッドID |
GPUは32スレッドを同時に実行するSIMD方式のため、
Warpの概念はCUDA最適化の基本
メモリ階層
CUDAの性能はメモリアクセスに大きく依存
主なメモリ階層
Global Memory (GPU DRAM)
Shared Memory (SM内部)
Register
速度の順序
Register > Shared Memory > Global Memory
そのためCUDA最適化の基本戦略は
Global → Shared → Register
という順序でデータを移動すること
今回のコードでも Shared Memory を利用
__shared__ half As[BLOCK_M][BLOCK_K];
__shared__ half Bs[BLOCK_K][BLOCK_N];
Shared Memory Tiling
GPUで行列演算を高速化する代表的なテクニックの Tiling
行列積
C = A × B
を小さなブロックに分割して計算
今回のコードでは次のタイルサイズを使用
BLOCK_M = 128
BLOCK_N = 128
BLOCK_K = 16
行列は以下のように分割
C matrix
+---------+---------+
| block00 | block01 |
+---------+---------+
| block10 | block11 |
+---------+---------+
Tensor Core
Tensor Core は行列演算専用のハードウェアユニット
Tensor Core は次の演算を高速に実行
D = A × B + C
これは FMA (Fused Multiply Add) と呼ばれる。
Turing GPU の Tensor Core は
16 × 16 × 16
サイズの行列演算を1命令で処理可能。
Turing GPUでは
- 1 SM に 8 Tensor Core
- Tensor Core は FP16 FMA を高速実行
という構成
WMMA (Warp Matrix Multiply Accumulate)
CUDAからTensor Coreを利用するAPIが
WMMA
CUDAでは
nvcuda::wmma
という名前空間で提供されている
今回のコードでも使用
wmma::fragment<wmma::matrix_a, WMMA_M, WMMA_N, WMMA_K, half, wmma::row_major> a_frag;
wmma::fragment<wmma::matrix_b, WMMA_M, WMMA_N, WMMA_K, half, wmma::row_major> b_frag;
wmma::fragment<wmma::accumulator, WMMA_M, WMMA_N, WMMA_K, float> c_frag;
fragment は Tensor Core 用のレジスタ行列
Tensor Core 演算の流れ
Tensor Core演算は次の流れで実行
Global Memory
↓
Shared Memory
↓
WMMA fragment
↓
Tensor Core
↓
Accumulator
今回のコードでも以下の3ステップで実行。
wmma::load_matrix_sync(a_frag,&As[warpRow][k],BLOCK_K);
wmma::load_matrix_sync(b_frag,&Bs[k][warpCol],BLOCK_N);
wmma::mma_sync(c_frag,a_frag,b_frag,c_frag);
Warp Tiling
Tensor Coreの性能を最大化するために 計算を以下の階層で分割
block
└ warp
└ tensor core
今回のコードでは
- block = 128 × 128
- warp = 64 × 64
- tensor core = 16 × 16
という構造で計算を分割
Mixed Precision (混合精度)
Tensor Coreでは Mixed Precision が使用される。
FP16 × FP16 → FP32
今回のコードでも
half *a
half *b
float *c
という構成
これは
- FP16で高速計算
- FP32で精度保持
を両立するための設計
CUDA Kernel Execution
Kernelの起動は次の形式
tensor_gemm<<<blocks,threads>>>(dA,dB,dC);
今回の設定は
threads = 128
block = 128×128 tile
つまり
1 block = 4 warps
という構成
まとめ
今回実装したCUDA Tensor Core GEMMコードを理解するために、 以下の概念が重要
- GPU階層構造 (SM)
- Warp execution
- メモリ階層 (Global / Shared)
- Tiling
- Tensor Core
- WMMA API
- Mixed Precision
- Warp tiling
これらを組み合わせることで GPUでCPUより高速な行列演算を実現可能
評価に用いる行列-行列積(GEMM)
こちらでOpenMPで評価している行列-行列積と同じ演算を行います。
Kernel部分の実装部分
全ての要素に値を入れた密な行列AとBを入力として行列の積（General Matrix Multiply）を計算し、結果として行列Cを出力する処理を考えます。例を示すと以下のような処理です。
前提条件として「行列Ａの列数＝行列Ｂの行数」となっていることが挙げられます。得られる結果の行列Ｃの大きさは、「行×列」＝「行列Ａの行数×行列Ｂの列数」です。
簡易速度評価
演算精度や、評価環境など厳密に条件をそろえて比較等しておりません。簡易的なテストだと思ってみてください。
実装環境
CUDA
makefile:
matrix_gemm: ../src/matrix_gemm.cu
nvcc ../src/matrix_gemm.cu -O2 -arch=sm_75 -o matrix_gemm時間測定
GPU上でのカーネル実行時間を測定しています。
CPU-GPU間のメモリ転送やCPU側の処理時間は含めません。
float total_ms=0.0f;
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start);
// kernel launch
tensor_gemm<<>>(dA,dB,dC);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float ms;
cudaEventElapsedTime(&ms,start,stop);
total_ms+=ms;
cudaEventDestroy(start);
cudaEventDestroy(stop); cudaEventRecord()で、GPUのコマンドキューへイベントを追加しています。
cudaEventRecord()のstartと、stopの間にメモリ転送を含めると転送時間を含む時間の測定をすることも可能です。
OpenCL
makefile:
matrix_gemm: ../src/matrix_gemm.cpp
g++ ../src/matrix_gemm.cpp -O2 -lOpenCL -mcmodel=medium -o matrix_gemm時間測定
double total_host_ms = 0.0; // CPU上の時間
double total_device_ms = 0.0; // GPU内部の時間
cl_event event;
auto start = std::chrono::high_resolution_clock::now(); // CPU測定START
err = clEnqueueNDRangeKernel(queue, kernel, 2, NULL, global, local, 0, NULL, &event);
clFinish(queue);
auto end = std::chrono::high_resolution_clock::now(); // CPU測定END
total_host_ms += std::chrono::duration(end - start).count();
// Device時間取得 (GPU)
cl_ulong start_time, end_time;
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &start_time, NULL);
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(cl_ulong), &end_time, NULL);
total_device_ms += (end_time - start_time) * 1.0e-6;
clReleaseEvent(event); 時間測定
CPUで実行するSerial版や、OpenMP版と比べると、GPUで実行しているCUDAや、OpenCL版のほうが高速に演算できています。
CUDAと比較してOpenCLでそれほど性能を出せていません。もう少し改良の余地はありそうです。
参考：
