インテルプロセッサー向け並列化プログラミング
― 対称行列の三重対角化

インテルソフトウェア開発製品による並列化プログラミング (1)

2009/03/02 00:00

ポスト

インテルスレッディング・ビルディング・ブロックを用いた並列化

　インテルコンパイラーの自動並列化機能を用いた並列化では、プログラムの70行目と71行目の二重ループの並列化を行うことができませんでした。72行目でa[k*n+j]の値を書き換える際、k*n+jの値はこの二重ループの実行時に同じ値になることはありませんが、コンパイラーが同じ値にならないことを判断できないため自動並列化が行えません。

　例えば、jがマイナスまたはn以上の値を取り、k*n+jの値がループ内で同じになった場合には、並列化によりプログラムの実行順序が変更されたことで、プログラムの実行結果が実行ごとに違うものになるため、コンパイラーがすべての条件下で並列化が安全に行われると判断ができない場合には自動並列化を行わないのです。

　ここで、インテルで提供しているテンプレート・ライブラリーである、インテルスレッディング・ビルディング・ブロック（以下、インテル TBB）を使って並列化してみます。コンパイラーの自動並列化が不可能だった二重ループ部分を並列化しても、計算結果は正しいことは分かっているので、この部分を並列化することが可能です。また、ここではインテル C++コンパイラー ver11.0の新機能であるラムダ関数についても紹介します。

　インテル TBBを用いて並列化をしたバージョンをサンプルファイルのtridiag-tbb-lambda.ccに示します。

　WindowsでTBBのアプリケーションをコンパイルする場合、デバッグを利用しないときは/MD、デバッグを利用する場合は/MDdオプションを指定しないと、実行時にランタイムライブラリのエラーで異常終了するので、今回は/MDを指定しました。

　また、本プログラムではC++0xの規格に含められたラムダ関数を利用したので、C++0xの規格に対応するための/Qstd:c++0xオプションも指定しました。ラムダ関数のみならず、コード内でC++0xの機能を使用したい場合、このオプションを指定し忘れるとコンパイルエラーとなります。ここで紹介したコンパイラーオプションをまとめると表3になります。

表3　インテル TBBを使うときに有用なコンパイラーオプション

オプション	意味
/MD	マルチスレッドのライブラリを利用する。
/MDd	/MDのデバッグバージョン。
/Qstd:c++0x	C++0xの機能を有効にする。

　ここで、tridiag.ccの70行目の二重ループがどのように並列化されたのかをみてみます。

図4　tridiag.cc 70行目の2重ループ

for (k=i+1; k < n; k++){
    for (j= k; j < n; j++){
        a[k*n+j] += - u[j]* q[k] - q[j]*u[k];
    }
}

　図5にインテル TBBで並列化したループを示します。図4の外側のループを並列化するために、tbb::parallel_forを利用しました。

図5　図4のループをインテル TBBで並列化したループ

tbb::parallel_for(tbb::blocked_range<int>(i+1,n),
    [a,u,q,n](const tbb::blocked_range<int> &r){
        for (int k = r.begin(); k != r.end(); k++){
            for (int j = k; j < n; j++){
                a[k*n+j] += -u[j]*q[k] -q[j]*u[k];
            }
        }
    },tbb::auto_partitioner());

　parallel_forの第一引数として、tbb::blocked_rangeを渡しています。図4の外側ループはint型でi+1からnまでの範囲を回るので、ここではtbb::blocked_range<int>(i+1,n)となります。

　parallel_forの第二引数のBodyとして、C++のラムダ関数を渡しています。ラムダ関数はインテル C++ コンパイラー11.0以降で対応しているもので、1度だけしか使用されない関数の処理を単純化し、STLやインテルTBBのようなスレッド化を容易にするループ構造を使用してテンプレート・ライブラリーを活用するものです。[]がラムダ関数をあらわしています。この中に、外部からキャプチャする変数を書きます。もし何も書かなければ、外部の変数を一切使わないということになります。外部の変数をすべて値渡しでキャプチャするときには「=」、参照渡しでキャプチャする場合には「&」を利用できます。ここでは、必要な変数だけ値渡しでキャプチャしています。

　ラムダ関数の中では、図4の外側ループが、i+1からnまでの範囲を回っていたのに対し、インテルTBBのランタイムライブラリーがループの範囲を適当に分割して、tbb::blocked_rangeとして渡します。渡された範囲に対してループの処理を行えばよいので、r.begin()からr.end()までの範囲をこのループで繰り返すことになります。

　parallel_forの三番目の引数はループをどのように分割するかを指定します。ここでは、tbb::auto_partitionerを用いています。これによって、スレッド間の負荷をモニターしながら自動的に分割します。

　もう一つの例として25行目のループを見てみます（図6）。このループはtmpに結果を連続して足していくループで、reductionと呼ばれるタイプのループです。このループをインテル TBBで並列化する場合にはtbb::parallel_reductionを利用します。この場合は、複数のメソッドを定義する必要があり、ラムダ関数を利用することができないので、クラス（またはストラクト）を定義する必要があります。

図6　trigiag.cc 25行目のループ

tmp = 0.0;
for (j = i+1; j < n; j++) {
    tmp += a[i*n+j]*a[i*n+j];
}

　ここではDdotストラクトを定義しました。ストラクトには分割する場合のコンストラクタ、各範囲に対し処理をするoperator()、結果を合わせるためのメソッドjoinが必要です。

　まず、Ddotストラクトには内部変数として足していった値を保存するためのvalueと配列を参照するためのフィールドa、およびbを定義します。コンストラクタでは、これらのフィールドを適切な値に設定します。特に分割する場合のコンストラクタでは、忘れずに適切な値を設定します。operator()では元のループボディのループインデックスを変更したプログラムとなります。最後にjoinでは、結果を足しこむ必要があります。このストラクトを利用して、並列に処理を行う関数をparallel_ddotとして定義しました。parallel_ddotでは配列a,bと処理する範囲をtbb::blocked_range<int>として受け取り、足しこんだ結果を戻り値とします。この関数を元のループのあった場所から呼び出しますが、このときループの範囲がi*n+i+1からi*n+nとなることに注意が必要です。

図7　インテル TBBで並列化した例

struct Ddot{
    double value;
    double *a;
    double *b;
    Ddot(double* a, double* b) : value(0.0),a(a), b(b){};
    Ddot(Ddot& d, tbb::split) { value = 0;a=d.a; b=d.b;}
    void operator()(const tbb::blocked_range<int>& r){
        double temp = value;
        for (int i = r.begin(); i < r.end(); i++){
            temp += a[i]* b[i];
        }
        value = temp;
    }
    void join(Ddot& rhs) { value += rhs.value;}
};

double parallel_ddot(double *a, double *b, const tbb::blocked_range<int> &r){
    Ddot ddot(a,b);
    tbb::parallel_reduce(r,ddot);
    return ddot.value;
}

tmp = parallel_ddot(a, a, tbb::blocked_range<int>(i*n+i+1, i*n+n,100));

　これを実行した結果、実行時間は13.4秒となりました。

次のページ
インテルマス・カーネル・ライブラリーの利用