CUDA9以降の新しいWarpShuffle命令の挙動について

Volta GPUからWARPの挙動が変わった（32単位で実行**されなくなった**）ことをうけて、WarpShuffle系の命令が変更され、新たにマスク変数が導入された。参考になる資料は色々あるのだが、ちょっとピンとこなかった部分もあるので実験してみた。結果を共有しておく。
結果に基づく推測が含まれていることもあり、間違いがあったら申し訳ない。


参考資料

* [Using CUDA Warp-Level Primitives | NVIDIA Developer Blog](https://devblogs.nvidia.com/using-cuda-warp-level-primitives/)
* [Programming Guide :: CUDA Toolkit Documentation](https://docs.nvidia.com/cuda/archive/8.0/cuda-c-programming-guide/index.html#warp-shuffle-functions) ※CUDA8のWarpShuffleの説明
* [Programming Guide :: CUDA Toolkit Documentation](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#warp-shuffle-functions) ※CUDA10.1のWarpShuffleの説明

実験に使ったGPUはPascal=GeForce GTX1080とVolta=TeslaV100。


# WarpShuffle命令の基本（旧WarpShuffle命令）
WarpShuffle命令は、本来は共有（参照）できないはずの他スレッド（ただし同じWarp内に限る）のローカル変数の値を参照するための命令。共有メモリ（SharedMemory、GlobalMemory）を使うよりも高速な実行が期待できる。
例えば従来（CUDA10.1でもまだ利用はできるが、関数が古いよとコンパイラに警告される）のWarpShuffle命令で`tmp=__shfl_down(tmp, 1, warpSize)`を実行すれば、ローカル変数tmpの値が一つずつ前にずれる。

```c++:shuffle1.cu
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  tmp = __shfl_down(tmp, 1, warpSize);
  B[tid] = tmp;
}


int main(int argc, char **argv)
{
  int i, N;
  float *A, *B;
  float *dA, *dB;

  N = 32;
  A = (float*)malloc(sizeof(float)*N);
  B = (float*)malloc(sizeof(float)*N);

  for(i=0;i<N;i++){
        A[i] = (float)(i+1);
        B[i] = 0.0f;
  }

  cudaMalloc((void**)&dA, sizeof(float)*N);
  cudaMalloc((void**)&dB, sizeof(float)*N);

  printf("A\n");
  for(i=0; i<N; i++){
        printf(" %2.0f", A[i]);
  }
  printf("\n");
  printf("B (before)\n");
  for(i=0; i<N; i++){
        printf(" %2.0f", B[i]);
  }
  printf("\n");

  cudaMemcpy(dA, A, sizeof(float)*N, cudaMemcpyHostToDevice);
  cudaDeviceSynchronize();
  gpukernel<<<1,32>>>(dB, dA);
  cudaDeviceSynchronize();
  cudaMemcpy(B, dB, sizeof(float)*N, cudaMemcpyDeviceToHost);

  printf("B (after)\n");
  for(i=0; i<N; i++){
        printf(" %2.0f", B[i]);
  }
  printf("\n");

  cudaFree(dA); cudaFree(dB);
  free(A); free(B);
  return 0;
}
```

```
# CUDA-10.1でコンパイル（PascalでもVoltaでも動くようにarch/codeオプションを指定）
$ nvcc -O3 -arch=compute_60 -code=sm_60,sm_70 -o shuffle1 shuffle1.cu
shuffle1.cu(10): warning: function "__shfl_down(float, unsigned int, int)"
/mnt/nfs/packages/x86_64/cuda/cuda-10.1/bin/../targets/x86_64-linux/include/sm_30_intrinsics.hpp(207): here was declared deprecated ("__shfl_down() is deprecated in favor of __shfl_down_sync() and may be removed in a future release (Use -Wno-deprecated-declarations to suppress this warning).")

ptxas /tmp/tmpxft_0000fa84_00000000-5_shuffle1.ptx, line 40; warning : Instruction 'shfl' without '.sync' is deprecated since PTX ISA version 6.0 and will be discontinued in a future PTX ISA version
```
```
# CUDA-10.1のパスをLD_LIBRARY_PATH参照した状態で実行
$ ./shuffle1
A
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
B (before)
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
B (after)
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 32
```
配列Aのi番目にあった値が配列Bのi-1番目に入っており、「1つずつ前にずれた」ことがわかる。スレッド番号31については、値が更新されなかった。（一周回ってスレッド番号0から持ってきたりはせず、変わらぬまま。これはProgramming Guideに掲載されている通りの挙動である。）
もちろんPascalでもVoltaでも同じ結果である。


# 新しいWarpShuffle命令への置き換え
新しいWarpShuffle命令を使って書き換えるのは簡単で、新しく増えた第1変数（mask）に0xffffffffを指定すれば良い。
つまり、先のshuffle1.cuにおいて、```tmp = __shfl_down(tmp, 1, warpSize);```を```tmp = __shfl_down_sync(0xffffffff, tmp, 1, warpSize);```に変更すれば良い。GPUカーネルだけ書けばこのようになる。実行結果は、もちろんPascalでもVoltaでも、shuffle1.cuと一緒である。

```c++:shuffle2.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  tmp = __shfl_down_sync(0xffffffff, tmp, 1, warpSize);
  B[tid] = tmp;
}
```


# maskは何を意味するのか？：その1（基本）
CUDA 10.1のProgramming Guideを読む限り、maskに対応するスレッドのみが対応するWarpShuffle命令を実行するようだ。
そこで、maskに異なる値を指定して実行してみた。

```c++:shuffle3.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  tmp = __shfl_down_sync(0, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle4.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  tmp = __shfl_down_sync(0x00ff00ff, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

結果は、なんといずれもshuffle1.cuとまったく変わらなかった。
maskとはいったいなんなのか？
Programming Guideには「maskに対応するスレッドが走っていない場合に結果が不定」らしいことは書かれているが、納得できない結果である。


# maskは何を意味するのか？：その2（mask指定に問題がありそうな分岐）
流石にmaskが意味を成さないというのは気持ちが悪いため、あえてWarp内での分岐を含む意地悪な……というか、正直言って問題があるコードで実験してみた。

```c++:shuffle5.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid<16)tmp = __shfl_down_sync(0xffffffff, tmp, 1, warpSize);
  if(tid>=16)tmp = __shfl_down_sync(0xffffffff, tmp, 2, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle6.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid>=16)tmp = __shfl_down_sync(0xffffffff, tmp, 2, warpSize);
  if(tid<16)tmp = __shfl_down_sync(0xffffffff, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle7.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid<16)tmp = __shfl_down_sync(0, tmp, 1, warpSize);
  if(tid>=16)tmp = __shfl_down_sync(0, tmp, 2, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle8.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid>=16)tmp = __shfl_down_sync(0, tmp, 2, warpSize);
  if(tid<16)tmp = __shfl_down_sync(0, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle9.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid<16)tmp = __shfl_up_sync(0xffffffff, tmp, 1, warpSize);
  if(tid>=16)tmp = __shfl_up_sync(0xffffffff, tmp, 2, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle10.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid>=16)tmp = __shfl_up_sync(0xffffffff, tmp, 2, warpSize);
  if(tid<16)tmp = __shfl_up_sync(0xffffffff, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

shuffle5はWarpの前半と後半で分岐してmask=0xffffffffでシャッフル。前半はdelta=down幅1、後半はdelta=down幅2。shuffle6はshuffle5の順番入れ替え版。
shuffle7はWarpの前半と後半で分岐してmask=0でシャッフル。前半はdelta=down幅1、後半はdelta=down幅2。shuffle8はshuffle7の順番入れ替え版。
shuffle9とshuffle10はshuffle5とshuffle6のdownをupに変更したもの。
**いずれも問題がありそうなmaskの使い方をしているコードである。**
実行結果出力から「B (after)」だけ抽出すると以下の通り。**ただし途中の縦棒(|)は分岐の切れ目をわかりやすくするためにあとで追加したものなので注意されたし。**

```
A
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 | 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Pascal
$ ./shuffle5
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle6
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle7
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle8
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle9
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 0  0 17 18 19 20 21 22 23 24 25 26 27 28 29 30
$ ./shuffle10
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 0  0 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Volta
$ ./shuffle5
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle6
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle7
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 | 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
$ ./shuffle8
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 | 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
$ ./shuffle9
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
$ ./shuffle10
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
```

実行結果に違いが生じた。maskは仕事をしていないわけではなかったようだ。

Pascalの結果は一見maskの影響を受けていないように見えるが、分岐を跨ぐ部分（downのtid=15やupのtid=16,17）の結果が常に0になってしまっている。PascalはWarp単位でスケジューリングされるため、分岐があっても同じShuffle命令を触っており、maskがどうなっていようがシャッフル演算自体は実行されている、が、分岐方向が異なっているスレッドからは値をもらえない……みたいな動きをしているのだろうか？Warp外を参照してしまう場合とは挙動が異なる点は注意が必要かもしれない。

Voltaの結果はmaskの影響を受けているようで、mask=0では全ての結果が0になってしまっている。一方でPascalと異なり分岐を跨ぐ部分の結果が0になってしまっていない。これはこれで、よく考えてみると細かい実行タイミングのズレによってはシャッフル演算実行後のtmpの値をもらってきたりしないのか気になる。（例えばshuffle5のtid=15は17ではなく19になってしまわないのか？）

やや気になる点がいくつかあったが、PascalにしてもVoltaにしても、多少Warp内の分岐があったとしても「ある程度それらしい」シャッフル演算結果が得られる（エラー終了したりはしない）ようであることがわかった。


# maskは何を意味するのか？：その3（正しいと思われる分岐）
maskの正しい使い方を考えてコードを更新してみた。分岐対象のスレッドと同じmaskのみが指定されているコードである。

```c++:shuffle25.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid<16)tmp = __shfl_down_sync(0x0000ffff, tmp, 1, warpSize);
  if(tid>=16)tmp = __shfl_down_sync(0xffff0000, tmp, 2, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle26.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid>=16)tmp = __shfl_down_sync(0xffff0000, tmp, 2, warpSize);
  if(tid<16)tmp = __shfl_down_sync(0x0000ffff, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle27.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid<16)tmp = __shfl_down_sync(0xffff0000, tmp, 1, warpSize);
  if(tid>=16)tmp = __shfl_down_sync(0x0000ffff, tmp, 2, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle28.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid>=16)tmp = __shfl_down_sync(0x0000ffff, tmp, 2, warpSize);
  if(tid<16)tmp = __shfl_down_sync(0xffff0000, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle29.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid<16)tmp = __shfl_up_sync(0x0000ffff, tmp, 1, warpSize);
  if(tid>=16)tmp = __shfl_up_sync(0xffff0000, tmp, 2, warpSize);
  B[tid] = tmp;
}
```

```c++:shuffle30.cu
__global__ void gpukernel(float *B, float *A)
{
  int tid = threadIdx.x;
  float tmp = A[tid];
  if(tid>=16)tmp = __shfl_up_sync(0xffff0000, tmp, 2, warpSize);
  if(tid<16)tmp = __shfl_up_sync(0x0000ffff, tmp, 1, warpSize);
  B[tid] = tmp;
}
```

25,26,29,30は、それぞれ5,6,9,10に対して適切なmask、つまりtid<16には0x0000ffff、それ以外には0xffff0000を適用したもの。27,28は25,26のmaskをひっくり返したもの、つまりmaskが間違っている例である。

実行結果は次の通り。

```
A
  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 | 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

Pascal
$ ./shuffle25
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle26
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle27
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle28
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle29
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 0  0 17 18 19 20 21 22 23 24 25 26 27 28 29 30
$ ./shuffle30
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 0  0 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Volta
$ ./shuffle25
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle26
  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16  0 | 19 20 21 22 23 24 25 26 27 28 29 30 31 32 31 32
$ ./shuffle27
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 | 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
$ ./shuffle28
  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0 | 0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
$ ./shuffle29
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 0  0 17 18 19 20 21 22 23 24 25 26 27 28 29 30
$ ./shuffle30
  1  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 | 0  0 17 18 19 20 21 22 23 24 25 26 27 28 29 30
```

25,26,29,30の結果はPascalもVoltaも同様となった。
Warpをはみ出す範囲については値は変わらず、分岐を跨ぐ部分は0になっている。
27と28についてはPascalとVoltaで異なる結果となった。Voltaではmaskがあわない部分については演算を行っていないような結果であるのに対して、Pascalはmaskを見ないで演算しているような結果である。

**こうしてみると、やはりshuffle3（mask=0のとき）のVoltaの挙動だけが大きく外れてしまっている気がする。実はバグなのだろうか……？**

最後に、Voltaで```-arch=compute_60```を指定し、Pascalのように1Warp=32threads単位で動くモードにした場合の動作も確認してみたのだが、特に結果は変わらなかった。
もしかしたら一部Pascalと同じ結果になるのだろうかと思ったのだが、そんなことはなかったぜ……。


# 余談：それで、性能に寄与する話なの？
中の人はHPC民につき、元々はVolta専用のコードを書くことで性能が上がる余地を考えていた。
Voltaで実行したいプログラムの中に、32スレッドではなく16スレッドとか8スレッドとかが1グループとなって動いて結果をシャッフルで集約しているカーネルがあったのでどうかと思ったのだが、同一WARP内の各グループは基本的に同じ長さのループ（ズレるとしてもたかだか1）なので、意味がなさそうだ。（もしかしたら、適切なタイミングでsyncを行う、とかで性能向上する余地がある……かも？）
この場合は「Cooperative Groups」というヤツを活用するべきかもしれない。これはこれで、思っているような効果があるものなのかどうかがピンときてないんだけどね……。