HPCメモ

HPC(High Performance Computing)に関連したりしなかったりすることのメモ書き

最適化の限界

昨日ついったーで面白い実験をみかけたので、ちょっと考察してみました。

lpha-z.hatenablog.com

元記事の方は、JuliaとCで同じ処理を行うプログラム書いて性能を比較されてたんですが、 私としてはJuliaとの比較ではなくHWスペックに対してC言語の実験結果が異様に遅いことに興味をひかれました。

使われていたプログラムは、1000要素の倍精度の配列の値を10000000足し合わせるというもので、 Intel(R) Xeon(R) CPU E5-2603 v3 @ 1.60GHzで測定されたとのことです。

とりあえず、メモリネックなのか演算ネックなのか確認したいので、HWのB/Fを計算してみましょう。

こちらにあるようにE5-2603 V3はHaswell世代の6コアのCPUですが、今回は1コアしか使ってないようなのでHW的には1コア分の性能だけを見ます。

ark.intel.com

L1Dキャッシュは、32KBなので、1000要素の配列なら全データがオンキャッシュで動作します。初回アクセス時にキャッシュミスは発生しますが、 1千万回も繰り返しているので性能に与える影響はほぼ無視できるはずです。Haswel世代のL1Dキャッシュは、load時は64Bytes/cycleとのことなので、この値をメモリアクセス性能とします*1

演算性能は、1コアで256bit幅のFMA2基ですが今回は加算しか演算が無いのでただの加算器が2つと考えて、8Flop/cycleとします。

したがって、この場合のHWのB/F値は8となります。

一方、プログラム側のB/F値は8Byteロードに対して1演算なので8B/Fとなります。

ということは、このプログラムは1コア内の加算器を限界まで使えるので、理論上の性能の限界は2x4x1.6 = 12.8GFlops となります。 演算数から時間を逆算すると、1000x10000000/(12.8x10003) = 0.78 sec ですね。

こちらのツイートで示されている、Juliaでまっとうな書き方をした時の性能が0.89secということなので、もし同じようなCPUでの測定結果であれば 少なくともこのケースに限ってはJuliaの最適化は優秀と言えそうです。

さて、C言語の方ですが、元記事ではgccで-O3でコンパイルして18.8secという結果だったそうです。

アセンブリも示されているのでまずそちらを見ると、加算命令としてaddsd命令が使われていることがわかります。 つまり、SIMD演算を使っておらず8Byteづつ加算している状態です。 また、twitterで元の記事を書かれた方から指摘をいただいて気付いたんですが、配列の先頭から順に加算するため、 せっかく2ポートで同時に加算命令が発行できるのにこの機能も使えません。 したがって、さきほど計算した理論性能と比較すると、SIMDが使えないことで1/4さらに命令発行が1命令/cycleとなることで1/2になり合わせて1/8の性能しか使えません。

そうすると実行時間で考えるとだいた6.4秒くらいですね。 実測値は、18.8秒なのでさらに3倍近く遅いのですが、おそらくこれは、ADDSD命令のlatencyとThroughputの差と思われます。

これまで1cycle毎に命令が行なわれる*2前提で考えてきたので、1演算あたり1cycleとしていましたが、先行する加算命令の結果を書かれたレジスタから値を読んで後続の加算命令が行なわれているので、ここはlatencyである3cycleが1演算あたり必要な状態と考えるべきです。*3 したがって、さらに3倍遅くて19秒くらい(なんかちょっと多いな・・・)という計算で辻褄が合いそうです。

最適化

さて、ここまで、性能低下の要因を色々とみてきましたが原因は全部同じで、「1要素づつしか加算していない」ということにつきます。*4 じゃ、どうすりゃ良いのかというと

for(i =0; i<1000; i++){
 s += a[i];
}

というような最内ループをストリップマイニングして

for(i=0; i<1000; i+=8){
 s0 += a[i];
 s1 += a[i+1];
 s2 += a[i+2];
 s3 += a[i+3];
 s4 += a[i+4];
 s5 += a[i+5];
 s6 += a[i+6];
 s7 += a[i+7];
}
s=s0+s1+s2+s3+s4+s5+s6+s7

という形に書き換えるようにコンパイラに指示すれば良いわけです。 gccの場合はその名もずばりの-floop-strip-mineというのがあるので、これと(SIMD化のために)-march=haswellを指定してコンパイルすればたぶん最適化してくれるでしょう。

実際に(全体をコンパイルすると読むのが面倒なので)次のような2重ループの関数を作って、-Sでアセンブリを吐いてみます。

#include <stddef.h>
double sum(double* a){
  double s=0;
  for(size_t j=0; j<10000000;j++){
    for(size_t i =0; i<1000; i++){
       s += a[i];
    }
  }
  return s;
}

-O3だけの時

> docker run -it --rm -w /usr/src -v ${PWD}:/usr/src gcc:9.1.0 gcc -S -O3 test.c -o test_O3.s
> cat test_O3.s
        .file   "test.c"
        .text
        .p2align 4
        .globl  sum
        .type   sum, @function
sum:
.LFB0:
        .cfi_startproc
        movq    %rdi, %rsi
        movl    $10000000, %ecx
        pxor    %xmm0, %xmm0
        leaq    8000(%rdi), %rdx
.L2:
        movq    %rsi, %rax
        .p2align 4,,10
        .p2align 3
.L3:
        movsd   (%rax), %xmm1
        addq    $16, %rax
        addsd   %xmm0, %xmm1
        movsd   -8(%rax), %xmm0
        addsd   %xmm1, %xmm0
        cmpq    %rax, %rdx
        jne     .L3
        subq    $1, %rcx
        jne     .L2
        ret
        .cfi_endproc
.LFE0:
        .size   sum, .-sum
        .ident  "GCC: (GNU) 9.1.0"
        .section        .note.GNU-stack,"",@progbits

-march=haswellと-floop-strip-mineをつけると

> docker run -it --rm -w /usr/src -v ${PWD}:/usr/src gcc:9.1.0  gcc -S -O3 -march=haswell -floop-strip-mine test.c -o test_O3_haswell_stripmine.s
> cat test_O3_haswell_stripmine.s
        .file   "test.c"
        .text
        .section        .text.unlikely,"ax",@progbits
.LCOLDB1:
        .text
.LHOTB1:
        .p2align 4
        .globl  sum
        .type   sum, @function
sum:
.LFB0:
        .cfi_startproc
        movq    %rdi, %rax
        vxorpd  %xmm1, %xmm1, %xmm1
        leaq    8000(%rdi), %rsi
        xorl    %ecx, %ecx
.L3:
        movq    %rax, %rdx
.L2:
        vmovupd (%rdx), %xmm2
        vmovupd (%rdx), %ymm3
        addq    $32, %rdx
        vaddsd  %xmm1, %xmm2, %xmm1
        vunpckhpd       %xmm2, %xmm2, %xmm2
        vextractf128    $0x1, %ymm3, %xmm0
        vaddsd  %xmm2, %xmm1, %xmm1
        vaddsd  %xmm0, %xmm1, %xmm1
        vunpckhpd       %xmm0, %xmm0, %xmm0
        vaddsd  %xmm0, %xmm1, %xmm1
        cmpq    %rdx, %rsi
        jne     .L2
        incl    %ecx
        cmpl    $10000000, %ecx
        je      .L6
        jmp     .L3
        .cfi_endproc
        .section        .text.unlikely
        .cfi_startproc
        .type   sum.cold, @function
sum.cold:
.LFSB0:
.L6:
        vmovapd %xmm1, %xmm0
        vzeroupper
        ret
        .cfi_endproc
.LFE0:
        .text
        .size   sum, .-sum
        .section        .text.unlikely
        .size   sum.cold, .-sum.cold
.LCOLDE1:
        .text
.LHOTE1:
        .ident  "GCC: (GNU) 9.1.0"
        .section        .note.GNU-stack,"",@progbits

(・3・)あるぇ～？

3オペランド方式のvバージョンになっていますが、あいかわらずaddsdが使われているようです。

ちなみに、icc 19でデフォルトオプションだとこんな感じになります。

# mark_description "-std=c11 -S";
    .file "test.c"
    .text
..TXTST0:
.L_2__routine_start_sum_0:
# -- Begin  sum
    .text
# mark_begin;
       .align    16,0x90
    .globl sum
# --- sum(double *)
sum:
# parameter 1: %rdi
..B1.1:                         # Preds ..B1.0
                                # Execution count [0.00e+00]
    .cfi_startproc
..___tag_value_sum.1:
..L2:
                                                          #2.22
        xorl      %edx, %edx                                    #5.5
        pxor      %xmm0, %xmm0                                  #3.11
        pxor      %xmm1, %xmm1                                  #3.11
        xorl      %eax, %eax                                    #4.3
                                # LOE rdx rbx rbp rdi r12 r13 r14 r15 eax xmm0 xmm1
..B1.2:                         # Preds ..B1.4 ..B1.1
                                # Execution count [6.25e+05]
        movaps    %xmm1, %xmm10                                 #3.11
        movl      %eax, %ecx                                    #4.3
        movaps    %xmm10, %xmm9                                 #3.11
        movaps    %xmm9, %xmm8                                  #3.11
        movaps    %xmm8, %xmm7                                  #3.11
        movaps    %xmm7, %xmm6                                  #3.11
        movaps    %xmm6, %xmm5                                  #3.11
        movsd     (%rdi,%rdx,8), %xmm2                          #6.13
        movaps    %xmm5, %xmm4                                  #3.11
        unpcklpd  %xmm2, %xmm2                                  #6.13
        movaps    %xmm4, %xmm3                                  #3.11
                                # LOE rdx rbx rbp rdi r12 r13 r14 r15 eax ecx xmm0 xmm1 xmm2 xmm3 xmm4 xmm5 xmm6 xmm7 xmm8 xmm9 xmm10
..B1.3:                         # Preds ..B1.3 ..B1.2
                                # Execution count [6.25e+08]
        addl      $16, %ecx                                     #4.3
        addpd     %xmm2, %xmm10                                 #6.8
        addpd     %xmm2, %xmm9                                  #6.8
        addpd     %xmm2, %xmm8                                  #6.8
        addpd     %xmm2, %xmm7                                  #6.8
        addpd     %xmm2, %xmm6                                  #6.8
        addpd     %xmm2, %xmm5                                  #6.8
        addpd     %xmm2, %xmm4                                  #6.8
        addpd     %xmm2, %xmm3                                  #6.8
        cmpl      $10000000, %ecx                               #4.3
        jb        ..B1.3        # Prob 99%                      #4.3
                                # LOE rdx rbx rbp rdi r12 r13 r14 r15 eax ecx xmm0 xmm1 xmm2 xmm3 xmm4 xmm5 xmm6 xmm7 xmm8 xmm9 xmm10
..B1.4:                         # Preds ..B1.3
                                # Execution count [6.25e+05]
        addpd     %xmm9, %xmm10                                 #3.11
        addpd     %xmm7, %xmm8                                  #3.11
        addpd     %xmm5, %xmm6                                  #3.11
        addpd     %xmm3, %xmm4                                  #3.11
        addpd     %xmm8, %xmm10                                 #3.11
        addpd     %xmm4, %xmm6                                  #3.11
        addpd     %xmm6, %xmm10                                 #3.11
        movaps    %xmm10, %xmm2                                 #3.11
        incq      %rdx                                          #5.5
        unpckhpd  %xmm10, %xmm2                                 #3.11
        addsd     %xmm2, %xmm10                                 #3.11
        addsd     %xmm10, %xmm0                                 #3.11
        cmpq      $1000, %rdx                                   #5.5
        jb        ..B1.2        # Prob 99%                      #5.5
                                # LOE rdx rbx rbp rdi r12 r13 r14 r15 eax xmm0 xmm1
..B1.5:                         # Preds ..B1.4
                                # Execution count [7.63e-02]
        ret                                                     #9.10
        .align    16,0x90
                                # LOE
    .cfi_endproc
# mark_end;
    .type   sum,@function
    .size   sum,.-sum
..LNsum.0:
    .data
# -- End  sum
    .data
    .section .note.GNU-stack, ""
# End

なんかめっちゃアンロールしてますね。 -qopt-report=5で出力される最適化レポートを見ると

ループの開始 test.c(5,5)
   リマーク #25444: ループの入れ子の交換: ( 1 2 ) --> ( 2 1 )
   リマーク #15542: ループ はベクトル化されませんでした: 内部ループがすでにベクトル化されています。   [ test.c(5,5) ]

   ループの開始 test.c(4,3)
      リマーク #15305: ベクトル化のサポート: ベクトル長 2
      リマーク #15399: ベクトル化のサポート: アンロールファクターが 8 に設定されます。
      リマーク #15309: ベクトル化のサポート: 正規化されたベクトル化のオーバーヘッド 0.525
      リマーク #15301: 変換されたループ がベクトル化されました。
      リマーク #15475: --- ベクトルのコストサマリー開始 ---
      リマーク #15476: スカラーのコスト: 6 
      リマーク #15477: ベクトルのコスト: 2.500 
      リマーク #15478: スピードアップの期待値: 2.390 
      リマーク #15488: --- ベクトルのコストサマリー終了 ---
      リマーク #25015: ループの最大トリップカウントの予測=625000
   ループの終了
ループの終了

あ、入れ替えられてる・・・そらそうかorz

*1:もちろんここではキャッシュにしかアクセスしてませんが要するにB/F値の分子という意味で捉えてください

*2:正確に言うとリタイアするですが

*3:要するにバブルでパイプラインがすっかすかの状態

*4:3つの壁のうちの1つ、命令レベルの並列性の壁ですね