はじめに

前回の記事で、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の配列を用意し、例えば (x₀, x₁, x₂, x₃) と、粒子のx座標のみをvector float にパックしたStructure Of Array (SOA)　フォームとが考えられますが、計算効率を考えると、SOAフォームが最適です。ただしSOAフォームを取る場合は、粒子数が4で割り切れない場合は、余りの処理をループの外に別途書く必要があります。今回の記事ではそのあたりはサボっておりますので、注意が必要です。

では、早速コードをみてみましょう。

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による手作業の並列化に勝るコードを書くことができることがわかりました。

次に、SPEの数を変化させた場合の処理時間の変化を見てみます。ただし、粒子数8192, Step時間300 の条件では処理が軽すぎて、SPEの数を変えてもあまり処理時間に変化は見られなかったので、少し負荷を増やし、

• 粒子数: 1048576 個
• Step時間: 1024 sec
• ハードウエア: IBM QS22

とします。

今回はIBM BladeCenter QS22 で検証を行っているので、SPEは最大16個使うことができます。 使用SPEの数は、環境変数CSS_NUM_SPUSで決定されます。例えばSPEの数を6個使用したい場合は

$export CSS_NUM_SPUS=6

とコマンドを打ちます。特に環境変数を設定していない場合は、デフォルトとしてSPEが8個使用されます。さて、この環境変数CSS_NUM_SPUSを1から16まで変えながら処理時間を測定した結果を図2に示します。

ここで、Speed up factor とは、SPEが1個の場合の処理時間を1とした場合の処理速度向上の度合いで、具体的にはSPEが1個の場合の処理速度をSPEを変化させた場合の処理速度で割った値です。SPEの数の増加にしたがって、綺麗に処理時間が減少しているのがわかります。Euler粒子シミュレーションプログラムは、粒子間の相互作用がなく、完全な並列演算ができるため、もともと計算資源の数にしたがってリニアにパフォーマンスが出やすいのですが、この結果はCellSsによる各種オーバーヘッドが非常に少ないことを示唆しています。

まとめ

今回の記事では、Intrinsics を使用し、Euler粒子シミュレーションプログラムを高速化することに成功しました。自動並列化ソリューションは、性能面がいま一つで、実際は実用的でないものが多い印象がありますが、CellSsでは、SPEコードの中身を制約なくいじることができるので、手作業の並列化コードに近い性能を引き出すことができます。実際、上記コードをループアンローリング、ソフトウエアパイプライニング、等の手法を用いて、さらに高速化することも可能だと思います。是非、チャレンジしてみてください。

参考文献

学ぶために

• CellSs の開発チームが執筆した IBM Journal, CellSs: Programming the Cell/B.E. made easier, には、CellSsのコンセプトと性能が簡潔にまとめられています。
  
  
• Jonathon Bartlett の『Cell BE プロセッサーでのハイパフォーマンス・アプリケーションのプログラミング』シリーズ (developerWorks 2007年1月) には、PS3でのLinux入門, PS3上のSPEプログラミング, an SPUの紹介, SPU パフォーマンスプログラミング, C/C++ でのSPU プログラミング, and DMA転送における効率的なバッファー管理が収録されています。
  
  
• IBM Semiconductor Solutions Technical LibraryのCell Broadband Engine documentation セクションでは、有益なマニュアルや仕様書などをダウンロードできます。
  
  
• 全ての Cell/B.E. 情報 -- 関連記事、ディスカッションフォーラム、ダウンロードサイトなど -- は IBM developerWorks のCell Broadband Engine resource center: your definitive resource for all things Cell/B.E.にあります。
  
  

製品や技術を入手するために

• Balcelona Super Computing Center のCellSs サイトはこちらです 。
  
  
• コンパイル済みのCellSs が収録されているYellow Dog Linux 6.1 は、こちらからフリーでダウンロード可能です。
  
  
• Yellow Dog LinuxのPlayStation 3への導入については、ZDNetに詳細な解説記事が掲載されています。
  
  
• Cell/B.E. を手に入れるにはこちらから: Contact IBM about custom Cell-based or custom-processor based solutions.
  
  

議論するために

• IBMのエキスパートやdeveloperWorksの実践者達から素早い回答を得るには Cell Broadband Engine discussion forumに参加してください。
  
  
• The Cell Broadband Engine/Power Architecture notebook はブログベースの情報サイトで、ニュースや、2つの教育的な資料集: フォーラムでの面白い質問や議論の熱いトピックスをまとめた"Forum watch"や、"Infobomb" シリーズ (Cell/B.E.ドキュメント群から、短くて、正確で、特定のタスクに注目できて、斜め読みに適した"ボム"が集められている) が利用できます。
  
  

著者について

浅原明広は、Cell/B.E.の発表当初から関連の開発業務に携わってきました。多数のハードウエア、多数のプログラミングモデルが混在し、マルチコアプロセッサ業界は今が旬。市況はさんざんですが、エンジニアにとっては今年は非常に面白い年になるのでは？と思っています。少々遅れましたが、写真は今年の干支、牛です。

この記事を評価する

コメント

上に戻る

目次

• はじめに
• Intrinsics
• まとめ
• 参考文献
• 著者について
• コメント