配列表記 (アレイ・ノーテーション) による SIMD 並列処理

この記事は、インテル® デベロッパー・ゾーンに掲載されている「SIMD Parallelism using Array Notation」の日本語参考訳です。

この記事は、APL や Fortran 90 の配列表記を使ってみたいと考えたことがある C/C++ 開発者、そして配列表記をまだご存知ない開発者向けに、インテル® Parallel Composer の機能である C/C++ で利用可能な配列表記について説明します。

背景

以前、私は並列プログラミング向けの Three Layer Cake パターンに関する論文を執筆しました。パターンとは、近年のマルチコアチップを最大限に活用するため、プログラムを構造化する方法です。次の 2 つのレイヤーで構成されます。

• fork-join: 複数のハードウェア・スレッドを利用します。

• SIMD: SIMD 命令を利用します。

インテル® Parallel Composer に含まれるコンパイラーは C++ を拡張し、この 2 つのレイヤーを直接サポートしています。この拡張はインテル® Cilk™ Plus と呼ばれ、次の機能を利用できます。

• fork-join 並列を指定する Cilk 表記

• SIMD 並列を指定する配列表記 (アレイ・ノーテーション)

ここでは、Seismic Duck (http://home.comcast.net/~arch.robison/seismic_duck.html) カーネルの例を使って配列表記について説明し、Cilk 表記については別のブログで説明します。この 2 つの表記は独立しています。配列表記は、インテル® スレッディング・ビルディング・ブロック (インテル® TBB) など他のスレッド化ツールを使用している場合や、高速なシリアルコードを記述する際にも役立ちます。

配列表記の概要

配列表記の拡張は、APL や Fortran 90 形式の配列表記を連想させます。例えば、次の配列表記について考えてみます。

a[index:count]

これは、index から始まる count 個の要素からなる配列セクションを表しています。   この配列セクションに対して直感的な方法でスカラー演算を適用できます。スカラー間および配列セクション間の操作も可能です。スカラーは明白な方法で (APL や Fortran 90 のように) 拡張されます。次に例を示します。

z[i:n] = x[i:n];      // x[i..i+n-1] を z[i..i+n-1] へコピーする 

z[i:n] = 2*x[i+1:n];  // z[i..i+n-1] を x[i+1..i+n] の対応する要素の 2 倍に設定する 

u[i:m][j:n] += 1;     // 2 次元の m x n 配列セクションの部分要素 [i][j]をインクリメントする

セクション表記では、[index:count:stride] 形式の配列、リダクション、短縮形も使えます。ここでは、配列表記の概要を示すことが目的であるため、これらについては取り上げません。詳細は、コンパイラーのドキュメントを参照してください。

例

別のブログで、Seismic Duck の地震波伝搬が「タイル」 (キャッシュに収まる小さな部分配列) の更新にどのように依存しているかを述べましたが、以下は実行時間の大半を占めるスカラーコードです。均一の係数 A と B を使ってタイルを更新します。

float a = 2*A[iFirst][jFirst]; 

float b = B[iFirst][jFirst]; 

for( int i=iFirst; i<iLast; ++i ) { 

    for( int j=jFirst; j<jLast; ++j ) { 

         Vx[i][j] += a*(U[i][j+1]-U[i][j]); 

         Vy[i][j] += a*(U[i+1][j]-U[i][j]); 

         U[i][j] += b*((Vx[i][j]-Vx[i][j-1])+(Vy[i][j]-Vy[i-1][j])); 

    } 

}

コンパイラーの自動ベクトル化により SIMD コードが生成されないスカラーループの速度を向上するため、SSE 組込み関数を用いて主要なループを記述し、一度に実行される計算が 1 つから 4 つになるようにコードを変更しました。以下は変更後のコードです。

#define CAST(x) (*(__m128*)&(x)) /* アライメントされたロードまたはストア */ 

#define LOAD(x) _mm_loadu_ps(&(x)) /* アライメントされていないロード */ 

#define ADD _mm_add_ps 

#define MUL _mm_mul_ps 

#define SUB _mm_sub_ps 

... 

__m128 a = CAST(A[iFirst][jFirst]); 

a = ADD(a,a); 

__m128 b = CAST(B[iFirst][jFirst]); 

for( int i=iFirst; i<iLast; ++i ) { 

    for( int j=jFirst; j<jLast; j+=4 ) { 

        __m128 u = CAST(U[i][j]); 

        CAST(Vx[i][j]) = ADD(CAST(Vx[i][j]),MUL(a,SUB(LOAD(U[i][j+1]),u))); 

        CAST(Vy[i][j]) = ADD(CAST(Vy[i][j]),MUL(a,SUB(CAST(U[i+1][j]),u))); 

        CAST(U[i][j]) = ADD(u,MUL(b,ADD(SUB(CAST(Vx[i][j]),LOAD(Vx[i][j-1])),SUB(CAST(Vy[i][j]),CAST(Vy[i-1][j]))))); 

    } 

}

この変更によりコードが複雑になりましたが、ほかの部分のロジックは jLast-jFirst が 4 の倍数になることを保証しているため、これは単純なケースと言えます。そうでなければ、倍数で割り切れない反復の余りを処理するため、コードがさらに難解なものになったでしょう。

この例は、コンパイラーの自動でベクトル化 により SIMD 命令に変換されるため、実際には明示的な SSE 組込み関数を記述する必要ありません。最近のコンパイラーで試してみたところ、自動でベクトル化されました  (ただし、2008 年にリリースされたバージョンのコンパイラーでは自動でベクトル化されませんでした)。ここでは、オプティマイザーを支援するため、配列 Vx、Vy、U をポインターではなく、ファイルスコープの静的配列として宣言しています。これはオブジェクト指向プログラミングではありませんが、コンパイラーはエイリアシングが存在しないこと、つまりベクトル化の妨げになるループ伝搬依存が存在しないことを証明できます。このような方法でオプティマイザーを支援するのは、常に現実的であるとは限りません。さらに、配列表記はそれ自体が洗練されています。ここでは、一例としてアルゴリズムのカーネル部分を使用します。

インテル® Cilk™ Plus の配列表記は、開発者の意図 (「SIMD 並列化」) をより明確にコンパイラーに伝えます。以下のコードは、配列表記を利用した例です。

int i = iFirst;

int j = jFirst;

size_t m = iLast-iFirst;

size_t n = jLast-jFirst;

float a = 2*A[i][j];

float b = B[i][j];

Vx[i:m][j:n] += a*(U[i:m][j+1:n]-U[i:m][j:n]);

Vy[i:m][j:n] += a*(U[i+1:m][j:n]-U[i:m][j:n]);

U[i:m][j:n] += b*((Vx[i:m][j:n]-Vx[i:m][j-1:n])+(Vy[i:m][j:n]-Vy[i-1:m][j:n]));

最後の 3 行と、次のオリジナルの forall ループ (別のブログ) を比較してください。

forall i, j { 

    Vx[i][j] += (A[i][j+1]+A[i][j])*(U[i][j+1]-U[i][j]); 

    Vy[i][j] += (A[i+1][j]+A[i][j])*(U[i+1][j]-U[i][j]); 

} 

forall i, j { 

    U[i][j] += B[i][j]*((Vx[i][j]-Vx[i][j-1]) + (Vy[i][j]-Vy[i-1][j])); 

}

配列表記により、明確に並列化の意図を伝えることができます。この例では i、j、m、n をセットアップするコードを追加していますが、これは前述の単純な C++ for ループの例の名残です。コードのほかの部分では等価な {i, j, m, n} から {iFirst, iLast, jFirst, jLast} を計算しています。残りのコードでも配列表記を使うのであれば、{iFirst, iLast, jFirst, jLast} をすべて削除し、代わりに {i,j,m,n} を使用することができます。

一時部分配列の削除はオプティマイザーに任せています。例えば、次の部分配列表記について考えてみます。

a*(U[i:m][j+1:n]-U[i:m][j:n]);

これは、概念的には “-” および “*” 演算の結果を格納するため 2 つの一時部分配列を生成します。インテル® コンパイラーは、このような一時部分配列を削除します   (これに関しては、Fortran 90 で長年の実績があります)。ただし、一部の一時部分配列が残ったとしても、並列化は明確です。コンパイラーは、シリアル for ループの依存性解析から並列化を推測する必要はありません。

簡潔に表記できるのは良いことですが、パフォーマンスはどうでしょうか?   インテル® コンパイラーでコンパイルし、インテル® Core™2 Quad プロセッサー・ベースのシステムで実行してみたところ、配列表記バージョンのほうが、手動で変更した SSE バージョンよりも速いという結果になりました。[パフォーマンスは異なることがあります。最適化に関する注意事項を参照してください。] 原因を調査するため、オブジェクト・コードを検証したところ、Vx、Vy、U の更新は個別のループで行ったほうが、1 つの大きなループで行うよりもパフォーマンスが良いことが分かりました。手動で変更した SSE バージョンでも、この 3 つの配列の更新を個別のループで行うようにしたところ、パフォーマンスが向上しました。この例で、配列表記を利用した場合と同様のパフォーマンスを手動でも実現できたことを嬉しく思います。

まとめ

配列表記は、SIMD 並列化を簡潔に表現することができます。ほかのコンパイラーでも配列表記がサポートされることを期待します。インテル® Cilk™ Plus の fork-join 並列化を Seismic Duck に適用する方法を紹介した別のブログもご覧ください。

コンパイラーの最適化に関する詳細は、最適化に関する注意事項を参照してください。