OpenMPのリダクション演算を自分でする方法

質問

OpenMPでプログラムを並列化しているのですが、以下の内積を並列化する際にreduction演算をしなくてはなりません

for( i=0; i<N; i++ ){
　r += x[i] * y[i];
}

↓

#pragma omp for reduction(+:r)
for( i=0; i<N; i++ ){
　r += x[i] * y[i];
}

このように、for文に対してreduction節を入れてスレッド並列化するのが普通ですが、ここでBLASというライブラリを使いたいので、reductionを自分でしたいと考えています

int INC = 1;
int n = i_end - i_start;
for( i=i_start; i<i_end; i++ ){
　r += ddot_( &n, &(x[i_start]), &INC, &(y[i_start]), &INC );
}
/* reduction here */

0～N-1までの配列の要素を、自分でスレッドに分けて、i_start～i_endまでの要素の内積を各スレッドがddot_関数を使って計算します
その後、reductionの計算をする必要があるのですが、一体どうすればreduction演算ができるのか分かりません
どなたか分かる方、ご教授ください

なお、スレッド並列版のBLASは使わない予定です
あくまで、自分でスレッド並列化を行ない、その中で逐次版BLASを呼び出すという方向で考えています

回答No.2

ki073

No.1のお礼欄について

threadが32有ったときに、順番に足すと31回加算をしないといけない。隣どおしの加算でも16+8+4+2+1=31と回数は同じだと思いますが、ちょっと眠いので頭が回らんので勘違いをしているかも。確かに並列化の効果は出るかも知れませんが、現実的にはどうなのでしょうか。

それはともかく、この場合だけを考えると計算の負荷は分担する配列要素数に比例するので、thread数をあまり増やさずにCore数と同じにしれば一番効率が良いように思いますがいかがでしょうか。この計算はコンパイラの最適化によるベクトル化が働きますしprefetchも効きやすいので意外と2番目の方法が健闘しそうな気がしますが、

補足

隣同士の加算は同時に行えるので、
32threadsで
16 + 8 + 4 + 2 + 1
ではなく、threadあたり
1 + 1 + 1 + 1 + 1 = 5
の5回の計算になります

もちろん、コア数と同じにはしますが、反復法の中で行う処理なので、多数回呼び出されることからも、こういう細かい高速化も考慮している次第です

あと、やっぱりBLASは必須です
BLAS使うと全然速度違いますし

回答No.1

ki073

reductionで問題ないように見えますが、どうしても使いたくないのならatomicかcriticalを使えば計算できます。
http://www2.cc.kyushu-u.ac.jp/scp/system/library/OpenMP/openmp0209.pdf
の21ページからを参考にしてください。
r += ddot_( &n, &(x[i_start]), &INC, &(y[i_start]), &INC );
で全体をcriticalにしたら、全然並列処理になりませんので、一旦private変数に代入して、加算だけをcriticalにしてください。

3つとも同じ答えがでるつもりで書いているのでしょうか？最後のはこのままでも答えが合わないような気がしますので、考えてみてください。

お礼

最後のはrがprivate変数だと思ってください
その後にreductionでprivate変数のrから総和を求める感じです
一応、現段階ではmasterで総和を取ってるのですが、masterにしろcriticalにしろ、その部分では処理時間のオーダーがスレッド数になってしまいます
各スレッドが隣同士と足しあわせながらやればオーダーがlogになるので、できればそういう方法がないかなぁと思い質問させていただいた次第です
少なくとも、MPIのReduce通信はそういう風にしてますし、OpenMPでもそういう風にできないかと