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

同カテゴリーの次の記事

インテル® Cilk™ Plus の配列表記 (アレイ・ノーテーション) 入門

この記事は、インテル® デベロッパー・ゾーンに掲載されている「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 に適用する方法を紹介した別のブログもご覧ください。

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

関連記事

  • ガイド付き自動並列化ガイド付き自動並列化 この記事は、インテル® デベロッパー・ゾーンに掲載されている「Guided Autoparallelism」の日本語参考訳です。 インテル® Composer XE Linux* 版には、ガイド付き自動並列化 (GAP) 機能が含まれており、Fortran と C/C++ […]
  • インテル® Advisor ユーザー向けベクトル化リソースインテル® Advisor ユーザー向けベクトル化リソース この記事は、インテル® デベロッパー・ゾーンに公開されている「Vectorization Resources for Intel® Advisor Users」の日本語参考訳です。 この記事の PDF 版はこちらからご利用になれます。 インテル® Advisor 2016 には、Fortran、C および C++ […]
  • インテル® Cilk™ Plus の配列表記 (アレイ・ノーテーション) 入門インテル® Cilk™ Plus の配列表記 (アレイ・ノーテーション) 入門 この記事は、インテル® デベロッパー・ゾーンに掲載されている「Getting Started with Intel® Cilk™ Plus Array Notations」の日本語参考訳です。 はじめに配列表記 (アレイ・ノーテーション) は、通常の C/C++ […]
  • インテル® MIC アーキテクチャーの準備インテル® MIC アーキテクチャーの準備 この記事は、インテル® デベロッパー・ゾーンに掲載されている「Preparing for the Intel® Many Integrated Core Architecture」の日本語参考訳です。 インテル® MIC アーキテクチャー向けのコンパイラー手法 インテル® […]
  • ベクトル異形関数でのプロセッサー ID 指定ベクトル異形関数でのプロセッサー ID 指定 インテル® C/C++ および Fortran コンパイラーのバージョン 16.0 以降では、プログラマーがスカラー関数に対応するベクトル関数を明示的に記述することを可能にする、ベクトル異形関数 (Vector Variant Function) をサポートしています。 Windows* および Linux* […]