はじめに
「MATLABでGPUコンピューティング」
の記事です。
担当は@cyan_rej55です。
どんな内容？
GPU計算を利用するためにはParallel Computing Toolboxが必要となります。
また、Compute Capability 2.0以上を満たすCUDAに対応したnvidia社のGPUが必要です。
こちらで対応しているかどうかを確認できます。
なお、本記事では実行環境としてMATLAB R2015bを利用します。
GPUコンピューティングとは？
今回のテーマであるGPUコンピューティングとは、GPU(Graphics Processing Units)を利用した計算のことを指します。 GPUはその名の通りグラフィックス周りの演算に特化したプロセッサユニットです。グラフィックスの演算では単純な演算を大量に行うことが多く、GPUではそれに特化した形でハードウェアが構成されています。
GPUはグラフィックスの処理の用途で用いられていたのですが、大量の演算を一度に行う並列演算の能力を他の分野に応用するGPGPU(General-Purpose Computing on GPU)が近年流行っています。 応用分野としては、画像処理などのグラフィックス処理はもちろんとして、物理シミュレーションや機械学習など多岐に渡ります。
これらを扱うにはGPUでの処理を記述するプログラムを作る必要がありますが、多くの場合はCUDAを利用する必要があります。しかし、CUDAによる実装ではGPU上のメモリ管理などをすべて自分で行う必要があるなど、非常に手間がかかります。そこで今回はMATLABを利用して手軽にGPUコンピューティングに触れてみよう！というのが今回の趣旨となります。
準備：自分のGPUデバイスを調べよう
MATLAB上で自分の使用しているGPUデバイスの情報を調べるには、gpuDeviceを利用します。MATLABのコマンドウィンドウ上で
>> gpuDevice
と入力すると自分の環境を調べることができます。複数GPU環境を利用されている場合は、gpuDevice(index)という形でインデックスを引数で与えることで利用するGPUを選択して調べられます。 今回は話が複雑になるのでシングルGPU環境という前提で話を進めます。
基本編：gpuArrayを使おう
GPU上で計算を行う場合には、GPUのメモリ上にデータを転送する必要があります。そのための関数として gpuArrayがあります。また、逆にGPU上のメモリを拾ってくるときにはgatherという関数を利用します。
GPU上での計算の基本的な流れは、
- gpuArrayでGPUにデータを転送
- 演算を実行
- gatherでGPU上のデータをワークスペースに転送
それでは、実際に動かしてみましょう。今回は固有値計算を例にしてみたいと思います。
こちらのスクリプトを実行してみます。私の環境では以下のようになりました。
>> gputest
経過時間は 14.211752 秒です。
経過時間は 6.273133 秒です。
このように、高速化されていることがわかるかと思います。先ほどのコードのNの値を増やせば増やすほどその効果を実感できると思います。ただし、Nが少ない場合には逆にGPUのほうが遅くなることも多々あります。この点には注意してください。
このように、gpuArrayで確保したデータに対して通常の関数を適用するだけでGPUの恩恵を受けることができます。gpuArrayのデータに対応した関数一覧はこちらに掲載されています。
応用編：arrayfunを使って並列計算
ここまで、gpuArrayとビルトイン関数を利用して簡単にGPUを利用する方法を説明しました。 しかし、ビルトイン関数をそのまま使えないなど、自分で関数を作りたいという要求もあるかと思います。 そこで、arrayfunという関数を利用して並列計算を行ってみたいと思います。
[B1, B2, ..., Bm] = arrayfun(fun, A1, A2, ..., An)は、n個の引数からm個の値を返す関数funを、A1, A2, ... , Anの各要素に適用してB1, B2, ..., Bmに返す関数です。 反復法的に書くと、
for i = 1:k
[B1(i), B2(i), ..., Bm(i)] = fun(A1(i), A2(i), ..., An(i))
end
という処理をまとめてやってくれるものになります。今回は例として、行列の各要素に値が0.5を超えるものは1に、超えないものは0にするという処理をやってみたいと思います。
（注意：この処理の場合はarrayfunを使うよりもずっと簡単なやり方がありますが、今回はあくまでも例ということでarrayfunで実行します。）
上記2つのスクリプトファイルを用意した上で実行してみます。
>> gputest2
経過時間は 15.774822 秒です。
経過時間は 2.594181 秒です。
やはり高速化されていますね。これで個々の配列要素に対してまとめてGPUで演算を行うという処理ができるようになりました。
実践編：MATLABで2次元移流拡散方程式を解く
今回は、GPUの利用例のひとつとして2次元の移流拡散方程式を数値的に計算するプログラムを作ってみようと思います。なお、2次元移流拡散方程式を解く際にはステンシル計算というものを行う必要がありますが、上記のarrayfunでステンシル計算を行う場合はちょっとしたテクニックが必要です。それについては話がややこしくなるので詳細は私の記事を参考にしてください。
移流拡散方程式は、 $$ \frac{\partial \phi}{\partial t} + \nabla \cdot (c\phi) = D\nabla^2 \phi $$ で表されます。この方程式は、物理量$\phi$が速度$c$で移流し、拡散係数（ここでは簡単のため定数）$D$のもとで拡散する様子を表すものです。特に今回は速度場が定常状態であると考えて、$c$の$x$方向成分を定数$u$、$y$方向成分を定数$v$とすると、 $$ \frac{\partial \phi}{\partial t} + u\frac{\partial \phi}{\partial x} + v\frac{\partial \phi}{\partial y} = D\left( \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} \right) $$ のように書き表すことも可能です。
今回は、有限差分法を用いて2次元の移流拡散方程式を解きます。離散化手法としては、移流項$\nabla \cdot (c\phi)$を1次の風上差分法、拡散項$D\nabla^2 \phi$を中心差分法で離散化します。また、時間発展はオイラー法で行います。こちらの離散化の詳細については長くなるので省略します。 そして、境界条件は周期境界条件を仮定します。
ステンシル計算では、配列の各要素に対して、その近傍の要素の値を利用して配列のデータをアップデートします。先ほどのarrayfunを利用する場合、近傍の要素の値の参照をうまく行うために入れ子関数を利用します。入れ子関数は関数の中に関数を定義する手法で、親関数で定義された変数を利用することが可能になります。ということで、移流拡散方程式を解く関数の中に、ステンシル計算を行う関数と境界条件処理を行う関数を定義するという形で実装します。
では、こちらを実行してみましょう。 パラメータ等を指定したら、adv_dif_2d_runを実行すれば計算が行われ、結果がアニメーションで表示されます。 私の実行環境では以下の画像のような結果になりました。
このように、GPUを使うと大幅に高速化されていることがわかります。
まとめ
今回は、MATLABでGPUコンピューティングを行う方法をご紹介しました。いかがでしたでしょうか。
「GPUを使って高速化したい！でもCUDAを使うほど面倒なことはしたくない！」という方には非常にオススメです。基本的なビルトイン関数ではGPUに対応しているので、画像処理やフーリエ変換など基本的な処理なら基本編の内容だけで簡単にこなせます。少しマニアックなことをやりたい場合には応用編の内容が必要になります。
ただし、極限までチューニングをして高速化を極めたい場合にはCUDAを使わざるを得ないかと思います。CUDAを使う場合は結構労力をがいるので、どこまでやりたいかを考えて、そのときどきで適切なものを使うのがよろしいかと思います。
私の記事が皆様のお役に立てれば幸いです。
それでは、よいMATLABライフを！
