CUDA で演算の高速化を体験したく、まず重い演算プログラムを作ってみようと思いまして。
とりあえず、行列積が簡単ですぐに作ることができそうな重い処理かなと思い Visual C++ Express 2008 上で C++ プログラムを書いてみようかと。ちなみ、行列積の復習ということで、以下のサイト様を参考にしました。
先の行列の行と、後の行列の列の各要素をそれぞれ掛けて、足す、というアレです。高校生の時に習いましたね。
その前に、以下のサイト様を参考に、Excel 2007 の MMULT() 関数で 500x500 の行列積を求めてみました。
ま、もちろん数値自体は意味のないダミーな値（というか、コピペで全セル同じような値）。
結果なんですが、quitada は計算速度よりかは（Excel だと遅いのはわかっているので）、計算中の CPU 使用率に着目。当方のマシンはクアッドコアなんですが、計算中の Excel の CPU 使用率は 25% 程度にはりついてました。すなわち、1 コア分しか使ってないわけですね。Excel 2007 からマルチコア対応しており、当方の Excel も 4 スレッドで動作する設定になっているんですが、MMULT() 関数自体はシングルスレッドな関数のようですね。ちょっとがっかり。
ということでついでにマルチコアマシン上で CPU パワーを効率的に使えるよう、Winodws のマルチスレッドプログラミングにも挑戦。
さて、本題ですが、ダミーの行列を作って行列積を求めて計算時間を計測するようなプログラムを書くことになるのですが、まずは行列の C での表現方法の考察。結論からすると、ニ次元配列になるんですが、行と列の要素数を引数で渡して可変にできるようにしようと思ったらちょっとつまづきました。Java 的にやろうとしたらできません！ということで、以下のサイト様を参考にしました。
C/C++でポインタによる多次元配列を連続したメモリ領域に作成する
なるほど、ポインタの登場なんですね。またこの方法を使うと、連続したメモリ領域にデータを確保できるということでアクセス効率が良さそうです。
次に、行列積のマルチスレッド化ですが、基本アイディアは下図の通り。
左側の行列の行ベクトル毎に、担当するスレッドを割り当ててやるわけです。すると、計算結果を配列に書き込む際にもスレッド競合が起きないので、Mutex ロックとかも不要なわけで、マルチコアマシンであればコア数分だけスレッドを起動して分担して計算すれば、CPU パワーを使いきって単純にパフォーマンスの向上が狙えるというわけで。
Windows におけるマルチスレッドプログラミングに関しては、以下のサイト様が参考になりました。
なお、スレッド関数に int 型の引数を渡したかったので以下のように _beginthreadex を使ってスレッドを生成しました。太字のところに着目していただきたいんですが、void * でキャストすると良いようで。
HANDLE *hThread = new HANDLE[threads];
for (int i = 0 ; i < threads ; i++) {
hThread[i] = (HANDLE)_beginthreadex(NULL, 0, threadFunc, (void *)i, 0, NULL);
}
対応する、スレッド関数は以下のイメージ。
unsigned __stdcall threadFunc(void *p)
{
int i = (int) p;
:
なお、上述の for ループでスレッドを起動した後、全スレッドの計算が終了してから次の処理へ移動したい、という作りなんで、それまで待つ必要があります。それに対しては、WaitForMultipleObjects 関数を使いました。こちらに関しては、詳しくは以下のサイト様でお勉強しました。
Microsoft Visual C++ マルチスレッドプログラミング
参考までに、サンプルコードを本ブログに添付しておきます。
MatrixProduct.cpp 
C++ なのにほとんど C 的ですが、Visual C++ 2008 Express とかで新規プロジェクトを作って、Win32 コンソールアプリケーションのテンプレートからこのコードを所定の CPP ファイルに貼りつければコンパイル・実行可能なはず。なお、本プログラムでは実行スレッド数を指定可能です。計算が終了すると、最後に計算が何秒かかったか出力されるんですが、マルチコアマシンでスレッド数を変更して実行すると CPU パワーを使いきって計算時間が高速化される様子が見られます。でも、1 コアでやった場合と 4 コアでやった場合とで 4 倍の差がでるかというとでないですけど。
＃ま、スレッド起動のタイムラグとか、メモリアクセスのレイテンシとか、そもそも計算ロジックがいけてないとか色々理由はあると思いますけど。
ちなみに、各要素の数値は倍精度（double）ではなくてあえて単精度（float）にしております。これは、現行アーキテクチャの NVIDIA の GPU では、倍精度演算は遅いとわかっているので、CUDA に移植した際のパフォーマンス向上インパクトがなさそうだと思ったので。ということで、以下のサイト様にあるインテルアーキテクチャの構造上の理由のためか、コンパイル時にウォーニングがでますが、気にしない、気にしない。
ともあれ、時間のかかる行列積をやろうとすると、数千x数千の行列を指定する必要があり、そうするとメモリを大量に食います。CUDA に移植の際、ホストからデバイスへのデータ転送がボトルネックとなりそうな予感。
