cuda
並行縮小（例えば、配列を合計する方法）
サーチ…
備考
並列縮小アルゴリズムは、通常、要素の配列を結合して単一の結果を生成するアルゴリズムを指します。このカテゴリに属する典型的な問題は次のとおりです。
- 配列内のすべての要素を集計する
- 配列内で最大値を見つける
一般に、並列リダクションは任意のバイナリ連想演算子 、すなわち(A*B)*C = A*(B*C)適用することができます。このような演算子*を使用すると、パラレルリダクションアルゴリズムは配列の引数を繰り返しペアごとにグループ化します。各ペアは他のペアと並列に計算され、1ステップで全体のアレイサイズが半分になります。このプロセスは、単一の要素のみが存在するまで繰り返される。
演算子が結合的であることに加えて可換性 （すなわち、 A*B = B*A ）である場合、アルゴリズムは異なるパターンで対になることができる。理論的には違いはありませんが、実際にはより良いメモリアクセスパターンが得られます。
すべての連想演算子が可換性 - マトリックス乗算などを取るわけではありません。
交換可能演算子のための単一ブロック並列化
CUDAの並列削減の最も簡単なアプローチは、タスクを実行するために単一のブロックを割り当てることです。
static const int arraySize = 10000;
static const int blockSize = 1024;
__global__ void sumCommSingleBlock(const int *a, int *out) {
int idx = threadIdx.x;
int sum = 0;
for (int i = idx; i < arraySize; i += blockSize)
sum += a[i];
__shared__ int r[blockSize];
r[idx] = sum;
__syncthreads();
for (int size = blockSize/2; size>0; size/=2) { //uniform
if (idx<size)
r[idx] += r[idx+size];
__syncthreads();
}
if (idx == 0)
*out = r[0];
}
...
sumCommSingleBlock<<<1, blockSize>>>(dev_a, dev_out);
これは、データサイズがそれほど大きくない場合（約数千の要素）、最も実現可能です。これは、通常、削減がより大きなCUDAプログラムの一部である場合に発生します。入力がblockSize最初から一致した場合、最初のforループを完全に削除することができます。
最初のステップでは、スレッドよりも多くの要素がある場合、完全に独立して追加することに注意してください。問題がblockSizeに減少した場合にのみ、実際の並列減少がトリガされます。乗算、最小値、最大値などの他の任意の可換性連想演算子にも同じコードを適用できます。
アルゴリズムは、例えば、ワープレベルの並列縮小を使用することによって、より高速にすることができることに留意されたい。
非可換性演算子のための単一ブロック並列化
非可換的な演算子の並列的な削減は、可換的なバージョンと比べて少し複雑です。この例では単純化のために整数以上の加算を使用しています。これは、例えば、実際には非可換である行列乗算で置き換えることができます。このとき、0は乗算の中立要素、すなわち恒等行列に置き換えなければならないことに注意してください。
static const int arraySize = 1000000;
static const int blockSize = 1024;
__global__ void sumNoncommSingleBlock(const int *gArr, int *out) {
int thIdx = threadIdx.x;
__shared__ int shArr[blockSize*2];
__shared__ int offset;
shArr[thIdx] = thIdx<arraySize ? gArr[thIdx] : 0;
if (thIdx == 0)
offset = blockSize;
__syncthreads();
while (offset < arraySize) { //uniform
shArr[thIdx + blockSize] = thIdx+offset<arraySize ? gArr[thIdx+offset] : 0;
__syncthreads();
if (thIdx == 0)
offset += blockSize;
int sum = shArr[2*thIdx] + shArr[2*thIdx+1];
__syncthreads();
shArr[thIdx] = sum;
}
__syncthreads();
for (int stride = 1; stride<blockSize; stride*=2) { //uniform
int arrIdx = thIdx*stride*2;
if (arrIdx+stride<blockSize)
shArr[arrIdx] += shArr[arrIdx+stride];
__syncthreads();
}
if (thIdx == 0)
*out = shArr[0];
}
...
sumNoncommSingleBlock<<<1, blockSize>>>(dev_a, dev_out);
最初のwhileループでは、スレッドより多くの入力要素がある限り実行されます。各反復では、1回の縮小が実行され、結果はshArr配列の前半に圧縮されます。後半は新しいデータで満たされます。
すべてのデータがgArrからロードされると、2番目のループが実行されます。さて、結果は圧縮されなくなり__syncthreads()追加の__syncthreads() ）。各ステップにおいて、スレッドnは、 2*n番目の能動要素にアクセスし、それを2*n+1番目の要素で加算する。
この単純な例をさらに最適化するには、ワープレベルの削減や共有メモリバンクの競合の除去など、さまざまな方法があります。
交換可能演算子のためのマルチブロック並列化
ブロックが通信に制限されているため、CUDAの並列削減に対するマルチブロックアプローチは、シングルブロックアプローチと比較して追加の課題を提起します。考え方は、各ブロックに入力配列の一部を計算させ、最後にすべての部分的な結果をマージする最後のブロックを持たせることです。そのために、2つのカーネルを起動して、暗黙的にグリッド全体の同期ポイントを作成することができます。
static const int wholeArraySize = 100000000;
static const int blockSize = 1024;
static const int gridSize = 24; //this number is hardware-dependent; usually #SM*2 is a good number.
__global__ void sumCommMultiBlock(const int *gArr, int arraySize, int *gOut) {
int thIdx = threadIdx.x;
int gthIdx = thIdx + blockIdx.x*blockSize;
const int gridSize = blockSize*gridDim.x;
int sum = 0;
for (int i = gthIdx; i < arraySize; i += gridSize)
sum += gArr[i];
__shared__ int shArr[blockSize];
shArr[thIdx] = sum;
__syncthreads();
for (int size = blockSize/2; size>0; size/=2) { //uniform
if (thIdx<size)
shArr[thIdx] += shArr[thIdx+size];
__syncthreads();
}
if (thIdx == 0)
gOut[blockIdx.x] = shArr[0];
}
__host__ int sumArray(int* arr) {
int* dev_arr;
cudaMalloc((void**)&dev_arr, wholeArraySize * sizeof(int));
cudaMemcpy(dev_arr, arr, wholeArraySize * sizeof(int), cudaMemcpyHostToDevice);
int out;
int* dev_out;
cudaMalloc((void**)&dev_out, sizeof(int)*gridSize);
sumCommMultiBlock<<<gridSize, blockSize>>>(dev_arr, wholeArraySize, dev_out);
//dev_out now holds the partial result
sumCommMultiBlock<<<1, blockSize>>>(dev_out, gridSize, dev_out);
//dev_out[0] now holds the final result
cudaDeviceSynchronize();
cudaMemcpy(&out, dev_out, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(dev_arr);
cudaFree(dev_out);
return out;
}
1つは理想的には、GPU上のすべてのマルチプロセッサを完全占有するのに十分なブロックを起動することです。この数を超えると、特に配列内に要素があるほど多くのスレッドを起動すると、生産性が低下します。そうすることで、未処理のコンピューティングパワーはそれ以上上がることはありませんが、非常に効率的な最初のループを使用することはできません。
最後のブロックガードの助けを借りて、単一のカーネルを使って同じ結果を得ることも可能です：
static const int wholeArraySize = 100000000;
static const int blockSize = 1024;
static const int gridSize = 24;
__device__ bool lastBlock(int* counter) {
__threadfence(); //ensure that partial result is visible by all blocks
int last = 0;
if (threadIdx.x == 0)
last = atomicAdd(counter, 1);
return __syncthreads_or(last == gridDim.x-1);
}
__global__ void sumCommMultiBlock(const int *gArr, int arraySize, int *gOut, int* lastBlockCounter) {
int thIdx = threadIdx.x;
int gthIdx = thIdx + blockIdx.x*blockSize;
const int gridSize = blockSize*gridDim.x;
int sum = 0;
for (int i = gthIdx; i < arraySize; i += gridSize)
sum += gArr[i];
__shared__ int shArr[blockSize];
shArr[thIdx] = sum;
__syncthreads();
for (int size = blockSize/2; size>0; size/=2) { //uniform
if (thIdx<size)
shArr[thIdx] += shArr[thIdx+size];
__syncthreads();
}
if (thIdx == 0)
gOut[blockIdx.x] = shArr[0];
if (lastBlock(lastBlockCounter)) {
shArr[thIdx] = thIdx<gridSize ? gOut[thIdx] : 0;
__syncthreads();
for (int size = blockSize/2; size>0; size/=2) { //uniform
if (thIdx<size)
shArr[thIdx] += shArr[thIdx+size];
__syncthreads();
}
if (thIdx == 0)
gOut[0] = shArr[0];
}
}
__host__ int sumArray(int* arr) {
int* dev_arr;
cudaMalloc((void**)&dev_arr, wholeArraySize * sizeof(int));
cudaMemcpy(dev_arr, arr, wholeArraySize * sizeof(int), cudaMemcpyHostToDevice);
int out;
int* dev_out;
cudaMalloc((void**)&dev_out, sizeof(int)*gridSize);
int* dev_lastBlockCounter;
cudaMalloc((void**)&dev_lastBlockCounter, sizeof(int));
cudaMemset(dev_lastBlockCounter, 0, sizeof(int));
sumCommMultiBlock<<<gridSize, blockSize>>>(dev_arr, wholeArraySize, dev_out, dev_lastBlockCounter);
cudaDeviceSynchronize();
cudaMemcpy(&out, dev_out, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(dev_arr);
cudaFree(dev_out);
return out;
}
カーネルは、例えば、ワープレベルの並列縮小を使用することによって、より高速にすることができることに留意されたい。
非可換演算子のためのマルチブロック並列化
並列リダクションに対するマルチブロックアプローチは、シングルブロックアプローチと非常によく似ています。グローバル入力配列はセクションに分割し、それぞれを単一のブロックで縮小する必要があります。各ブロックからの部分的な結果が得られたとき、最終的な結果を得るために、最終的な1つのブロックがそれらを減らす。
-
sumNoncommSingleBlock詳細については、シングルブロック削減の例を参照してください。 -
lastBlockはそれに到達する最後のブロックだけを受け入れます。これを避けたい場合は、カーネルを2つの別々の呼び出しに分割することができます。
static const int wholeArraySize = 100000000;
static const int blockSize = 1024;
static const int gridSize = 24; //this number is hardware-dependent; usually #SM*2 is a good number.
__device__ bool lastBlock(int* counter) {
__threadfence(); //ensure that partial result is visible by all blocks
int last = 0;
if (threadIdx.x == 0)
last = atomicAdd(counter, 1);
return __syncthreads_or(last == gridDim.x-1);
}
__device__ void sumNoncommSingleBlock(const int* gArr, int arraySize, int* out) {
int thIdx = threadIdx.x;
__shared__ int shArr[blockSize*2];
__shared__ int offset;
shArr[thIdx] = thIdx<arraySize ? gArr[thIdx] : 0;
if (thIdx == 0)
offset = blockSize;
__syncthreads();
while (offset < arraySize) { //uniform
shArr[thIdx + blockSize] = thIdx+offset<arraySize ? gArr[thIdx+offset] : 0;
__syncthreads();
if (thIdx == 0)
offset += blockSize;
int sum = shArr[2*thIdx] + shArr[2*thIdx+1];
__syncthreads();
shArr[thIdx] = sum;
}
__syncthreads();
for (int stride = 1; stride<blockSize; stride*=2) { //uniform
int arrIdx = thIdx*stride*2;
if (arrIdx+stride<blockSize)
shArr[arrIdx] += shArr[arrIdx+stride];
__syncthreads();
}
if (thIdx == 0)
*out = shArr[0];
}
__global__ void sumNoncommMultiBlock(const int* gArr, int* out, int* lastBlockCounter) {
int arraySizePerBlock = wholeArraySize/gridSize;
const int* gArrForBlock = gArr+blockIdx.x*arraySizePerBlock;
int arraySize = arraySizePerBlock;
if (blockIdx.x == gridSize-1)
arraySize = wholeArraySize - blockIdx.x*arraySizePerBlock;
sumNoncommSingleBlock(gArrForBlock, arraySize, &out[blockIdx.x]);
if (lastBlock(lastBlockCounter))
sumNoncommSingleBlock(out, gridSize, out);
}
1つは理想的には、GPU上のすべてのマルチプロセッサを完全占有するのに十分なブロックを起動することです。この数を超えると、特に配列内に要素があるほど多くのスレッドを起動すると、生産性が低下します。そうすることで、未処理のコンピューティングパワーはそれ以上上がることはありませんが、非常に効率的な最初のループを使用することはできません。
交換可能演算子のための単一ワープ並列リダクション
場合によっては、より大きなCUDAカーネルの一部として、縮小を非常に小さいスケールで実行する必要があります。たとえば、入力データに正確に32の要素、つまりワープ内のスレッドの数があるとします。そのようなシナリオでは、縮小を実行するために単一のワープを割り当てることができます。ワープが完璧な同期で実行される場合、ブロックレベルの削減と比較して、多くの__syncthreads()命令を削除することができます。
static const int warpSize = 32;
__device__ int sumCommSingleWarp(volatile int* shArr) {
int idx = threadIdx.x % warpSize; //the lane index in the warp
if (idx<16) shArr[idx] += shArr[idx+16];
if (idx<8) shArr[idx] += shArr[idx+8];
if (idx<4) shArr[idx] += shArr[idx+4];
if (idx<2) shArr[idx] += shArr[idx+2];
if (idx==0) shArr[idx] += shArr[idx+1];
return shArr[0];
}
shArrは、好ましくは共有メモリ内の配列である。この値は、ワープ内のすべてのスレッドで同じにする必要があります。場合sumCommSingleWarp複数の経糸によって呼び出され、 shArr経糸（各ワープ内で同じ）との間で異なるべきです。
引数shArrは、配列上の操作が指示されているところで実際に実行されるように、 volatileとしてマークされます。さもなければ、 shArr[idx]反復的な割り当てはレジスタへの割り当てとして最適化され、最終的なアサイメントのみがshArrに対する実際のストアとshArr 。そのようなことが起こると、即時の割り当ては他のスレッドには見えず、結果が不正確になります。非constをconstパラメータとして渡すときと同じように、通常の不揮発性配列をvolatileの引数として渡すことができることに注意してください。
縮小後にshArr[1..31]内容を気にかけなければ、コードをさらに単純化することができます：
static const int warpSize = 32;
__device__ int sumCommSingleWarp(volatile int* shArr) {
int idx = threadIdx.x % warpSize; //the lane index in the warp
if (idx<16) {
shArr[idx] += shArr[idx+16];
shArr[idx] += shArr[idx+8];
shArr[idx] += shArr[idx+4];
shArr[idx] += shArr[idx+2];
shArr[idx] += shArr[idx+1];
}
return shArr[0];
}
このセットアップでは、多くのif条件を削除しました。余分なスレッドはいくつかの不要な追加を実行しますが、生成する内容についてはもう気にしません。ワープはSIMDモードで実行されるので、実際には何もしないことで時間を節約することはできません。一方、これらのif文の本体は非常に小さいため、条件の評価には比較的長い時間がかかります。 shArr[32..47]に0を埋め込むと、最初のif文も削除できます。
ワープレベルの削減は、ブロックレベルの削減を強化するためにも使用できます。
__global__ void sumCommSingleBlockWithWarps(const int *a, int *out) {
int idx = threadIdx.x;
int sum = 0;
for (int i = idx; i < arraySize; i += blockSize)
sum += a[i];
__shared__ int r[blockSize];
r[idx] = sum;
sumCommSingleWarp(&r[idx & ~(warpSize-1)]);
__syncthreads();
if (idx<warpSize) { //first warp only
r[idx] = idx*warpSize<blockSize ? r[idx*warpSize] : 0;
sumCommSingleWarp(r);
if (idx == 0)
*out = r[0];
}
}
引数&r[idx & ~(warpSize-1)]は基本的にr + warpIdx*32です。これにより、 r配列を効果的に32要素のチャンクに分割し、各チャンクは別々のワープに割り当てられます。
非可換オペレータのための単一ワープ並列リダクション
場合によっては、より大きなCUDAカーネルの一部として、縮小を非常に小さいスケールで実行する必要があります。たとえば、入力データに正確に32の要素、つまりワープ内のスレッドの数があるとします。そのようなシナリオでは、縮小を実行するために単一のワープを割り当てることができます。ワープが完璧な同期で実行される場合、ブロックレベルの削減と比較して、多くの__syncthreads()命令を削除することができます。
static const int warpSize = 32;
__device__ int sumNoncommSingleWarp(volatile int* shArr) {
int idx = threadIdx.x % warpSize; //the lane index in the warp
if (idx%2 == 0) shArr[idx] += shArr[idx+1];
if (idx%4 == 0) shArr[idx] += shArr[idx+2];
if (idx%8 == 0) shArr[idx] += shArr[idx+4];
if (idx%16 == 0) shArr[idx] += shArr[idx+8];
if (idx == 0) shArr[idx] += shArr[idx+16];
return shArr[0];
}
shArrは、好ましくは共有メモリ内の配列である。この値は、ワープ内のすべてのスレッドで同じにする必要があります。場合sumCommSingleWarp複数の経糸によって呼び出され、 shArr経糸（各ワープ内で同じ）との間で異なるべきです。
引数shArrは、配列上の操作が指示されているところで実際に実行されるように、 volatileとしてマークされます。さもなければ、 shArr[idx]反復的な割り当てはレジスタへの割り当てとして最適化され、最終的なアサイメントのみがshArrに対する実際のストアとshArr 。そのようなことが起こると、即時の割り当ては他のスレッドには見えず、結果が不正確になります。非constをconstパラメータとして渡すときと同じように、通常の不揮発性配列をvolatileの引数として渡すことができることに注意してください。
shArr[1..31]最終的な内容を気にせず、 shArr[32..47]に0をshArr[32..47]ことができれば、上記のコードを単純化することができます：
static const int warpSize = 32;
__device__ int sumNoncommSingleWarpPadded(volatile int* shArr) {
//shArr[32..47] == 0
int idx = threadIdx.x % warpSize; //the lane index in the warp
shArr[idx] += shArr[idx+1];
shArr[idx] += shArr[idx+2];
shArr[idx] += shArr[idx+4];
shArr[idx] += shArr[idx+8];
shArr[idx] += shArr[idx+16];
return shArr[0];
}
このセットアップでは、命令の半分を構成するif条件をすべて削除しました。余分なスレッドはいくつかの不要な追加を実行し、最終的な結果に最終的に影響を与えないshArrセルに結果を格納します。ワープはSIMDモードで実行されるので、実際には何もしないことで時間を節約することはできません。
レジスタのみを使用したシングルワープ並列リダクション
通常、グローバルまたは共有配列で縮小が実行されます。ただし、CUDAカーネルの一部として、非常に小さなスケールで縮小を実行すると、1回のワープで実行できます。これが起こると、Keppler以上のアーキテクチャ（CC> = 3.0）では、ワープシャッフル機能を使用して共有メモリをまったく使用しないようにすることができます。
たとえば、ワープ内の各スレッドが単一の入力データ値を保持しているとします。すべてのスレッドには合計32個の要素があり、合計する必要があります（または他の連想操作を実行する）
__device__ int sumSingleWarpReg(int value) {
value += __shfl_down(value, 1);
value += __shfl_down(value, 2);
value += __shfl_down(value, 4);
value += __shfl_down(value, 8);
value += __shfl_down(value, 16);
return __shfl(value,0);
}
このバージョンは、可換性と非可換性の両方の演算子で動作します。