本記事ではCUDAカーネル関数を作成し、行列転置など5つの問題を例に、GPUの多数のコアを利用して大規模問題の解決を加速します。
0.1 行列転置
0.1.1 アルゴリズムの流れ
転置する行列をGPUメモリに格納
転置後の行列を格納するためのGPUメモリを確保
行列転置を実行するCUDAカーネル関数を定義。このカーネル関数は、スレッドブロックとスレッドグリッドの概念を使用して行列のすべての要素を処理。各スレッドブロックでは、共有メモリを使用してデータを処理。最後に、結果をグローバルメモリに書き戻す
CUDAカーネル関数を呼び出して行列転置を実行
転置後の行列をGPUメモリからホストメモリにコピー
GPUメモリを解放
0.1.2 コード実装
従来のコードをベースに、共有メモリを利用してグローバルメモリへのアクセスを最適化し、パディング操作によってバンクコンフリクトを回避します。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
const int N = 10000 ;
const int TILE_DIM = 32 ;
const int grid_size_x = ( N + TILE_DIM - 1 ) / TILE_DIM ;
const int grid_size_y = grid_size_x ;
const dim3 block_size ( TILE_DIM , TILE_DIM );
const dim3 grid_size ( grid_size_x , grid_size_y );
__global__ void transpose ( const float * A , float * B , const int N ){
__shared__ float S [ TILE_DIM ][ TILE_DIM + 1 ];
int bx = blockIdx . x * TILE_DIM ;
int by = blockIdx . y * TILE_DIM ;
//順序読み込み
int nx1 = bx + threadIdx . x ;
int ny1 = by + threadIdx . y ;
if ( nx1 < N && ny1 < N )
S [ threadIdx . y ][ threadIdx . x ] = A [ ny1 * N + nx1 ];
__syncthreads ();
//順序書き込み
int nx2 = bx + threadIdx . y ;
int ny2 = by + threadIdx . x ;
if ( nx2 < N && ny2 < N )
B [ nx2 * N + ny2 ] = S [ threadIdx . x ][ threadIdx . y ];
}
0.2 行列乗算
0.2.1 カーネル関数の設計
共有メモリの使用 : 2つの元の行列の一部をタイルとして抽出し、共有メモリに配置することで、グローバルメモリへのアクセスを減らし、プログラムのパフォーマンスを向上させます。
非正方行列の処理 : 非正方行列の乗算では、共有メモリ内のタイルサイズを調整し、境界チェックを追加して範囲外アクセスを防ぐ必要があります。
0.2.2 コード実装
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
__global__ void matrixmul ( const float * A , const float * B , float * C , const int N , const int M , const int K ){ //N*M x M*K
__shared__ float A_S [ TILE_DIM ][ TILE_DIM ];
__shared__ float B_S [ TILE_DIM ][ TILE_DIM ];
int tx = threadIdx . x ;
int ty = threadIdx . y ;
int bx = blockIdx . x ;
int by = blockIdx . y ;
int row = by * TILE_DIM + ty ;
int col = bx * TILE_DIM + tx ;
float value = 0 ;
for ( int ph = 0 ; ph < M / TILE_DIM + 1 ; ph ++ ) {
if ( row < N && ph * TILE_DIM + tx < M )
A_S [ ty ][ tx ] = A [ row * M + ph * TILE_DIM + tx ];
else
A_S [ ty ][ tx ] = 0.0 ;
if ( col < K && ph * TILE_DIM + ty < M )
B_S [ ty ][ tx ] = B [( ph * TILE_DIM + ty ) * K + col ];
else
B_S [ ty ][ tx ] = 0.0 ;
__syncthreads ();
for ( int k = 0 ; k < TILE_DIM ; k ++ )
value += A_S [ ty ][ k ] * B_S [ k ][ tx ];
__syncthreads ();
}
if ( row < N && col < K )
C [ row * K + col ] = value ;
}
0.3 ヒストグラム計算
共有メモリの使用 : 共有メモリに一時配列を確保し、グローバルメモリへのアクセスを減らしてプログラムのパフォーマンスを向上させます。
グリッド間ループ : データ量が大きい場合、1つのスレッドが複数のデータを処理する必要があります。
アトミック操作 : 共有メモリとグローバルメモリの両方で、競合を避けるためにアトミック操作を使用する必要があります。
1
2
3
4
5
6
7
8
9
10
11
12
13
__global__ void histo_kernel ( unsigned char * buffer , long size , unsigned int * histo ){
__shared__ unsigned int temp [ 256 ];
temp [ threadIdx . x ] = 0 ;
__syncthreads ();
int i = threadIdx . x + blockIdx . x * blockDim . x ;
int offset = blockDim . x * gridDim . x ;
while ( i < size ){
atomicAdd ( & temp [ buffer [ i ]], 1 );
i += offset ;
}
__syncthreads ();
atomicAdd ( & ( histo [ threadIdx . x ]), temp [ threadIdx . x ]);
}
0.4 リダクション和
グリッド間ループ : 大規模データ（>100万）の場合、GPUの多数のコアでもすべてのスレッドを処理できないため、1つのスレッドが複数のデータを処理する必要があります。
共有メモリの最適化 : 各スレッドの複数要素の和を共有メモリに格納することで、グローバルメモリへのアクセスを減らすことができます。
インターリーブペアリング : インターリーブペアリング法を使用してリダクションを行うことで、スレッドの分岐を緩和できます。
アトミック操作 : 結果をグローバルメモリに累積する際に競合が発生する可能性があるため、アトミック操作が必要です。
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
__global__ void _sum_gpu ( int * input , int count , int * output )
{
__shared__ int sum_per_block [ BLOCK_SIZE ];
int temp = 0 ;
for ( int idx = threadIdx . x + blockDim . x * blockIdx . x ;
idx < count ; idx += gridDim . x * blockDim . x
)
temp += input [ idx ];
sum_per_block [ threadIdx . x ] = temp ;
__syncthreads ();
//**********共有メモリ集約ステージ***********
for ( int length = BLOCK_SIZE / 2 ; length >= 1 ; length /= 2 )
{
int double_kill = - 1 ;
if ( threadIdx . x < length )
double_kill = sum_per_block [ threadIdx . x ] + sum_per_block [ threadIdx . x + length ];
__syncthreads ();
if ( threadIdx . x < length )
sum_per_block [ threadIdx . x ] = double_kill ;
__syncthreads ();
} //ブロックごとの部分和はsum_per_block[0]に格納
if ( threadIdx . x == 0 ) atomicAdd ( output , sum_per_block [ 0 ]);
}
0.5 TOP K問題
TOP K問題とは、データセット内で最大または最小のK個の要素を見つける問題です。CUDAリダクションを利用して効率的に解決できます。
以下はCUDAリダクションを使用してソートとK個の最大/最小要素の選択を実装する詳細な手順です：
値とインデックスを含むタプル型を定義（ソート後に元のデータを復元するため）
元のデータを走査してタプルに格納
タプルにリダクション操作を適用し、上位K個の最大/最小値のインデックスを取得
元のデータを復元してインデックスでソートし、上位K個の最大/最小値を取得
具体的な実装は以下の通りです：
データをGPUメモリにコピー
1
2
float * d_data ;
cudaMemcpy ( d_data , h_data , size , cudaMemcpyHostToDevice );
データをタプルに格納
1
2
3
4
5
6
7
8
9
typedef struct {
float value ;
int index ;
} Tuple ;
Tuple * d_tuples ;
int threadsPerBlock = 256 ;
int blocksPerGrid = ( n + threadsPerBlock - 1 ) / threadsPerBlock ;
initializeTuples <<< blocksPerGrid , threadsPerBlock >>> ( d_data , d_tuples , n );
タプルにリダクション操作を適用し、上位K個の最大/最小値のインデックスを取得
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
int * d_indices ;
kReduceKernel <<< blocksPerGrid , threadsPerBlock >>> ( d_tuples , d_indices , n , k );
__global__ void kReduceKernel ( Tuple * input , int * output , int n , int k ) {
extern __shared__ Tuple shared [];
int tid = threadIdx . x ;
int i = blockIdx . x * blockDim . x + threadIdx . x ;
shared [ tid ] = ( i < n ) ? input [ i ] : Tuple { 0 , 0 };
__syncthreads ();
for ( int s = blockDim . x / 2 ; s > 0 ; s >>= 1 ) {
if ( tid < s )
shared [ tid ] = ( shared [ tid ]. value > shared [ tid + s ]. value ) ? shared [ tid ] : shared [ tid + s ];
__syncthreads ();
}
if ( tid == 0 )
output [ blockIdx . x ] = shared [ 0 ]. index ;
}
CPU上で元のデータを復元してインデックスでソートし、上位K個の最大/最小値を取得
1
2
3
4
5
6
cudaMemcpy ( h_indices , d_indices , size , cudaMemcpyDeviceToHost );
for ( int i = 0 ; i < k ; ++ i ) {
int index = h_indices [ i ];
h_result [ i ] = h_data [ index ];
}
std :: sort ( h_result , h_result + k );
Tim
支付宝
贝宝
微信