OpenMPのすすめ - あさっては曇りのち晴れ

いよいよ冬本番、大学の近くの店でもメニューにカキフライが加わっていました。かなり久しぶりに食べましたが、やはり美味ですね。
さて、最近気分がいまいちぱっとしないのですが、ぱっとしないついでに、Jacobi法にOpenMPで並列化をかましてみました。一般的にはGaussSeidelの方が収束が速いと言われていますが、Jacobi法の方は簡単に並列化できます。計算の大部分が並列化できるので、Core 2 Duoに2スレッド計算させると、ちゃんと倍近い性能が出ます。すばらしい。
このOpenMP、数値計算をかじった方ならご存じの通り、とってもお手軽です。一番簡単なのは

#pragma omp parallel for //OpenMPのおまじない
for(int i=0; i<1000; ++i){
	vector[i] = 0.0;
}

みたいなのでしょうか。omp.hをincludeして、gcc4.2以降なら-fopenmpオプションをつけてコンパイルしてあげればできあがりです。最近のLinuxとかだとデフォルトで使えるかと思います。これだけじゃああまりに芸がないですが、

double sum = 0.0
//sumに足し込むときにみんなで一斉にやると混乱するので、おまじないを追加
#pragma omp parallel for reduction(+:sum)
for(int i=0; i<1000; ++i){
	sum += vector[i];
}

とかやれば使えそうな雰囲気に。あとは、これは無意味な例ですが、

int temp;
//privateにしないと、例えばi=0の担当スレッドがtempを0にし、
//i=1の担当スレッドがその0をvector[1]に入れてしまったりということがおきる
#pragma omp parallel for private(temp)
for(int i=0; i<1000; ++i){
	temp = i;
	vector[i] = temp;
}

あたりを組み合わせれば、これだけで結構使い物になります。
まあ、コンパイラまかせの自動並列化が一番安全で楽ですが、ひとつ内側のループを並列化したいときなど、少し考えて、人間側でおまじないを追加してやるといい結果が出ることもあります。ちょっとした計算をさせているのだけれど遅い、という方、検討してみてはいかがでしょうか。
さらに並列化ついでに、個人的に特に意味もなくあこがれのあったSSEに挑戦。どうやらdoubleだと2個いっぺんに計算できるらしいとか聞きつけて、intrinsicなんたらを書くとアセンブラを書かなくてもいいらしいうえコンパイラがよしなに最適化してくれるらしい、とやってみましたが、見事撃沈。実行時間はほとんど変わらずでした。キャッシュなども意識しつつ、きっちりと考えて使わないと性能は出ないようです。結論、こっちは素人が手を出すべきではないらしい・・・