 | レベル: 中級 浅原 明広 (asahara@fixstars.com), Cell/B.E. Developer, 株式会社FIXSTARS 技術開発本部長 / AKD48
2008年 11月 21日 前回の記事では、OpenMPそのものについての解説をするとともに、ベクトルと行列の掛け算をサンプルプログラムとして
"XL C/C++ Alpha Edition Single Source Compiler version 0.9"(XLC SSC)における、スレッド生成にかかるオーバーヘッドを定量的に見積もりました。第三回となるこの記事では、スレッドの数と演算データ量を変えた場合に、どのぐらいの演算処理能力が期待できるのか実際に測定してみましょう。また、後半では第 1 回の記事にも登場したEuler 法シミュレーションの最適化にチャレンジします。
はじめる前に
本記事はIBM SDK 3.0 Extra Package に付属の XLC SSC について述べたもの
です。
XLC SSC はIBM SDK 3.1でα版を卒業し、製品版として登場する見込みです。
スレッド数と高速化効率
こんにちは。すっかり気温が下がって朝方には寒ささえ感じるようになりましたが、ご機嫌いかがでしょうか？「XLCで行こう」の第三回目です。早速ですが、まずは前回の記事に引き続いて"行列とベクトルの掛け算"を題材にして、スレッド数と高速化効率について調べていきたいと思います。
前回の記事のリスト3で掲載したサンプルコードを再度使います。前回同様、行列のサイズを少しずつ変えながら、データ処理時間を計測しますが、ここでSpeed-up Factor を導入しましょう。ここでSpeed-up Factor は、OpenMP非対応コンパイラ(この場合、XLC MA) でコンパイルしてできた逐次処理プログラムの処理時間を、OpenMP対応コンパイラでコンパイルしてできた並列処理プログラムの処理時間で割ったものとします。もし、スレッド分割によるオーバーヘッドが無く、スレッド間の通信も行われない理想的な場合には、Speed-up Factor は、2threadの場合は2.0、4threadの場合は4.0と、スレッド数にリニアに増加していくはずです。結果を図1に示します。縦軸がSpeed-up Factor、横軸Memory Footprint は前回の記事同様、行列となる配列の大きさです。
図 1. スレッド数とSpeed-up Factor
見てわかるように、思ったよりも複雑なグラフプロファイルになっています。また、性能的には、やや残念な結果になってしまいました。4threadまで増やしても、非OpenMPコンパイラによる結果を超えられず、8thread (つまりCell/B.E.丸々一個分)使用してようやく1.5倍の速度向上が得られます。16threadでも、約3倍の速度向上に留まってしまっています。
並列化によるSpeed-up Factorが思ったより伸びなかった理由の一つは、題材としているサンプルコードのデータ量に対する演算量が小さいからと考えられます。今、N行N列の正方行列を入力にしたとすると、このコードでの浮動小数点演算量は、N×N×2 ～ O(N^2) です。一方、Load/Store 命令発行数もN×N + N + N ～ O(N^2)のオーダーです。このようなケースは、XLC SSCに限らずCell/B.E.ではなかなか性能が出にくいことが多いのです。対策としては出来るだけ演算量を大きくすることです。例えば、ベクトル×行列ではなく、行列×行列の場合であれば、浮動小数点演算量は2N×N×N～O(N^3) ですが、Load/Store命令は変わらず N×N + N×N + N×N ～ O(N^2) のオーダーなので、Nの値にもよりますが計算量の方がLoad/Store命令よりもはるかに大きく、Cell/B.E.の性能を引き出しやすくなります。
ともあれ、図1を見てわかることは、Memory Footprint が1MB以上にならないとOpenMPの恩恵には預かれないということです。前回の記事でのオーバーヘッドに関する考察に比べ、OpenMPが効果的に働く閾値がもう一桁上になっていることに注意する必要があるでしょう。
Euler 法シミュレーション最適化
XLC SSCにおける並列化の特性がだんだんわかってきたところで、以前紹介したEuler 法シミュレーションの最適化をやってみましょう。まずはその前に、OpenMPの最適化(一部XLC SSC特有の話も含みます)について、いままでに紹介した話と一部重複するところもありますが、ここにまとめておきます。
不必要なbarrier を避ける
スレッド間の変数の同期をとるためなどbarrierは必須の機能ですが、必要の無い部分にbarrierを入れてしまうと当然ながら先に演算を終えたスレッドが待たされてしまうために並列化効率が落ちます。OpenMPでは、意図せぬbarrierが暗黙のうちに挿入されてしまう場合があります。例えば、#pragma omp for構文は何も指定がなければ対象ループの完了時にbarrierが入ります。また、#pragma omp parallel構文で指定された並列化ブロックを抜ける際にも暗黙のbarrierが入ります。
そこで、nowait節の登場です。nowait節を#pragma omp for構文に付与すると、そのような暗黙のbarrierが行われません。例を挙げましょう。以下は、配列aとbの積とcとdの積を求め、それぞれの和をとるコードです。
リスト1. nowait節を使用したサンプルコード
‥‥
#pragma omp parallel default(none) \
shared(n,a,b,c,d,result) private(i)
{
#pragma omp for nowait
for(i=0;i<n;i++) a[i] *= b[i];
#pragma omp for nowait
for(i=0;i<n;i++) c[i] *= d[i];
#pragma omp barrier
#pragma omp for nowait
for(i=0;i<n;i++) result[i] = a[i] + c[i];
}
‥‥
|
ここでは、a[i] *= b[i];とc[i] *= d[i];の演算に依存関係は無いので、barrierで同期を取る必要はありません。一方、result[i] = a[i] + c[i];の演算ではa[i]とc[i]が新しい値になっていることが必須ですから、直前に#pragma omp barrier構文を明示的に入れています。また、result[i] = a[i] + c[i];の演算にもnowait節が付与されていますが、これはどのみち#pragma omp parallelによる並列化ブロックの完了に伴い、暗黙のbarrierが入ってくるからです。
以上のように、for構文には片っ端からnowait節をつけてしまって、barrierが必要なところは明示的に#pragma omp barrier構文を使って同期をとるようにすると、意図しない同期による処理速度の低下を防ぐことができます。
では、以上を踏まえた上で第一回の記事で取り上げたEuler法シミュレーションに最適化を施してみましょう。
最適化後のソースコードをリスト2に示します。
リスト2. 最適化されたEuler法シミュレーションのソースコード
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <malloc_align.h>
#include <free_align.h>
#ifdef _OPENMP
#include <omp.h>
#else
#define omp_get_num_threads() 0
#define omp_get_max_threads() 0
#define omp_set_num_threads(x) x
#endif
/*
* 4-D vector definition.
*/
typedef struct {
float x, y, z, w;
} vec4D;
struct timeval tv1, tv2;
int elapsed_microsecs;
float elapsed_time;
int main(int argc, char *argv[])
{
int i,m;
int m_max;
float time;
float xpos;
float ypos;
int nthreads;
vec4D *pos; // particle positions
vec4D *vel; // particle velocities
float *inv_mass; // inverse mass of the particles
float dt_inv_mass __attribute__ ((aligned (16)));
// current force being applied to the particles
vec4D force __attribute__ ((aligned (16)));
float dt __attribute__ ((aligned (16))) = 1.0f;
int end_of_time,particles;
if(argc!=4){
printf("usage: >euler <thread num> <end_of_time> <particles>\n");
return 1;
}
nthreads = atoi(argv[1]);
end_of_time = atoi(argv[2]);
particles = atoi(argv[3]);
printf("Th=%d, end_of_time=%d, particles=%d\n",nthreads,end_of_time,particles);
/*
* Memory allocation
*/
pos = (vec4D *)_malloc_align(sizeof(vec4D)*particles,7);
if(pos == NULL){
printf("pos malloc error\n");
return 1;
}
vel = (vec4D *)_malloc_align(sizeof(vec4D)*particles,7);
if(pos == NULL){
printf("vel malloc error\n");
return 1;
}
inv_mass = (float *)_malloc_align(sizeof(float)*particles,7);
if(pos == NULL){
printf("inv_mass malloc error\n");
return 1;
}
//initialize
force.x = 0.0;
force.y = 0.0;
force.z = -9.80; // [m/s]:gravity constant
srand(111);
for (i=0; i<particles; i++) {
pos[i].x = 0.0;
pos[i].y = 0.0;
pos[i].z = 0.0;
vel[i].x = ((float)rand()/(RAND_MAX + 1.0)*2.0 - 1.0);
vel[i].y = ((float)rand()/(RAND_MAX + 1.0)*2.0 - 1.0);
vel[i].z = ((float)rand()/(RAND_MAX + 1.0)*2.0 - 1.0);// [m/s]
inv_mass[i] = 0.001; // [kg]
}
m_max = (int)(end_of_time / dt);
omp_set_num_threads(nthreads);
nthreads = omp_get_max_threads();
gettimeofday(&tv1,NULL);
#pragma omp parallel default(none) \
firstprivate(force,nthreads,inv_mass,dt,particles,end_of_time) \
private(i,time,dt_inv_mass)\
shared(pos,vel)
{
#pragma omp for nowait
// For each particle
for (i=0; i<particles; i++) {
dt_inv_mass = dt * inv_mass[i];
// For each step in time
for(time=0.0;time<end_of_time;time += dt){
// Compute the new position and velocity as acted upon by the force f.
pos[i].x = vel[i].x * dt + pos[i].x;
pos[i].y = vel[i].y * dt + pos[i].y;
pos[i].z = vel[i].z * dt + pos[i].z;
vel[i].x = dt_inv_mass * force.x + vel[i].x;
vel[i].y = dt_inv_mass * force.y + vel[i].y;
vel[i].z = dt_inv_mass * force.z + vel[i].z;
}
}
}
gettimeofday(&tv2,NULL);
elapsed_microsecs =
(int)(tv2.tv_sec - tv1.tv_sec) * 1000000 +
(int)(tv2.tv_usec - tv1.tv_usec);
elapsed_time = (float) elapsed_microsecs * 0.000001;
printf("Nthread = %d, Elapsed time = %f \n",nthreads,elapsed_time);
_free_align(pos);
_free_align(vel);
_free_align(inv_mass);
return (0);
}
|
どうでしょうか。先にあげた最適化のTipに従い、配列のアラインをとり、並列化ブロックにおけるprivate変数を増やし、#pragma omp for構文にnowait節をつけました。また、dt_inv_mass = dt * inv_mass[i];の演算を内部ループの外に出しています。では、気になる演算性能の測定結果を図2に示しましょう。測定は粒子数(particles) 10000個、シミュレーション時間(end_of_time) 300秒で行っています。
図 2. スレッド数とSpeed-up ratio
第一回で掲載したEuler法のコードと比較すると、どのスレッド数であっても100倍程度の演算速度向上が得られています。ちなみに、変更箇所で最も高速化に寄与したのは配列のアラインでした。図には、参考としてCell SDK 3.0に含まれている最適化されたEuler法のコードによる演算速度も掲載しています。ここでSDK版の演算結果は、ppu-gccとspu-gccによる一般的なCell/B.E.プログラミングコードで、SIMD化、ループアンローリング等、Cell/B.E.ではおなじみの最適化手法を適用済みのものです。スレッド数が8個の場合の値のみ掲載しています。
さすがに、手動で丹念に最適化されたSDKコードとはまだ2倍程度の差がありますが、XLC SSCもまだversion 0.9のベータ版ですから、今後のさらなる改善に期待しましょう。ただし、開発工数の面から考えると、現時点でもXLC SSCに大きなアドバンテージがあります。一つの目安としてコードの行数をとりあげると、XLC SSC版のEuler法シミュレーションコードは、値の初期化部分を除けばたかだか130行程度なのに対し、Cell SDK版はppuコード、spuコード合わせて300行近くあるのですから。
Speed-up Factor も見てみましょう。こちらは少し興味深い結果が得られています。図3をご覧下さい。
図 3. Euler法シミュレーションプログラムのSpeed-up Factor
図3ではスレッド数は16に固定、x軸に粒子数(particles)、y軸にSpeed-up Factor (XLC MAで生成された逐次処理プログラムの処理速度との比)を取っています。3本のグラフは、シミュレーション時間(end_of_time) を100秒、1000秒、10000秒と変化させたものです。前節のベクトル-行列の掛け算と異なり、Euler法シミュレーションでは、16スレッドの際に6倍のSpeed-up Factor を得ることが出来ています。
グラフの粒子数8000付近でSpeed-up Factor が急激に高くなる不連続点が存在しています。これはXLC SSCのコードがこの付近で急激に改善されているというわけではなくて、逆にXLC MAで生成されたシングルスレッドプログラムがparticle数8000付近で急に遅くなっているようです。この原因のはっきりしたところはよくわかりませんが、PPEのキャッシュ容量とループで処理するデータサイズが何らかの要因になっているかもしれません。ちなみに、このグラフで粒子数をさらに増やして200000個あたりまで測定を行ってみましたが、Speed-up Factorはだいたい6.5前後で、特に不連続点は見つかりませんでした。また、シミュレーション時間(end_of_time)の方も100000sec辺りまで増加させて測定してみましたが、end_of_time が10000sec以上ではグラフのプロファイルにはほとんど変化が見られませんでした。
まとめ
今回は前回に引き続きベクトルと行列の掛け算のサンプルコードと、Euler法シミュレーションのコードを使って処理速度がデータサイズやスレッド数に伴ってどのように変化するかを見てみました。
また、Euler法シミュレーションでは、適切な最適化を施すことで、前々回の記事のコードと比べると10倍以上の速度改善が見られました。
Cell SDKで提案されているプログラミングモデルには、Accelerated Library Framework (ALF)もあります。こちらはMultiple program multiple data (MPMD) プログラミングモジュールで構成されるフレームワークで、やはり開発工数を削減することができます。詳しくはdeveloper Works内にあるALFに関する記事 を読んでいただきたいのですが、これら3つのCell/B.E.プログラミングモデル、1. 通常のppu用spu用2つのコンパイラによるコーディング、2. ALFによるコーディング、3. XLC SSCによるコーディング、それぞれの立ち位置は大まかには図4のようなイメージになります。つまり、XLC SSCでは、パフォーマンスではデュアルソースコンパイラによる手動の最適化に劣るものの、開発工数に関しては大幅な削減を期待することができるでしょう。
今後製品版が登場すればパフォーマンスに関してもさらなる改善がなされていると思います。この記事を参考にして、お手軽に並列化処理を行うことができるXLC SSCを、是非試してみてください。
図 4. 3つのCell/B.E.プログラミングモデルの違い
参考文献 学ぶために
製品や技術を入手するために
議論するために
著者について | 
| | 浅原明広(Akihro Asahara, Ph.D.) は、Cell/B.E. 関連の開発業務に携わっているエンジニアです。私事ですが、2008年9月30日に日本IBM を卒業し、株式会社FIXSTARS に移りました。会社は変わりましたが、今後もCell/B.E.に関わる仕事を続けていくことに変わりはありません。またこのdeveloper Works での執筆活動も今までどおり続けていく予定ですので、これからもどうぞよろしくお願いします。左の写真はわかりにくいですが、海に潜るの図です。その心は、8SPEと1PPEのハチ・ワンにダイブするという・・・ |
記事の評価
| |
