XeonPhiでOpenCLやってみる（ベクトル命令出す編

この記事は、CUDA & OpenCL Advent Calendar 2014の10日目の記事です

## 前回のあらすじ
XeonPhiでのOpenCLの出力に成功した我々。
さて、では次なる課題として、OpenCLを使うと、XeonPhiの特徴ともいえる、512bit SIMD命令は出力できるのか、という点に着目してみようと思います

## OpenCL カーネルコードのコンパイル
昨日のカレンダーの記事のネタがそのままなので、すこし心苦しい感じもなきにしもあらずではあるけれど。
OpenCLデバイス上で実行されるコードを、OpenCLの世界ではOpenCLカーネルと呼称します。
カーネルはホストから呼ばれるエントリポイントとなる関数と、それが呼び出す関数群として定義されます。
エントリポイントとなるカーネルについては、__kernel修飾子を頭につけることで定義することができます。

さて、ではこのカーネルをOpenCLデバイス上のバイナリとしてコンパイルするには、どのようにすればよいのでしょうか？
方法は、オンラインコンパイルとオフラインコンパイルの二つがあります。

### オンラインコンパイル
OpenCLは専用のAPIとして、プログラム実行時にカーネルコードのコンパイルを可能とするものを用意しています。
昨日のカレンダーの記事はオンラインコンパイル、いわゆるJITコンパイラを使用する際にどのようにしてカーネルをプログラム中に保持するか、という問題についての一つの解です。
実行時にビルドするという特性上、ユーザーのプラットフォームに依らないという利点があります。ただし、実行時にビルドしなきゃいけないため、初回起動時に時間がかかります
(一応、各ベンダはコンパイルしたものをキャッシュする仕組みを用意はしてますが)

### オフラインコンパイル
ふつうのプログラムと同様、事前にビルドしてバイナリになっているものをロードして使う方法です。
想定する各アーキテクチャ、たとえばx86のCPUだったり、あるいはIris proだったり、はたまたNVIDIAのGPUだったり、AMDのGPUだったり…と、多種多様な環境を想定しなければなりませんが、コンパイル時間はかかることなく実行が可能です。
あと、文字列としてプログラム中に埋め込んだり、.clファイルで外部に置いたりしないため、カーネルコードの秘匿性が高まるという利点もあります。

## ベクトル命令出す編
ここから本編。
XeonPhiのベクトル命令、IMCI命令ですが、OpenCLの言語仕様的には、そういったベクトル型はどのように定義されているのでしょう。
OpenCLの世界では、__m128等に代表されるベクトル変数は、float4, float8...といった形でマッピングされます。
つまり、それらの型を使ってあげれば、コンパイラがよろしく判断してくれてうれしいことになる……のでは！？

さて、ここで一つ問題が。
オンラインコンパイルでどうやってアセンブリ見るのよ？
なんかメモリを吐けば確かできた気がしますが、そういう黒魔法は黒魔道士に任せて、まっとうな方法でいきましょう。

IntelさんはKernel builderというOpenCLカーネルのオフラインコンパイラ(という名の開発環境っぽくなってしまってもいる)を用意してくれています。
これを使いましょう。
使うカーネルはこちら
https://bitbucket.org/telmin/opencl_vec

まずはコマンドのヘルプを観ましょう

```
$ ioc64
sh: line 1:  2475 アボートしました/opt/intel/ism/bin/intel64/intelremotedialog > /dev/null 2>&1
No input parameters
usage: ioc64 [COMMAND] <ARGUMENT> [OPTIONS]
Kernel Builder for OpenCL API - compiler command line, version 1.5.0.92
Copyright (C) 2014 Intel Corporation.  All rights reserved.

COMMAND:
    -cmd=<command>                   - Command to be performed:
                                      'build' create executable IR from source
                                       code (default if none specified)
                                      'compile' create compiled object IR from
                                       source code
                                      'link' create executable IR / library
                                       from object IR and libraries
ARGUMENTS:
    -input=<input_file_path>         - Build the OpenCL code given in
                                       <input_file_path> (use with the 'build'
                                       & 'compile' commands)
    -binary="<binary_files_paths>"   - Link binary files, comma seperated if
                                       more than one (use with 'link' command)
    -version                         - Show tool version
    -help                            - Show available commands
OPTIONS:
    -device=<device_type>            - Set target device type:
                                       'cpu' for Intel CPU device (default)
                                       'gpu' for Intel(R) Graphics device
                                       'co' for Intel(R) Xeon Phi(TM) coprocessor device
    -targetos=<os>                   - Set target operating system if it is
                                       different from current:
                                       (supported in 32-bit version only)
                                       'android' (use with 'cpu' device only)
    -simd=<instruction_set_arch>     - Set target instruction set architecture
                                       (use with 'cpu' device only):
                                       'sse42' for Streaming SIMD Extension 4.2
                                       'avx' for Advanced Vector Extensions
                                       'avx2' for Advanced Vector Extensions 2
    -output[=<output_file_path>]     - Write the build log to <output_file_path>
    -asm[=<file_path>]               - Generate assembly code for Intel CPU
                                       & Intel(R) Xeon Phi(TM) coprocessor
    -llvm[=<file_path>]              - Generate llvm code
    -llvm-spir32[=<file_path>]       - Generate llvm spir code (32-bit)
    -llvm-spir64[=<file_path>]       - Generate llvm spir code (64-bit)
    -ir[=<file_path>]                - Generate intermediate binary file
    -spir32[=<file_path>]            - Generate SPIR (32-bit) binary file
    -spir64[=<file_path>]            - Generate SPIR (64-bit) binary file
    -bo[="<build_options>"]          - Add build options

```

多分、inputオプションでファイルを指定して、deviceオプションにcoを指定してあげればXeonPhi向けのコードがコンパイルされると思うんですよ

```
$ ioc64 -input=kernel.cl -device=co -asm=kernel.asm
sh: line 1:  2520 アボートしました/opt/intel/ism/bin/intel64/intelremotedialog > /dev/null 2>&1
No command specified, using 'build' as default
OpenCL Intel(R) Xeon Phi(TM) Coprocessor device was found!
Device name: Intel(R) Many Integrated Core Acceleration Card
Device version: OpenCL 1.2 (Build 8)
Device vendor: Intel(R) Corporation
Device profile: FULL_PROFILE
Compilation started
/home/telmin/Source/opencl_vec/kernel.cl:20:7: warning: expression result unused
/home/telmin/Source/opencl_vec/kernel.cl:44:7: warning: expression result unused
/home/telmin/Source/opencl_vec/kernel.cl:68:7: warning: expression result unused
Compilation done
Linking started
Linking done
Device build started
Device build done
Build started
Kernel <add_scalar> was successfully vectorized (16)
Kernel <add_vec> was successfully vectorized (16)
Kernel <mul_scalar> was successfully vectorized (16)
Kernel <mul_vec> was successfully vectorized (16)
Kernel <fma_scalar> was successfully vectorized (16)
Kernel <fma_vec> was successfully vectorized (16)
Kernel <gather_vec> was successfully vectorized (16)
Done.
Build succeeded!
```

やったーコンパイルされたよー
…なんかスゲー不穏な文字が見えましたが、まぁいいとして、一つずつカーネルを見ていきましょう。

### add

```cpp
__kernel void add_scalar(__global float* a, __global float* b, __global float* c, const unsigned int num)
{
  for(unsigned int i = 0; i < num; ++i) {
    c[i] = a[i] + b[i];
  }
}

__kernel void add_vec(__global float* a, __global float* b, __global float* c, const unsigned int num)
{
  unsigned int i;
  __global float16* v_a = (__global float16*)a;
  __global float16* v_b = (__global float16*)b;
  __global float16* v_c = (__global float16*)c;

  for(i = 0; i < num / 16; ++i) {
    v_c[i] = v_a[i] + v_b[i];
  }

  i *= 16;
  for(i; i < num; ++i) {
    c[i] = a[i] + b[i];
  }
}
```

こんなソースコードをコンパイルすると、何が出るかっていうと

```
# mark_begin;
# Threads 4
add_scalar:
..B1.1:                         # Preds ..B1.0 Latency 9
        movq      %rsi, %rdi                                    # c1
        movq      %rax, %r8                                     # c1
        vbroadcastss (%rdi), %zmm0{%k1}                         # c5
        movq      %rdx, %r9                                     # c5
        vaddps    (%r8){1to16}, %zmm0, %zmm1{%k1}               # c9
        lea       4(%rdi), %rsi                                 # c9
        lea       4(%r8), %rax                                  # c13
        lea       4(%r9), %rdx                                  # c13
        vpackstorelps %zmm1, (%r9){%k1}                         # c17
        decl      %ecx                                          # c17
        jne       ..B1.3        # Prob 0%                       # c21
# mark_end;

# mark_begin;
# Threads 4
add_vec:
..CL4:
..B2.3:                         # Preds ..B2.1 ..B2.3 Latency 29
        movq      %r9, %r10                                     # c1
        incl      %edi                                          # c1
        vmovaps   (%r10,%rax), %zmm0                            # c5
        vprefetch0 64(%r10,%rax)                                # c5
        vaddps    (%r10,%rdx), %zmm0, %zmm1                     # c9
        vprefetch1 320(%r10,%rax)                               # c9
        vprefetch0 64(%r10,%rdx)                                # c13
        lea       64(%r10), %r9                                 # c13
        vprefetch1 320(%r10,%rdx)                               # c17
        cmpl      %r8d, %edi                                    # c17
        vprefetche0 64(%r10,%rcx)                               # c21
        movb      %al, %al                                      # c21
        vprefetche1 320(%r10,%rcx)                              # c25
        movb      %al, %al                                      # c25
        vmovdqa32 %zmm1, (%r10,%rcx)                            # c29
        jb        ..B2.3        # Prob 0%                       # c29
                                # LOE rax rdx rcx rbx rbp r9 r12 r13 r14 r15 esi edi r8d
...
..B2.9:                         # Preds ..B2.9 ..B2.8 Latency 21
        movq      %rax, %rcx                                    # c1
        movq      %r8, %rdi                                     # c1
        vbroadcastss (%rcx), %zmm0{%k1}                         # c5
        movq      %rdx, %r9                                     # c5
        vaddps    (%rdi){1to16}, %zmm0, %zmm1{%k1}              # c9
        lea       4(%rcx), %rax                                 # c9
        lea       4(%rdi), %r8                                  # c13
        lea       4(%r9), %rdx                                  # c13
        vpackstorelps %zmm1, (%r9){%k1}                         # c17
        decl      %esi                                          # c17
        jne       ..B2.9        # Prob 0%                       # c21
                                # LOE rax rdx rbx rbp r8 r12 r13 r14 r15 esi k1
..CL6:
# mark_end;
```

x86asmわからん(震え声
かろうじてわかるのは、どっちもvaddpsが出ているってことですね。
しばらく眺めてみたけど%r10が何したいのかさっぱりわからん…お前何ものだよ…

### fma

さて、XeonPhiの演算能力を支える基礎的なものの一つとして、fmaの存在があげられます
このfmaははたして出るのでしょうか？

```cpp
__kernel void fma_scalar(__global float* a, __global float* b, __global float* c, const unsigned int num)
{
  for(unsigned int i = 0; i < num; ++i) {
    c[i] = a[i] * b[i] + c[i];
  }
}

__kernel void fma_vec(__global float* a, __global float* b, __global float* c, const unsigned int num)
{
  unsigned int i;
  __global float16* v_a = (__global float16*)a;
  __global float16* v_b = (__global float16*)b;
  __global float16* v_c = (__global float16*)c;

  for(i = 0; i < num / 16; ++i) {
    v_c[i] = v_a[i] * v_b[i] + v_c[i];
  }

  i *= 16;
  for(i; i < num; ++i) {
    c[i] = a[i] * b[i] + c[i];
  }
}
```

これをコンパイルすると

```
# mark_begin;
# Threads 4
fma_scalar:
..CL18:
..B5.3:                         # Preds ..B5.3 ..B5.2 Latency 21
        movq      %rsi, %rdi                                    # c1
        movq      %rax, %r8                                     # c1
        vbroadcastss (%rdi), %zmm1{%k1}                         # c5
        movq      %rdx, %r9                                     # c5
        vbroadcastss (%r8), %zmm0{%k1}                          # c9
        lea       4(%rdi), %rsi                                 # c9
        vfmadd213ps (%r9){1to16}, %zmm0, %zmm1{%k1}             # c13
        lea       4(%r8), %rax                                  # c13
        lea       4(%r9), %rdx                                  # c17
        decl      %ecx                                          # c17
        vpackstorelps %zmm1, (%r9){%k1}                         # c21
        jne       ..B5.3        # Prob 0%                       # c21
                                # LOE rax rdx rbx rbp rsi r12 r13 r14 r15 ecx k1
# mark_end;

# mark_begin;
# Threads 4
        .align    16,0x90
        .globl fma_vec
fma_vec:
..CL20:
..B6.3:                         # Preds ..B6.1 ..B6.3 Latency 29
        movq      %r9, %r10                                     # c1
        incl      %edi                                          # c1
        vmovaps   (%r10,%rax), %zmm1                            # c5
        vprefetch0 64(%r10,%rax)                                # c5
        vmovaps   (%r10,%rdx), %zmm0                            # c9
        vprefetch1 320(%r10,%rax)                               # c9
        vfmadd213ps (%r10,%rcx), %zmm0, %zmm1                   # c13
        vprefetch0 64(%r10,%rdx)                                # c13
        vprefetch1 320(%r10,%rdx)                               # c17
        lea       64(%r10), %r9                                 # c17
        vprefetch0 64(%r10,%rcx)                                # c21
        cmpl      %r8d, %edi                                    # c21
        vprefetch1 320(%r10,%rcx)                               # c25
        movb      %al, %al                                      # c25
        vmovdqa32 %zmm1, (%r10,%rcx)                            # c29
        jb        ..B6.3        # Prob 0%                       # c29
                                # LOE rax rdx rcx rbx rbp r9 r12 r13 r14 r15 esi edi r8d

# mark_end;
```

おぉ、なんかやっぱ出てるっぽい。
そしてやっぱりx86asmわからん…


### gather
XeonPhiといえばGather！
Haswellでは次回作にご期待くださいとか言われちゃってるGather！！！
あるとすっごい便利なんだけど速くなくていつもかなしくなるGather！！
これ、OpenCLだとどうやって書けばいいんだろう…
苦肉の策が以下

```cpp
__kernel void gather_vec(__global float* a, __global float* b)
{
  unsigned int index[16] = {15, 14, 13, 12, 11, 10, 9, 8,
                            7,  6,  5,  4, 3, 2, 1, 0};

  float16 vec;
  vec.s0 = a[index[0]];
  vec.s1 = a[index[1]];
  vec.s2 = a[index[2]];
  vec.s3 = a[index[3]];
  vec.s4 = a[index[4]];
  vec.s5 = a[index[5]];
  vec.s6 = a[index[6]];
  vec.s7 = a[index[7]];
  vec.s8 = a[index[8]];
  vec.s9 = a[index[9]];
  vec.sa = a[index[10]];
  vec.sb = a[index[11]];
  vec.sc = a[index[12]];
  vec.sd = a[index[13]];
  vec.se = a[index[14]];
  vec.sf = a[index[15]];

  __global float16* v_b = (__global float16*)b;
  *v_b = vec;
}
```

で、これをコンパイルすると

```
# mark_begin;
# Threads 4
        .align    16,0x90
        .globl gather_vec
gather_vec:
..B7.1:                         # Preds ..B7.0 Latency 137
        pushq     %rbp                                          #
        movq      %rsp, %rbp                                    #
        andq      $-64, %rsp                                    #
        pushq     %r12                                          # c1
        pushq     %r13                                          # c5
        pushq     %r14                                          # c9
        pushq     %r15                                          # c13
        pushq     %rbx                                          # c17
        movq      (%rdi), %r15                                  # c21
        movl      (%r15), %edx                                  # c25
        movl      32(%r15), %r8d                                # c25
        movl      4(%r15), %ecx                                 # c29
        movl      %edx, -16(%rsp)                               # c29
        movl      %ecx, -8(%rsp)                                # c33
        movq      8(%rdi), %rax                                 # c33
        movl      40(%r15), %esi                                # c37
        movq      %rax, -24(%rsp)                               # c37
        movl      %esi, -68(%rsp)                               # c41
        movl      %r8d, -60(%rsp)                               # c41
        movl      8(%r15), %r14d                                # c45
        movl      16(%r15), %r12d                               # c45
        movl      12(%r15), %r13d                               # c49
        movl      20(%r15), %r11d                               # c49
        movl      24(%r15), %r10d                               # c53
        movl      36(%r15), %edi                                # c53
        movl      28(%r15), %r9d                                # c57
        movl      44(%r15), %ebx                                # c57
...

```

延々とグルグルまわしてはセットするコードが出来上がりましたとさ
うーん、OpenCLでGatherを吐くコードはいまだ見つからず…

## まとめ
さて、OpenCLでベクトル命令が吐けるかどうかを検証してまいりましたが、簡単な命令であればあっさり吐いてくれそうですね！
ていうか、コンパイラのログを注視するとですね

```
Kernel <add_scalar> was successfully vectorized (16)
Kernel <add_vec> was successfully vectorized (16)
Kernel <mul_scalar> was successfully vectorized (16)
Kernel <mul_vec> was successfully vectorized (16)
Kernel <fma_scalar> was successfully vectorized (16)
Kernel <fma_vec> was successfully vectorized (16)
Kernel <gather_vec> was successfully vectorized (16)
```

とか言ってくれちゃってて、自動ベクトル化がめっちゃ働いている感じがするわけです。
しかしながら、自動ベクトル化したスカラコードを見ていると、プリフェッチ命令が出力されていないようにみえます。
XeonPhiはメモリアクセスのペナルティが本当に大きいため、どうにかしてデータをキャッシュに引っ張ってくる必要があります。ICCでコードを書いていると、本当にベクトルロードをするたびにプリフェッチ命令が出力されるという事態に遭遇します。
また、実際問題、キャッシュアクセスを意識するようにしただけで速度がガラッと変わるのはよくある話です。
そんなわけで、やっぱりちゃんとコンパイラに優しいコードを書くためにも、float16を多用してあげることがXeonPhiと仲良くなる第一歩ということです！

さて、じゃあそんなXeonPhiのOpenCL性能は？ というのは、まて次回ということで。

明日は@uchihashi_kさんです！