x86_64で以下をベクトル化するためにどの組み込み関数を使用しますか（ベクトル化が可能である場合）？
double myNum = 0;
for(int i=0;i<n;i++){
myNum += a[b[i]] * c[i]; //b[i] = int, a[b[i]] = double, c[i] = double
}
1184 次
5 に答える
8
これが私の試練であり、完全に最適化され、テストされています。
#include <emmintrin.h>
__m128d sum = _mm_setzero_pd();
for(int i=0; i<n; i+=2) {
sum = _mm_add_pd(sum, _mm_mul_pd(
_mm_loadu_pd(c + i),
_mm_setr_pd(a[b[i]], a[b[i+1]])
));
}
if(n & 1) {
sum = _mm_add_pd(sum, _mm_set_sd(a[b[n-1]] * c[n-1]));
}
double finalSum = _mm_cvtsd_f64(_mm_add_pd(
sum, _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1))
));
gcc -O2 -msse2これにより、 （4.4.1）を使用して非常に美しいアセンブリコードが生成されます。
お分かりのように、偶数nを使用すると、このループが高速になり、整列されcます。調整できる場合はc、に変更_mm_loadu_pdし_mm_load_pdて実行時間をさらに短縮します。
于 2010-03-02T01:54:03.710 に答える
4
ループを展開することから始めます。何かのようなもの
double myNum1 = 0, myNum2=0;
for(int i=0;i<n;i+=2)
{
myNum1 += a[b[ i ]] * c[ i ];
myNum2 += a[b[i+1]] * c[i+1];
}
// ...extra code to handle the remainder when n isn't a multiple of 2...
double myNum = myNum1 + myNum2;
うまくいけば、コンパイラが負荷を算術でインターリーブできるようになります。プロファイルを作成し、アセンブリを見て、改善があるかどうかを確認します。理想的にはコンパイラがSSE命令を生成しますが、実際にはそうではありません。
さらに展開すると、これが可能になる場合があります。
__m128d sum0, sum1;
// ...initialize to zero...
for(int i=0;i<n;i+=4)
{
double temp0 = a[b[ i ]] * c[ i ];
double temp1 = a[b[i+1]] * c[i+1];
double temp2 = a[b[i+2]] * c[i+2];
double temp3 = a[b[i+3]] * c[i+3];
__m128d pair0 = _mm_set_pd(temp0, temp1);
__m128d pair1 = _mm_set_pd(temp2, temp3);
sum0 = _mm_add_pd(sum0, pair0);
sum1 = _mm_add_pd(sum1, pair1);
}
// ...extra code to handle the remainder when n isn't a multiple of 4...
// ...add sum0 and sum1, then add the result's components...
（開始時と終了時の擬似コードについてお詫びします。重要な部分はループだったと思います）。それが速くなるかどうかはわかりません。これは、さまざまなレイテンシーと、コンパイラーがすべてをどれだけうまく再配置できるかに依存します。実際に改善があったかどうかを確認するために、前後にプロファイルを作成してください。
お役に立てば幸いです。
于 2010-02-28T11:28:16.880 に答える
2
Intelプロセッサは2つの浮動小数点演算を発行できますが、サイクルごとに1つの負荷がかかるため、メモリアクセスが最も厳しい制約になります。そのことを念頭に置いて、私は最初にパックロードを使用してロード命令の数を減らすことを目指し、便利であるという理由だけでパック演算を使用しました。それ以来、メモリ帯域幅の飽和が最大の問題である可能性があることに気付きました。ベクトル化を学習するのではなく、コードを高速化することがポイントである場合、SSE命令をいじくり回すのは時期尚早の最適化であった可能性があります。
SSE
のインデックスを想定せずに可能な限り少ない負荷ではb、ループを4回展開する必要があります。1つの128ビットロードはから4つのインデックスを取得しb、2つの128ビットロードはそれぞれから隣接するダブルのペアを取得し、必要な独立した64ビットロードcを収集します。aこれは、シリアルコードの4回の反復あたり7サイクルのフロアです。a（へのアクセスが適切にキャッシュされない場合は、メモリ帯域幅を飽和させるのに十分です）。4の倍数ではない反復の数を処理するなど、いくつかの厄介なことを省略しました。
entry: ; (rdi,rsi,rdx,rcx) are (n,a,b,c)
xorpd xmm0, xmm0
xor r8, r8
loop:
movdqa xmm1, [rdx+4*r8]
movapd xmm2, [rcx+8*r8]
movapd xmm3, [rcx+8*r8+8]
movd r9, xmm1
movq r10, xmm1
movsd xmm4, [rsi+8*r9]
shr r10, 32
movhpd xmm4, [rsi+8*r10]
punpckhqdq xmm1, xmm1
movd r9, xmm1
movq r10, xmm1
movsd xmm5, [rsi+8*r9]
shr r10, 32
movhpd xmm5, [rsi+8*r10]
add r8, 4
cmp r8, rdi
mulpd xmm2, xmm4
mulpd xmm3, xmm5
addpd xmm0, xmm2
addpd xmm0, xmm3
jl loop
インデックスを取得することは最も複雑な部分です。movdqa16バイトに整列されたアドレスから128ビットの整数データをロードします（Nehalemには、「整数」と「浮動」のSSE命令を混合するための遅延ペナルティがあります）。punpckhqdq上位64ビットを下位64ビットに移動しますが、より単純な名前の。とは異なり、整数モードですmovhlpd。32ビットシフトは汎用レジスタで実行されます。movhpd下部を乱すことなく、xmmレジスタの上部に1つのdoubleをロードします。これは、aパックされたレジスタに直接の要素をロードするために使用されます。
このコードは、上記のコードよりも明らかに高速であり、単純なコードよりも高速であり、すべてのアクセスパターンB[i] = iで、ナイーブループが実際に最速である単純なケースです。SUM(A(B(:)),C(:))また、Fortranの関数のようないくつかのことを試しましたが、これは基本的に単純なループと同等になりました。
4モジュールで4GBのDDR2-667メモリを搭載したQ6600（2.4Ghzで65 nm Core 2）でテストしました。メモリ帯域幅をテストすると、約5333 MB / sが得られるため、1つのチャネルしか表示されていないようです。Debianのgcc4.3.2-1.1、-O3 -Ffast-math -msse2 -Ftree-vectorize -std=gnu99でコンパイルしています。
テストのために、私はn100万になり、配列を初期化してa[b[i]]、c[i]両方とも等しく1.0/(i+1)、いくつかの異なるパターンのインデックスを使用します。1つaは100万個の要素で割り当てbてランダム順列に設定し、もう1つaは10M個の要素で割り当てて10回ごとに使用し、最後の1つは10M個の要素で割り当てて、1から9までの乱数を追加してa設定します。1回の呼び出しにかかる時間を計り、配列を呼び出してキャッシュをクリアし、各関数の1000回の試行を測定します。基準の内臓からのコード（特に、パッケージ内のカーネル密度推定器）を使用して、平滑化された実行時分布をプロットしました。b[i+1]b[i]gettimeofdayclflushstatistics
帯域幅
ここで、帯域幅に関する実際の重要な注意事項について説明します。2.4Ghzクロックで5333MB/sは、1サイクルあたり2バイト強です。私のデータは十分に長いので、何もキャッシュできません。すべてが失敗した場合、ループの実行時間を反復ごとにロードされる（16 + 2 * 16 + 4 * 64）バイトで乗算すると、システムの帯域幅はほぼ正確に〜5333MB/sになります。 。SSEなしでその帯域幅を飽和させるのは非常に簡単なはずです。a完全にキャッシュされたと仮定しても、読み取りbとc1回の反復で12バイトのデータが移動し、ナイーブはパイプライン処理を使用して3サイクルごとに新しい反復を開始できます。
完全なキャッシュに満たないものを想定するとa、算術演算と命令数のボトルネックがさらに少なくなります。私のコードのスピードアップのほとんどがより少ないロードを発行することによってもたらされ、より多くのスペースがの過去のキャッシュミスを追跡および推測するために自由になっているとしても、私は驚かないbでしょcうa。
より広いハードウェアはより大きな違いを生むかもしれません。DDR3-1333の3つのチャネルを実行しているNehalemシステムは、メモリ帯域幅を飽和させるために1サイクルあたり3 * 10667 / 2.66=12.6バイトを移動する必要があります。キャッシュに収まる場合、単一のスレッドでは不可能ですaが、64バイトでは、ベクトルのラインキャッシュミスがすぐに加算されます。キャッシュで欠落しているループ内の4つの負荷のうちの1つだけで、必要な平均帯域幅が16バイトになります。サイクル。
于 2010-03-05T13:21:46.837 に答える
0
短い答えはありません。長い答えはい、しかし効率的ではありません。整列されていない負荷を実行するとペナルティが発生し、あらゆる種類の利点が無効になります。b [i]の連続するインデックスが整列していることを保証できない限り、ベクトル化後のパフォーマンスが低下する可能性があります。
インデックスが何であるかを事前に知っている場合は、明示的なインデックスを展開して指定するのが最善です。テンプレートの特殊化とコード生成を使用して、同様のことを行いました。あなたが興味を持っているなら、私は共有することができます
あなたのコメントに答えるには、基本的に配列に集中する必要があります。すぐに試すのが最も簡単なのは、ループを2倍にブロックし、ローとハイを別々にロードしてから、通常のようにmm *_pdを使用することです。擬似コード：
__m128d a, result;
for(i = 0; i < n; i +=2) {
((double*)(&a))[0] = A[B[i]];
((double*)(&a))[1] = A[B[i+1]];
// you may also load B using packed integer instruction
result = _mm_add_pd(result, _mm_mul_pd(a, (__m128d)(C[i])));
}
関数名を正確に覚えていないので、再確認することをお勧めします。また、エイリアシングの問題がないことがわかっている場合は、ポインタでrestrictキーワードを使用します。これにより、コンパイラをより積極的にすることができます。
于 2010-02-28T05:00:40.740 に答える
0
配列インデックスの二重間接参照のため、これはそのままではベクトル化されません。ダブルスを使用しているため、特に最近のほとんどのCPUには2つのFPUがあるため、SSEから得られるものはほとんどまたはまったくありません。
于 2010-02-28T09:14:16.050 に答える