インテルTBBライブラリを用いた並列処理の実装

粒子法シミュレーションに本格的にハマっていくうちに、計算時間が気になりだしたので、せっかくCore2Duoのマシンなんだし（Quadでないのが悔やまれた！）、マルチスレッド化でもしてみようかと、インテルのTBBライブラリの、parallel_for機能を使ってみた。粒子法では、同時刻における個々の粒子の挙動は全て互いに独立しているので、並列化にはもってこいの題材でもある。それに、TBBは今やオープンソースだ。
ここから誰でもダウンロード取得できる。
並列化というと何やら面倒そうに聞こえるが、TBBにおいては、その考え方はstlの
for_eachを呼ぶのとほとんど同じで、コーディングする上でそれ以上のことを考える必要は一切ない。要はfor_eachのTBB版を呼べば、それが自動的にマルチスレッドで並列して走ってくれるのだ。
parallel_forは簡単に、stl::for_each のTBB版と考えることができる。stl::for_eachと同じく関数テンプレートとして実装されており、ループをまわしたい範囲の開始点と終了点と、各要素について何を呼び出すかを指定するためのFunction述語を受け取る点も、考え方としては全く同じだ。
まずはよく使い慣れた、for_eachのインタフェースがこれ。

// stl::for_each
template< class IterType, class Function >
Function for_each( IterType first, IterType last, Function func ){...}

そして、以下が、parallel_forのインタフェースだ。

// tbb::parallel_for
template<typename Range, typename Body>
void parallel_for( const Range& range, const Body& body ) {...}

stl::for_eachと大きく違う点は、stl::for_eachにおける一番最後のFunctionの引数が、ループで回る要素そのものを一つ受け取るUnaryFunctionではなく、Rangeという、要素のコンテナ的な型を受け取るようになっている点だ。
（これを見ればわかるように、for_each におけるFunction にあたるものは、
parallel_for では、Body という型名になっている。）
parallel_for のテンプレート引数型としての、BodyクラスやRangeクラスが満たさなければならない、コンセプトは、TBB conceptsに明記されている。（これは一式をダウンロードしたときについてくるドキュメントの中でも確認できる。）
しかし、心配は無用だ。
Rangeの標準実装としてのblocked_rangeテンプレート型が、TBB名前空間に標準で定義されている。そして、Bodyは、要するに、コピーコンストラクタとデストラクタと、Range& を受け取る関数呼び出し演算子operator()を実装していればいいのだ。それだけだ。
それでは実際の使い方を見ていくことにしよう。
VC9にて空のプロジェクトを作成し、以下をビルド実行すれば、ネイティブループと、TBBを使った場合とのパフォーマンスの違いが体感できるはずだ。
（勿論tbbヘッダ一式への追加のインクルードパス設定やtbb.libへのリンクパスの設定、実行時のtbb.dll へのパスが通っていることが必要だ。）

#include <iostream>
#include <time.h>
#include <vector>
#include <cmath>
using namespace std;

//TBB関連組み込みファイル
#include "tbb/task_scheduler_init.h"
#include "tbb/blocked_range.h"
#include "tbb/parallel_for.h"

#ifdef _DEBUG
#pragma comment( lib, "tbb_debug.lib" )
#else
#pragma comment( lib, "tbb.lib" )
#endif

class Node
{
public:
  Node( vector< float >* results ){
	m_results = results;
  }
  ~Node( void ){}
  void operator()( const tbb::blocked_range< int >& range ) const{
    for( int i=range.begin(); i!=range.end(); ++i ){
      m_results->operator[](i) = pow( (double)i, 5 );
    }// i
  }
private:
  vector< float >* m_results;
};

int main( int artc, char* argv[] )
{
  const int num = 1.0e8;

  vector< float > results;
  {
    results.resize( num );
  }

  // Native Solution
  clock_t time0 = clock();
  {
    for( int i=0; i<num; ++i ){
      results[i] = pow( (double)i, 5 );
    }// i
  }
  clock_t time1 = clock();
  //

  //  TBB Solution
  clock_t time2 = clock();
  {
    tbb::task_scheduler_init tbb;
    // parallel_for呼び出し
    tbb::parallel_for( tbb::blocked_range<int>( 0, num, 1000 ), Node( &results ) );
    // ↑0番目からnum番目までの範囲を、1000要素単位で重複なく区切り、
    // その区切られた範囲のそれぞれを各スレッドにて、
    //  Node型関数オブジェクトのoperator() の引数として呼ばせる。
    tbb.terminate();
  }
  clock_t time3 = clock();
  //
  cout << "native solution takes "
        << (time1-time0)/(float)CLOCKS_PER_SEC <<" sec"<<endl;
  cout << "TBB solution takes "
        << (time3-time2)/(float)CLOCKS_PER_SEC <<" sec"<<endl;
  return 0;
}

以上のサンプルでは1000要素を区切り単位としているが、ループの中の処理によって、どのような値を区切り単位にするかで、パフォーマンスは結構違ってくるようだ。あまり細かく区切りすぎてもいけないようだ。
ちなみに筆者の環境では、コアが2個しかないのだが、区切り数1000で、パフォーマンスは大体2倍になった。

インテル スレッディング・ビルディング・ブロック ―マルチコア時代のC++並列プログラミング