並列コンピューティングにおける自動ベクトル化は、自動並列化の特殊なケースであり、コンピュータプログラムが、一度に1組のオペランドを処理するスカラー実装から、複数のオペランドのペアに対して一度に1つの演算を処理するベクトル実装に変換されます。たとえば、特殊なスーパーコンピュータを含む現代の従来型コンピュータは、通常、次の4つの加算などの演算を同時に実行するベクトル演算を備えています( SIMDまたはSPMDハードウェアを介して)。
しかし、ほとんどのプログラミング言語では、多数の数値を順次加算するループを記述するのが一般的です。以下は、 C言語で記述されたそのようなループの例です。
for ( i = 0 ; i < n ; i ++ ) c [ i ] = a [ i ] + b [ i ];ベクトル化コンパイラは、このようなループをベクトル演算のシーケンスに変換します。これらのベクトル演算は、配列a、bおよびの要素のブロックに対して加算を実行しますc。自動ベクトル化は、コンピュータサイエンスにおける主要な研究テーマです。
初期のコンピュータは通常、1つの論理ユニットしか持たず、一度に1組のオペランドに対して1つの命令を実行していました。そのため、コンピュータ言語やプログラムは順次実行されるように設計されていました。しかし、現代のコンピュータは多くの処理を同時に実行できます。そこで、多くの最適化コンパイラは自動的にベクトル化を行い、逐次プログラムの一部を並列処理に変換します。
ループベクトル化は、各オペランドのペアに処理ユニットを割り当てることで、手続き型ループを変換します。プログラムは、処理時間の大部分をこのようなループ内で費やします。そのため、ベクトル化は、特に大規模なデータセットを扱う場合に、ループの処理速度を大幅に向上させることができます。ループベクトル化は、IntelのMMX、SSE、AVX、Power ISAのAltiVec、ARMのNEON、SVE、SVE2、およびRISC-Vのベクトル拡張命令セットに実装されています。
ベクトル化を妨げたり阻害したりする制約は数多く存在します。パイプライン同期やデータ移動のタイミングなどが原因で、ベクトル化によって実行速度が低下する場合もあります。ループ依存性分析は、ループ内の命令のデータ依存性に基づいて、ベクトル化可能なループを特定します。
自動ベクトル化は、ループ最適化やその他のコンパイル時最適化と同様に、プログラムの動作を正確に維持しなければならない。
誤った結果を防ぐためには、実行時にすべての依存関係を尊重する必要があります。
一般的に、ループ不変の依存関係と字句順方向の依存関係は容易にベクトル化でき、字句逆方向の依存関係は字句順方向の依存関係に変換できます。ただし、これらの変換は、すべてのステートメント間の依存関係が元の状態を維持するように、安全に行う必要があります。
循環依存関係は、ベクトル化された命令とは独立して処理されなければならない。
コンパイラが過度に慎重に依存関係を想定する場合があります。一部のコンパイラは、依存関係を無視するようにコンパイラに指示するivdepというディレクティブを提供しています。このようなディレクティブを誤って適用すると、誤った結果が生じます。[ 1 ]
ベクトル命令の実行中は、整数精度(ビットサイズ)を維持する必要があります。適切なベクトル命令は、内部整数のサイズと動作に基づいて選択する必要があります。また、異なる整数型が混在する場合は、精度を損なうことなく正しく昇格/降格を行うよう、特に注意が必要です。符号拡張(複数の整数が同じレジスタ内に格納されるため)や、シフト演算、またはキャリービットを考慮した演算を行う際には、特に注意が必要です。
デフォルトではIEEE-754の浮動小数点セマンティクスが維持されるため、ベクトル化が全く不可能になる場合がよくあります。配列の合計を計算する以下の単純なループを考えてみましょう。
float sum ( float * A , int n ) { float sum = 0 ; for ( int i = 0 ; i < n ; ++ i ) sum += A [ i ]; return sum ; }浮動小数点演算は結合法則を満たさないため(加算の順序を変更すると結果が変わるため)、コンパイラはループをベクトル化することはできません。コンパイラが明示的に再結合(演算を結合法則を満たすかのように並べ替える)を許可された場合にのみ、次のようなコードにコンパイルできます。[ 2 ] : §削減
float sum ( float * A , int n ) { float4 sum4 = { 0 }, * A4 = A ; int n4 = n / 4 ; for ( int i = 0 ; i < n4 ; ++ i ) { sum4 += A4 [ i ]; }float sum = sum4 [ 0 ] + sum4 [ 1 ] + sum4 [ 2 ] + sum4 [ 3 ]; for ( int i = n / 4 ; i < n ; ++ i ) { sum += A [ i ]; }合計を返す; }この場合、結果は元の値と若干異なりますが、漸近的な丸め誤差はほぼ同じです。これは削減の一例です。
IEEE-754規格を無視しても、大きな差異が生じる場合は、通常、プログラマーのミスを示しています。よくある原因としては、完全な再関連付けによって、補正加算アルゴリズムが機能しなくなることが挙げられます。OpenMPのreduction句のような、より制限の厳しいコマンドを使用することで、より的確な設定が可能になります。
プログラムのベクトル化を行うには、コンパイラのオプティマイザがまずステートメント間の依存関係を理解し、必要に応じてそれらを再配置する必要があります。依存関係がマッピングされたら、オプティマイザは実装命令を適切に配置し、適切な候補を複数のデータ項目を操作するベクトル命令に変換する必要があります。
最初のステップは、依存関係グラフを作成し、どのステートメントが他のどのステートメントに依存しているかを特定することです。これには、各ステートメントを調べ、そのステートメントがアクセスするすべてのデータ項目を特定し、配列アクセス修飾子を関数にマッピングし、すべてのステートメントにおけるすべてのアクセスの依存関係をチェックすることが含まれます。エイリアス分析を使用すると、異なる変数がメモリ内の同じ領域にアクセス(または交差)していることを検証できます。
依存関係グラフには、ベクトルサイズ以下の距離を持つすべてのローカル依存関係が含まれます。したがって、ベクトルレジスタが128ビット、配列型が32ビットの場合、ベクトルサイズは128/32 = 4となります。同じベクトル命令内で同時アクセスが発生しないため、その他の非循環依存関係によってベクトル化が無効になることはありません。
ベクトルのサイズが4つの整数と同じだと仮定します。
for ( i = 0 ; i < 128 ; i ++ ) { a [ i ] = a [ i -16 ]; // 16 > 4 なので無視しても安全a [ i ] = a [ i -1 ]; // 1 < 4 なので依存関係グラフ上に残る}グラフを使用することで、オプティマイザは強連結成分(SCC)をクラスタリングし、ベクトル化可能なステートメントをそれ以外のステートメントから分離することができます。
例えば、ループ内に3つのステートメントグループ(SCC1+SCC2、SCC3、SCC4)を含むプログラム断片を考えてみましょう。このうち、ベクトル化できるのは2番目のグループ(SCC3)のみです。最終的なプログラムには、各グループに対応する3つのループが含まれますが、ベクトル化されるのは真ん中のループのみです。オプティマイザは、ステートメントの実行順序を破ることなく最初のループと最後のループを結合することはできません。実行順序を破ると、必要な保証が無効になってしまうからです。
一見分かりにくい依存関係の中には、特定の慣用表現に基づいてさらに最適化できるものがある。
例えば、以下の自己データ依存性はベクトル化できます。なぜなら、右辺の値(RHS)の値が取得され、左辺の値に格納されるため、代入内でデータが変更されることはないからです。
a [ i ] = a [ i ] + a [ i + 1 ];スカラーによる自己依存性は、変数消去によってベクトル化できる。
ループベクトル化の一般的なフレームワークは、次の4つの段階に分けられます。
ベクトル化の中には、コンパイル時に完全にチェックできないものもあります。例えば、ライブラリ関数は、処理するデータが呼び出し元から提供される場合、最適化を無効にする可能性があります。しかし、このような場合でも、実行時最適化によってループを動的にベクトル化することは可能です。
この実行時チェックは前処理段階で行われ、可能であればベクトル化命令に処理の流れを誘導し、そうでなければレジスタまたはスカラー変数に渡される変数に応じて標準処理に戻ります。
以下のコードは外部パラメータに依存しないため、コンパイル時に容易にベクトル化できます。また、これらの変数はローカル変数であり、実行スタック内にのみ存在するため、他の変数と同じメモリ領域を占有しないことが言語によって保証されます。
int a [ 128 ]; int b [ 128 ]; // b を初期化for ( i = 0 ; i < 128 ; i ++ ) a [ i ] = b [ i ] + 5 ;一方、以下のコードにはメモリ位置に関する情報がありません。なぜなら、参照はポインタであり、ポインタが指すメモリ領域が重複する可能性があるからです。
void compute ( int * a , int * b ) { int i ; for ( i = 0 ; i < 128 ; i ++ , a ++ , b ++ ) * a = * b + 5 ; }aとbのアドレス、およびループ反復空間(128)をランタイム時に簡単にチェックするだけで、配列が重複しているかどうかが分かり、依存関係が明らかになります。(C99以降では、restrictキーワードでパラメータを修飾すると(ここでは)、 aとbが指すメモリ範囲が重複しないことがコンパイラに伝わり、上記の例と同じ結果が得られます。)int *restrict a, int *restrict b
既存のアプリケーションを動的に分析して、コンパイラのさらなる進歩や手動コードの変更によって活用できるSIMD並列処理の潜在的な可能性を評価するツールがいくつか存在する。[ 3 ]
例えば、2つの数値データベクトルを乗算するプログラムが挙げられます。スカラー的なアプローチは次のようになります。
for ( i = 0 ; i < 1024 ; i ++ ) c [ i ] = a [ i ] * b [ i ];これは、次のようなベクトル化が可能です。
for ( i = 0 ; i < 1024 ; i += 4 ) c [ i : i + 3 ] = a [ i : i + 3 ] * b [ i : i + 3 ];ここで、c[i:i+3]はc[i]からc[i+3]までの4つの配列要素を表し、ベクトルプロセッサは1つのベクトル命令で4つの演算を実行できます。4つのベクトル演算は1つのスカラー命令とほぼ同じ時間で完了するため、ベクトル方式は元のコードよりも最大4倍高速に実行できます。
コンパイラのアプローチには大きく分けて2種類あります。1つは従来型のベクトル化技術に基づくもので、もう1つはループ展開に基づくものです。
この手法は、従来型ベクトルマシンで用いられており、ループレベルでSIMD並列性を見つけて活用しようとするものです。主な手順は以下の2つです。
最初のステップでは、コンパイラはベクトル化を妨げる障害を探します。ベクトル化における大きな障害の一つは、ベクトルの長さよりも短い真のデータ依存性です。その他の障害としては、関数呼び出しや反復回数の少なさなどが挙げられます。
ループがベクトル化可能であると判断されたら、ループはベクトルの長さに基づいてストリップマイニングされ、ループ本体内の各スカラー命令が対応するベクトル命令に置き換えられます。以下に、上記の例を使用して、このステップにおけるコンポーネント変換を示します。
for ( i = 0 ; i < 1024 ; i += 4 ) for ( j = 0 ; j < 4 ; j ++ ) c [ i + j ] = a [ i + j ] * b [ i + j ];for ( i = 0 ; i < 1024 ; i += 4 ) { for ( j = 0 ; j < 4 ; j ++ ) tA [ j ] = A [ i + j ]; for ( j = 0 ; j < 4 ; j ++ ) tB [ j ] = B [ i + j ]; for ( j = 0 ; j < 4 ; j ++ ) tC [ j ] = tA [ j ] * tB [ j ]; for ( j = 0 ; j < 4 ; j ++ ) C [ i + j ] = tC [ j ]; }for ( i = 0 ; i < 1024 ; i += 4 ) { vA = vec_ld ( & A [ i ]); vB = vec_ld ( & B [ i ]); vC = vec_mul ( vA , vB ); vec_st ( vC , & C [ i ]); }この比較的新しい技術は、特にベクトル長が短い最新のSIMDアーキテクチャを対象としています。[ 4 ]ループを展開して基本ブロック内のSIMD並列性を増やすことはできますが、この技術はループではなく基本ブロック内のSIMD並列性を活用します。主な手順は次の2つです。
この手法における段階的な変換を示すために、同じ例を再度使用します。
for ( i = 0 ; i < 1024 ; i += 4 ) { sA0 = ld ( & A [ i + 0 ]); sB0 = ld ( & B [ i + 0 ]); sC0 = sA0 * sB0 ; st ( sC0 , & C [ i + 0 ]); ... sA3 = ld ( & A [ i + 3 ]); sB3 = ld ( & B [ i + 3 ]); sC3 = sA3 * sB3 ; st ( sC3 , & C [ i + 3 ]); }for ( i = 0 ; i < 1024 ; i += 4 ) { ( sA0 , sA1 , sA2 , sA3 ) = ld ( & A [ i + 0 : i + 3 ]); ( sB0 , sB1 , sB2 , sB3 ) = ld ( & B [ i + 0 : i + 3 ]); ( sC0 , sC1 , sC2 , sC3 ) = ( sA0 , sA1 , sA2 , sA3 ) * ( sB0 , sB1 , sB2 , sB3 ); st (( sC0 , sC1 , sC2 , sC3 ), & C [ i + 0 : i + 3 ]); }for ( i = 0 ; i < 1024 ; i += 4 ) { vA = vec_ld ( & A [ i ]); vB = vec_ld ( & B [ i ]); vC = vec_mul ( vA , vB ); vec_st ( vC , & C [ i ]); }ここで、sA1、sB1、...はスカラー変数を表し、vA、vB、vCはベクトル変数を表します。
ほとんどの自動ベクトル化商用コンパイラは従来型のループレベルのアプローチを使用していますが、IBM XL Compiler [ 5 ]は両方を使用しています。
ループ本体に if 文が存在すると、変数の複数の値をマージするために、すべての制御パスで命令を実行する必要があります。一般的なアプローチの 1 つは、コード変換のシーケンスを実行することです。述語化 → ベクトル化 (上記の方法のいずれかを使用) → ベクトル述語の削除 → スカラー述語の削除。[ 6 ]次のコードを例としてこれらの変換を示します。
for ( i = 0 ; i < 1024 ; i ++ ) if ( A [ i ] > 0 ) C [ i ] = B [ i ]; else D [ i ] = D [ i -1 ];for ( i = 0 ; i < 1024 ; i ++ ) { P = A [ i ] > 0 ; NP = ! P ; C [ i ] = B [ i ]; ( P ) D [ i ] = D [ i -1 ]; ( NP ) }ここで、(P)はステートメントを保護する述語を表す。
for ( i = 0 ; i < 1024 ; i += 4 ) { vP = A [ i : i + 3 ] > ( 0 , 0 , 0 , 0 ); vNP = vec_not ( vP ); C [ i : i + 3 ] = B [ i : i + 3 ]; ( vP ) ( NP1 , NP2 , NP3 , NP4 ) = vNP ; D [ i + 3 ] = D [ i + 2 ]; ( NP4 ) D [ i + 2 ] = D [ i + 1 ]; ( NP3 ) D [ i + 1 ] = D [ i ]; ( NP2 ) D [ i ] = D [ i -1 ]; ( NP1 ) }for ( i = 0 ; i < 1024 ; i += 4 ) { vP = A [ i : i + 3 ] > ( 0 , 0 , 0 , 0 ); vNP = vec_not ( vP ); C [ i : i + 3 ] = vec_sel ( C [ i : i + 3 ], B [ i : i + 3 ], vP ); ( NP1 , NP2 , NP3 , NP4 ) = vNP ; D [ i + 3 ] = D [ i + 2 ]; ( NP4 ) D [ i + 2 ] = D [ i + 1 ]; ( NP3 ) D [ i + 1 ] = D [ i ]; ( NP2 ) D [ i ] = D [ i -1 ]; ( NP1 ) }for ( i = 0 ; i < 1024 ; i += 4 ) { vP = A [ i : i + 3 ] > ( 0 , 0 , 0 , 0 ); vNP = vec_not ( vP ); C [ i : i + 3 ] = vec_sel ( C [ i : i + 3 ], B [ i : i + 3 ], vP ); ( NP1 , NP2 , NP3 , NP4 ) = vNP ; if ( NP4 ) D [ i + 3 ] = D [ i + 2 ]; if ( NP3 ) D [ i + 2 ] = D [ i + 1 ]; if ( NP2 ) D [ i + 1 ] = D [ i ]; if ( NP1 ) D [ i ] = D [ i -1 ]; }ベクトルコードでは、すべての制御パスで命令を実行する必要があるため、スカラーベースラインと比較してベクトルコードの速度が低下する主な要因の 1 つは、このことにあります。制御フローが複雑になり、スカラーコードでバイパスされる命令が増えるほど、ベクトル化のオーバーヘッドは大きくなります。このベクトル化のオーバーヘッドを減らすために、スカラー分岐がスカラー命令をバイパスするのと同様に、ベクトル分岐を挿入してベクトル命令をバイパスすることができます。[ 7 ]以下では、AltiVec 述語を使用して、これがどのように実現できるかを示します。
for ( i = 0 ; i < 1024 ; i ++ ) { if ( A [ i ] > 0 ) { C [ i ] = B [ i ]; if ( B [ i ] < 0 ) D [ i ] = E [ i ]; } }for ( i = 0 ; i < 1024 ; i += 4 ) { vPA = A [ i : i + 3 ] > ( 0 , 0 , 0 , 0 ); C [ i : i + 3 ] = vec_sel ( C [ i : i + 3 ], B [ i : i + 3 ], vPA ); vT = B [ i : i + 3 ] < ( 0 , 0 , 0 , 0 ); vPB = vec_sel (( 0 , 0 , 0 , 0 ), vT , vPA ); D [ i : i + 3 ] = vec_sel ( D [ i : i + 3 ], E [ i : i + 3 ], vPB ); }for ( i = 0 ; i < 1024 ; i += 4 ) { if ( vec_any_gt ( A [ i : i + 3 ], ( 0 , 0 , 0 , 0 ))) { vPA = A [ i : i + 3 ] > ( 0 , 0 , 0 , 0 ); C [ i : i + 3 ] = vec_sel ( C [ i : i + 3 ], B [ i : i + 3 ], vPA ); vT = B [ i : i + 3 ] < ( 0 , 0 , 0 , 0 ); vPB = vec_sel (( 0 , 0 , 0 , 0 ), vT , vPA ); if ( vec_any_ne ( vPB , ( 0 , 0 , 0 , 0 ))) D [ i : i + 3 ] = vec_sel ( D [ i : i + 3 ], E [ i : i + 3 ], vPB ); } }ベクトル分岐を含む最終コードには、注意すべき点が2つあります。1つ目は、vec_any_gtを使用することで、vPAの述語定義命令が外側のベクトル分岐の本体内にも含まれていることです。2つ目は、vPBの内側のベクトル分岐の収益性は、vPAのすべてのフィールドが偽値であるという条件の下で、vPBのすべてのフィールドが偽値である条件付き確率に依存していることです。
スカラーベースラインの外側の分岐が常に実行され、ループ本体のほとんどの命令がスキップされる例を考えてみましょう。上記の中間ケースでは、ベクトル分岐がないため、すべてのベクトル命令が実行されます。ベクトル分岐を含む最終コードでは、比較と分岐の両方がベクトルモードで実行されるため、スカラーベースラインよりもパフォーマンスが向上する可能性があります。
ほとんどのCおよびC++コンパイラでは、組み込み関数を使用してSIMDを手動で使用できますが、プログラマの労力、保守性、移植性が犠牲になります。一部の言語(GNU C、C++ std::experimental::simd、Rustなどstd::simd)には、適切なSIMD命令にコンパイルされるベクトルデータ型が含まれており、移植性が向上し、必要な労力が削減されます。[ 8 ]
別のアプローチはSPMDです。これは、一度に 1 つの要素のみを操作するように見えるプログラムを作成し、コンパイラに SIMD ベクトルの幅に合うように拡張させるものです。これはグラフィックスシェーダーで使用されているアプローチであり、最近では Intel IPSC などの CPU 指向のツールでも採用されています。ループ内で外部関数が呼び出される場合など、失敗してスカラー コードにフォールバックする可能性がある自動ベクトル化とは異なり、[ a ] SPMD は、組み込み関数やベクトルデータ型を手動で使用する場合と同様に、適用可能な場合は必ずベクトル コードが生成されることが保証されています。[ 9 ]