この記事は DeNA 23 新卒 Advent Calendar 2022 の24日目の記事です。

はじめに

突然ですがみなさんはGPUは知っているでしょうか。自分はGPUに関する研究をしているのですが、自分の研究を説明してGPUって何？と聞かれると少し悲しくなります。GPUとはGraphics Processing Unitの略であり、画像処理に特化した演算装置です。しかし、近年はGPGPUと言われる一般用途の計算にGPUを活用する動きが高まっています。2022年末現在、暗号資産用GPUの需要低下からGPUの価格も安定しています。また、GPUが付属しているゲーム用PCを購入する人の増加により、GPUは持っているが画像処理にしか使っていない人も多いです。

個人的にGPUは色々出来て面白いので、GPUプログラミングを出来る人がもっと増えれば、使われていないGPUをより有効的に活用できるのではないかと思い、今回記事を執筆しようと考えました。本記事ではGPUプログラミングの詳細部分にはあまり触れすぎず、読者にGPUプログラミングに興味を持ってもらうことを目的とします。また、今回対象とするGPUはNVIDIA社のGPUとします。

どこでGPUを使うのか

GPUを使えば計算はなんでも速くなると勘違いされがちですが、そういうわけではありません。得意なものと苦手なものがあります。

まず、並列処理可能な計算である必要があります。GPGPUでは計算処理をGPUに存在するコアを用いて並列に実行することで処理が高速化していますが、そもそも並列に実行できない計算の場合はGPUを使用しても意味がありません。計算した結果を次の計算で使用する場合などは、GPUを使用しても逐次処理となり、CPUと実行時間が変わらないことが多いです。むしろデータをGPUとやり取りする分遅くなる可能性もあります。

つぎに、同じデータを何回も使用する計算が得意です。計算量が多くても、データ量も多いとGPUとホストとのデータやり取りの時間がネックになってしまいます。少ないデータを何回も使用した簡単な計算がGPUの得意とする計算です。

よくGPUが得意とする計算の例として挙げられるのが行列積や連立方程式の解法であるヤコビ法などです。とはいえあくまで例なので、単純に3重ループなどが存在するプログラムなどを見たら、「GPU使えないかなあ」と考えてみるとGPUを使いたい場所が見つかるかもしれません。

GPUプログラミング環境

GPUプログラミングではNVIDIA社が開発した並列コンピューティングプラットフォームであるCUDA[1]を使うことが多いです。一般的にはC言語版のCUDA Cが使われていますが、Fortran版もあります。CUDA CはC言語の拡張となっているので、C言語が使用できれば逐次プログラムならすぐに作成することが可能です。しかし、並列計算を行うためにはGPUの構造を理解して記述を行う必要があります。

CUDAとは別にOpenAcc[2]というものも存在します。CUDAとは異なり、C言語やFortranで作成された既存のプログラムに対して指示文を挿入することでGPU上で計算コードを書き換えることが可能です。ただ高速化可能なプログラム部分は自分で判断する必要があり、挿入する指示文も自分で考える必要があります。CPUで書かれた既存コードをGPU版に移行する際にGPUコーディングのコストカットのために利用されることが多いです。高速化する場所や挿入する指示文によって高速化具合は変わりますが、単純なものであればCUDAで1から書いたものに近い性能が出ることもあります。

単純にGPUを使用するという意味であれば、Pythonで簡単にGPUを使用できる汎用配列計算ライブラリであるCuPy[3]などもあります（GPUプログラミングとは少し違いますが）。PyTorch[4]などの人工知能向けフレームワークでもGPUが使用されています。

細かな最適化などもしたい場合やプログラムがまだ存在しない場合は1からCUDAで記述、既存のプログラムがすでにある場合はOpenAccを使用、特定分野に特化したプログラムであればフレームワークやライブラリを使用するというのが多いパターンだと思います。ライブラリならともかく、OpenAccなどはGPUのことを理解していなければ十分な性能を発揮するのが難しいため、最初は1からGPUプログラミングをすることがおすすめです。

この記事ではCUDAを用いたGPUプログラミングを主に扱っていきます。

CUDAにおけるスレッドの階層構造

プログラムをCUDAによるGPU化をするためには、GPU上で動かしたい部分を決め、GPUカーネル関数と呼ばれる別関数にする必要があります。GPUカーネル関数を用いて計算処理を高速化するために多くのスレッドを利用し、そのスレッドたちをグリッドと呼びます。グリッドは複数のスレッドブロック（以下ブロック）で構成され、ブロックは複数のスレッドによって構成されています。カーネル関数の呼び出し時にブロックの数と1ブロックあたりのスレッド数（以下ブロックサイズ）を決めることで、全体のスレッド数が決定されます。

下記のコードの場合N/BS個のブロック×ブロックサイズ＝N個のスレッドでGPUカーネル関数が実行されます。

func <<<(N/BS,BS)>>>(A);

グリッド、ブロックはともに3次元的に配置されており、x方向、y方向、z方向それぞれにスレッドが存在します。多次元グリッド、多次元ブロックにしたい場合は下記のように記述することが多いです。

dim3 dimGrid(N/BS, N/BS);
dim3 dimBlock(BS, BS);
func<<<dimGrid, dimBlock>>>(A);

上記のコードの場合、N/BS×N/BS個のブロックが生成され、ブロックサイズがBS×BS、合計でN×N個のスレッドが生成されます。この場合BSは1スレッドブロックあたりのスレッド数ではなく、正方形に形成されたスレッドブロックの1辺の長さにすぎないため注意が必要です。

各次元方向のグリッドサイズ、ブロックサイズ、1ブロックあたりのスレッド数には制限があり、CUDAインストール時のサンプルプログラムにあるファイルを実行すると確認できます。特に注意しなければならないのは1ブロックあたりのスレッド数が1024であることです。1次元ブロックであれば単純にブロックサイズを1024以下にする必要がありますが、2次元などの場合はスレッドブロックの各辺を掛け合わせて1024以下にする必要があるので注意が必要です。

このブロック数とブロックサイズによって性能が多少変わることもあります。まず、全体でどれくらいのスレッド数を実行するのかを決めましょう。全体のスレッド数は困ったらデータ数と同じだけのスレッド数、もしくはそれ以上にしましょう。多くて困ることはあまりありません。

最適なブロックサイズを決めるのはどのようなGPUカーネル関数を記述しているか、シェアードメモリなどを用いた特定の最適化手法を利用するかなどで決まるため、難しい問題です。しかし、CUDAではスレッドブロックごとのスレッド数の最大数が1024となっており、これ以上大きくすることは出来ません。そのため１ブロックあたりのスレッド数はとりあえず1024にすることが多いです。経験的に最適なブロックサイズが1024であることはとても多いわけではありませんが、1024で遅くなりすぎることはあまりありません。しかし、データ数が1024に近い場合などはもう少し小さいブロックサイズのほうが適切かもしれません。その他のブロックサイズを決める要因として、各ブロックに用意されているブロック内のスレッドでデータを共有するためのメモリであるシェアードメモリのサイズなどがありますがこの記事では触れません。

CUDAプログラムの流れ

CUDAを用いたGPU化プログラムでの一連の流れを説明します。

メモリの確保
データの転送（ホストからGPU）
GPUカーネル関数の実行
データの転送（GPUからホスト）

メモリの確保やデータの転送はカーネル関数の実行に必要なデータサイズさえ分かれば難しくないため、GPUカーネル関数の内容とブロック数、ブロックサイズの決定さえ出来てしまえばあとは簡単です。GPU化したい部分を決め、並列して動くように記述し、データの転送をするというのがCUDAプログラムを作成する際の流れです。これらの流れを踏まえて、実際のプログラム例を確認していきます。

プログラム例（行列積）

今回はN次正方行列の行列積を例にC言語版とGPU版を比較していきます。まずは単純にC言語で作成していきます。関数内で配列を宣言すると、スタック領域に確保され、要素数が大きい場合はエラーがでるので注意です。

...
#define N 2048
...
int a[N][N],b[N][N],c[N][N];

int main(int argc, char *argv[]){  
  init_array(a);
  init_array(b);
  init_array2(c);
  
  eval_start();
  for(int i=0;i<N;i++){
    for(int j=0;j<N;j++){
      for (int k=0; k<N; k++)
        c[i][j] += a[i][k] *b[k][j];
    }
  }
  eval_end();

  output_time();
  return 0; 
}

実行時間

TIME = 50.258955 sec

CPU版では50秒近くかかっています。
次にCUDAプログラムを見ていきます。

...

#define N 2048
#define BS 32

int ha[N*N],hb[N*N],hc[N*N];

__global__ void transpose(int *a, int *b, int *c){
  int i;
  int idx = blockDim.x*blockIdx.x+threadIdx.x;
  int idy = blockDim.y*blockIdx.y+threadIdx.y;
  int id = idy*N + idx;
  c[id]=0;

  for(i=0;i<N;i++)
    c[id] += a[idy*N + i] * b[(i*N) + idx];
}

...

int main(int argc, char *argv[]){
  int *da,*db,*dc;
  
  cudaMalloc(&da,N*N*sizeof(int));
  cudaMalloc(&db,N*N*sizeof(int)); 
  cudaMalloc(&dc,N*N*sizeof(int));
  init_array(ha);
  init_array(hb);

  eval_start();
  cudaMemcpy(da, (int*)ha , N*N*sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(db, (int*)hb , N*N*sizeof(int), cudaMemcpyHostToDevice);
  dim3 dimGrid(N/BS, N/BS);
  dim3 dimBlock(BS, BS);
  transpose<<<dimGrid, dimBlock>>>(da, db, dc);
  cudaMemcpy((int*)hc, dc , N*N*sizeof(int), cudaMemcpyDeviceToHost);
  eval_end();

  output_time();
  cudaFree(da);
  cudaFree(db);
  cudaFree(dc);
  return 0; 
}

実行時間

TIME = 0.082640 sec

実行速度にして約608倍ほどになっています。

CUDA C で記述したソースコードの流れを追っていきましょう。

まず、GPUカーネル関数に記述する処理を考える必要があります。行列積の計算は逐次処理の場合は3重ループで実行され、行方向へ移動するループ、列方向へ移動するループ、1要素を計算するためのループの３つに分かれます。今回は行方向と列方向へのループを並列化し、計算結果の１要素の計算処理を１スレッドが担当することとします（この方針は行列積のGPU化において最適の方針ではありませんが、簡単化のためこのようにしています）。

結果として出力される行列はN×Nの正方行列であるため、必要となるスレッドもN×N個のスレッドです。今回の実装の場合はグリッドやブロックの次元数は１次元でも問題ないですが、２次元を使用しました。そのためグリッドサイズはN/BS×N/BS、１スレッドブロックあたりのスレッド数はBS×BSとしました。１ブロックあたりの最大スレッド数は1024であり、今回１ブロックあたりのスレッド数はBS×BSなのでBS=32としました。N=2048、2048/32=64であるため、64×64個のブロックが生成され、１つのブロックの中に32×32=1024個のスレッドが存在することになります。

1つのスレッドが担当する計算処理は逐次処理版の一番内側のループの計算処理ですが、ブロック分割する場合はインデックスが問題になります。CUDAではビルトイン変数としてブロックサイズやブロック、スレッドのインデックスが存在します。これらのビルトイン変数を用いて配列のインデックスを修正します。今回は行列を1次元配列に格納していますが、2次元配列であっても問題はありません。これでGPUカーネル関数は完成です。

あとはメモリの確保、データの転送をすれば終了です。cudaMallocを用いてN次正方行列を格納するためのメモリを確保した後、cudaMemcpyを用いてホストからGPU側へデータを転送します。グリッドサイズとブロックサイズを指定してGPUカーネル関数を起動した後、cudaFreeを用いて確保したメモリを解放すれば完成です。

ここから更にチューニングする場合はシェアードメモリの利用や、計算処理のさらなる分割などが必要となってきます。他にもワープとブランチダイバージェンスやストリーム処理など意識しなければならないものは多く、様々なCUDAプログラムで完璧なパフォーマンスチューニングをするには更に高度な知識が必要です。

おわりに

今回は行列積を例にGPUプログラミングを行ってみました。コストを考えるとあまり簡単すぎる問題は損かもしれませんが、数分かかるような何回も実行するプログラムがあればGPU化をしてみてはどうでしょうか。眠っているGPUが活躍したら僕は嬉しいです！

この記事を公開する際、同期の@xavifrogさんと社員の@silversさんにレビューしていただきました。ありがとうございました！

参考文献

[1]:What is CUDA?
[2]:OpenAcc
[3]:Cupy
[4]:PyTorch