 | レベル: 中級 浅原 明広 (asahara@fixstars.com), 株式会社FIXSTARS 技術開発本部長 / AKD48
2009年 03月 23日 前回の記事に引き続き、Cell Broadband Engine (Cell/B.E.) 用 並列化フレームワーク "Cell Superscalar (CellSs)" を取り上げます。
今回は、CellSsの最大の特徴であるIntrinsics との併用を行い、Cell/B.E.のパフォーマンスを最大限まで引き出してみます。
はじめに
前回の記事で、Cell Super Scalar (以下CellSs) は、同じ自動並列化コンパイラであるIBM XLC SSCと比べていくつもの利点があることを説明しました。その中でも「パフォーマンス」ということに絞った場合に、最大の利点となってくるのが「SPU組み込み関数(Intrinsics)」を利用できる、という点です。IBM XLCでは、残念ながらIntrinsics を利用することができないので、Cell/B.E.の特徴である高速なSingle Instruction Multiple Data (SIMD) 演算を行うことができず、性能を引き出すことはできません。
今回の記事では、前回CellSsで並列化したEuler粒子シミュレーションをさらにintrinsicを使ってSIMD化し、どのぐらいの性能が出せるか検証してみます。
Intrinsics
前回の記事で掲載した"リスト1"のコードをIntrinsics を使ってSIMD化します。SIMDの詳細はCell Programming Tutorial 等に譲りますが、同じデータ型を持つ複数の変数を一つの"ベクトル"にパックし、四則演算や論理演算等をベクトルの単位で一度に行う処理方法です。
SIMD化の際のデータの取り方は、この問題の場合、構造体vec4Dをそのまま一つのベクトルとみなし、一つのvector floatに粒子のx,y,z,w の4要素が入るArray Of Structure (AOS) フォームと、粒子数の要素を持つxの配列、yの配列、zの配列を用意し、例えば (x0, x1, x2, x3) と、粒子のx座標のみをvector float にパックしたStructure Of Array (SOA) フォームとが考えられますが、計算効率を考えると、SOAフォームが最適です。ただしSOAフォームを取る場合は、粒子数が4で割り切れない場合は、余りの処理をループの外に別途書く必要があります。今回の記事ではそのあたりはサボっておりますので、注意が必要です。
では、早速コードをみてみましょう。
リスト 1. Euler法のサンプルプログラム (CellSsによる並列化+SIMD化)
001: #include <stdio.h>
002: #include <stdlib.h>
003: #include <sys/time.h>
004: #include <spu_intrinsics.h>
005:
006: struct timeval tv1, tv2;
007: int elapsed_microsecs;
008: float elapsed_time;
009:
010: #define END_OF_TIME 300
011: #define PARTICLES 8192
012: #define BLOCK_SIZE 512
013:
014: // 4-D vector definition.
015: typedef struct {
016: float x, y, z, w;
017: } vec4D;
018:
019: // SOA
020: float pos_x[PARTICLES], pos_y[PARTICLES], pos_z[PARTICLES];
021: float vel_x[PARTICLES], vel_y[PARTICLES], vel_z[PARTICLES];
022: vec4D force;
023: float inv_mass[PARTICLES];
024: float dt = 1.0f;
025:
026: #pragma css task inout(pos_x,pos_y,pos_z,vel_x,vel_y,vel_z) input(inv_mass,force,dt)
027: void task_euler_simd
028: (float pos_x[BLOCK_SIZE], float pos_y[BLOCK_SIZE], float pos_z[BLOCK_SIZE],
029: float vel_x[BLOCK_SIZE], float vel_y[BLOCK_SIZE], float vel_z[BLOCK_SIZE],
030: float inv_mass[BLOCK_SIZE], vec4D force, float dt)
031: {
032: float dt_inv_mass __attribute__ ((aligned (16)));
033: vector float* pos_x_v;
034: vector float* pos_y_v;
035: vector float* pos_z_v;
036: vector float* vel_x_v;
037: vector float* vel_y_v;
038: vector float* vel_z_v;
039:
040: vector float force_x_v, force_y_v, force_z_v;
041: vector float dt_v, dt_inv_mass_v;
042:
043: pos_x_v = (vector float *)pos_x;
044: pos_y_v = (vector float *)pos_y;
045: pos_z_v = (vector float *)pos_z;
046:
047: vel_x_v = (vector float *)vel_x;
048: vel_y_v = (vector float *)vel_y;
049: vel_z_v = (vector float *)vel_z;
050:
051: force_x_v = spu_splats(force.x);
052: force_y_v = spu_splats(force.y);
053: force_z_v = spu_splats(force.z);
054:
055: dt_v = spu_splats(dt);
056:
057: // For each particle
058: for (int i=0; i<BLOCK_SIZE/4; i++) {
059: float dt_inv_mass = dt * inv_mass[i];
060: dt_inv_mass_v = spu_splats(dt_inv_mass);
061:
062: // For each step in time
063: for(float time=0.0;time<END_OF_TIME;time += dt){
064: // Compute the new position and velocity as acted upon by the force f.
065: pos_x_v[i] = spu_madd(vel_x_v[i], dt_v, pos_x_v[i]);
066: pos_y_v[i] = spu_madd(vel_y_v[i], dt_v, pos_y_v[i]);
067: pos_z_v[i] = spu_madd(vel_z_v[i], dt_v, pos_z_v[i]);
068:
069: vel_x_v[i] = spu_madd(dt_inv_mass_v, force_x_v, vel_x_v[i]);
070: vel_y_v[i] = spu_madd(dt_inv_mass_v, force_y_v, vel_y_v[i]);
071: vel_z_v[i] = spu_madd(dt_inv_mass_v, force_z_v, vel_z_v[i]);
072: }
073: }
074: }
075:
076: int main(int argc, char *argv[])
077: {
078: //initialize
079: force.x = 0.0;
080: force.y = 0.0;
081: force.z = -9.80; // [m/s]:gravity constant
082:
083: srand(111);
084:
085: for (int i=0; i<PARTICLES; i++) {
086: pos_x[i] = 0.0;
087: pos_y[i] = 0.0;
088: pos_z[i] = 0.0;
089:
090: vel_x[i] = ((float)rand()/(RAND_MAX + 1.0)*2.0 - 1.0);
091: vel_y[i] = ((float)rand()/(RAND_MAX + 1.0)*2.0 - 1.0);
092: vel_z[i] = ((float)rand()/(RAND_MAX + 1.0)*2.0 - 1.0);// [m/s]
093:
094: inv_mass[i] = 0.001; // [kg]
095: }
096:
097: gettimeofday(&tv1,NULL);
098: #pragma css start
099: for (int i=0; i<PARTICLES; i+= BLOCK_SIZE) {
100:
101: task_euler_simd(pos_x +i, pos_y +i, pos_z +i,
102: vel_x +i, vel_y +i, vel_z +i,
103: inv_mass+i, force, dt);
104: }
105: #pragma css finish
106:
107: gettimeofday(&tv2,NULL);
108:
109: elapsed_microsecs =
110: (int)(tv2.tv_sec - tv1.tv_sec) * 1000000 +
111: (int)(tv2.tv_usec - tv1.tv_usec);
112: elapsed_time = (float) elapsed_microsecs * 0.000001;
113:
114: printf("Elapsed Time = 0.000000 \n",elapsed_time);
115:
116: return (0);
117: }
|
前回のコードは86行だったので、SIMD化に伴い1.5倍ほど行数が増加しています。それでもPPE, SPE両方のコードを用意し、データのやり取りやスレッドの管理も書かなければならない、libspe2を用いたコードに比べればはるかに少ない行数です。では、行順にコードの解説をしていきましょう。
-
4行目:
#include <spu_intrinsics.h>
SPU Intrinsicsを使用するために必要です。
-
20行目:
float pos_x[PARTICLES], pos_y[PARTICLES], pos_z[PARTICLES];
この演算ではSOAフォームの方がSIMDに乗りやすいため、x,y,zのそれぞれの配列を個別に確保します。velocity の配列も同様にします。また、本来Cell/B.E.プログラミングではこの配列についてはアラインを明示的にとるべきですが、CellSsではアラインをとってもとらなくても速度に変化はありませんでした。おそらくSource-to-Sourceコンパイラでよしなにやってくれているのでしょう。
-
26行目:
#pragma css task inout(pos_x,pos_y,pos_z,vel_x,vel_y,vel_z) input(inv_mass,force,dt)
SIMD版のEuler粒子計算のカーネルになります。SIMD演算に慣れている方には問題ないでしょう。初めての方でも、spu_splats(引数のスカラー数を全ての要素に持ったベクトルを返す), spu_madd(第一引数に第二引数を掛け算し、さらに第三引数を足し算する)の2つしかIntrinsics が出てこないので、わかりやすいと思います。前回のコードの同一部分と見比べてみてください。
-
58行目:
for (int i=0; i<BLOCK_SIZE/4; i++) {
SIMD演算で4要素同時に演算しているので、ループの数は1/4になります。注意してください。
76行目以下のmain関数の中身はほとんど変更がなく、問題ないでしょう。さて、どのぐらいパフォーマンスが上がったか見てみましょう。測定条件は以下のとおりです。
- 粒子数: 8192 個
- Step時間: 300 sec
- ハードウエア: PLAYSTATION 3
- SPE数: 4個
図1に測定結果を示します。前回の記事で紹介したSIMD無しのCellSsコードと、IBM SDKに収録されている手動で最適化されたコードと比較しました。見てわかるように、SDK最適化コードより、さらに演算時間が10%程短縮短縮されています。手動で最適化されたコードには敵わないとおもっていたので、意外でした。理由を調べてみると、IBM SDKのコードは粒子のデータ構造がAOSのまま変更できないと仮定し、AOSのまま最適化を行っていることがあげられます。ともかく、CellSsのオーバーヘッドは大変少なく、工夫次第でlibspe2による手作業の並列化に勝るコードを書くことができることがわかりました。
図 1. SIMD化したCellSsコードとIBM SDKとの比較
次に、SPEの数を変化させた場合の処理時間の変化を見てみます。ただし、粒子数8192, Step時間300 の条件では処理が軽すぎて、SPEの数を変えてもあまり処理時間に変化は見られなかったので、少し負荷を増やし、
- 粒子数: 1048576 個
- Step時間: 1024 sec
- ハードウエア: IBM QS22
とします。
今回はIBM BladeCenter QS22 で検証を行っているので、SPEは最大16個使うことができます。
使用SPEの数は、環境変数CSS_NUM_SPUSで決定されます。例えばSPEの数を6個使用したい場合は
とコマンドを打ちます。特に環境変数を設定していない場合は、デフォルトとしてSPEが8個使用されます。さて、この環境変数CSS_NUM_SPUSを1から16まで変えながら処理時間を測定した結果を図2に示します。
ここで、Speed up factor とは、SPEが1個の場合の処理時間を1とした場合の処理速度向上の度合いで、具体的にはSPEが1個の場合の処理速度をSPEを変化させた場合の処理速度で割った値です。SPEの数の増加にしたがって、綺麗に処理時間が減少しているのがわかります。Euler粒子シミュレーションプログラムは、粒子間の相互作用がなく、完全な並列演算ができるため、もともと計算資源の数にしたがってリニアにパフォーマンスが出やすいのですが、この結果はCellSsによる各種オーバーヘッドが非常に少ないことを示唆しています。
図 2. SPE数を変化させた場合の実行時間の変化
まとめ
今回の記事では、Intrinsics を使用し、Euler粒子シミュレーションプログラムを高速化することに成功しました。自動並列化ソリューションは、性能面がいま一つで、実際は実用的でないものが多い印象がありますが、CellSsでは、SPEコードの中身を制約なくいじることができるので、手作業の並列化コードに近い性能を引き出すことができます。実際、上記コードをループアンローリング、ソフトウエアパイプライニング、等の手法を用いて、さらに高速化することも可能だと思います。是非、チャレンジしてみてください。
参考文献 学ぶために
製品や技術を入手するために
議論するために
著者について | 
| | 浅原明広は、Cell/B.E.の発表当初から関連の開発業務に携わってきました。多数のハードウエア、多数のプログラミングモデルが混在し、マルチコアプロセッサ業界は今が旬。市況はさんざんですが、エンジニアにとっては今年は非常に面白い年になるのでは？と思っています。少々遅れましたが、写真は今年の干支、牛です。 |
記事の評価
| |
