DEIM Forum 2014 A2-1

データ分割と協調的マージに基づく GPU 上の効率的ソートアルゴリズム
小澤 佑介†

天笠 俊之††

北川博之††† 筑波大学大学院システム情報工学研究科 〒 305–8573 茨城県つくば市天王台 1–1–1†† 筑波大学システム情報系 〒 305–8573 茨城県つくば市天王台 1–1–1
E-mail: †kyusuke@kde.cs.tsukuba.ac.jp, ††{amagasa,kitagawa}@cs.tsukuba.ac.jp
あら ま しソートはデータベースをはじめ，コンピュータサイエンスにおける基本的な処理の一つであり，その高速化は非常に重要である．そのため，近年急速に性能が向上している GPU (Graphics Processing Unit) を用いたソートの研究が数多く存在する．GPU を用いたソートアルゴリズムの代表の一つにマージソートベースのものがある．このアルゴリズムには，ソートするデータが巨大な場合に性能が劣化する，という問題がある．本稿では，この欠点を解消した新たなアルゴリズムを提案する．本アルゴリズムは，データを複数のバケットに分割してから，これらのバケットを GPU 上で協調的にソートすることで，高性能を達成する．
キ ーワ ード

GPU, 並列ソート, サンプルソート, マージソート

ソートはコンピュータサイエンスにおける基礎的な処理の一つであり，様々なアルゴリズムやアプリケーションに活用されている [8]．データベースにおいても，ソートは基本的な演算であり，その高速化は非常に重要である [6]．そのため，並列化を含めた多数の研究が存在する [3, 6]．中でも近年では，性能向上が目覚ましい GPU (Graphics Processing Unit) を用いた高速化が注目されている [2, 4, 5, 9, 11, 14, 15, 17]．

Sorting rate (Mkeys/sec)

1. は じ め に
600
500
400
300
200
100
32-bit64-bit0 16 17 18 19 20 21 22 23 24 25 26
2 2 2 2 2 2 2 2 2 2 2

Number of keys

GPU は元来グラフィックス処理用プロセッサであったが，現在は多数のコアを搭載する並列処理用プロセッサとして注目を図 1 Baxter によるマージソート [2] の性能．横軸は入力の要素数，集めている．このように，GPU を一般的な処理の高速化に利

縦軸は性能 (一秒間にどれだけソートできたか) を示している．

用する手法を GPU コンピューティングという [13]．GPU は
低消費電力かつ安価で，非常に高い理論ピーク性能を持つため，

ス回数も増加する．そのため，小さいデータではメモリアクセ

科学計算をはじめ様々な分野で利用が進んでいる．前述の通り，

ス回数も少なく，また GPU の並列度を十分に活かせていない

近年では GPU によるソートの研究も多数存在する．

ため，220 あたりまでは性能が向上していく．一方，データが現在までに提案されている GPU 上のソートアルゴリズムの中で代表的なものとして，基数ソートによるもの [11] と，大きくなってくると，メモリバンド幅がボトルネックとなり，性能が劣化してしまう．

Baxter のマージソート [2] の二つがある．基数ソートは 32-bit

本稿では，以上の欠点を解消した新たな GPU 上のソートアのキーをソートする場合には高速であるが，大きなキーをソー

ルゴリズムを提案する．本アルゴリズムは，はじめにデータをトする場合には大幅に性能が劣化することが知られている．ま複数のバケットに分割することで，後のマージの回数を削減す

た，基数ソートはデータのビット表現を利用するアルゴリズム

る．データの分割はマージと比較して，メモリアクセスが少なであるため，適用できるデータの種類が限られている．そこで

く実現できるため，性能向上が期待できる．また，複数のバケッ

本研究では，より汎用的なアルゴリズムである後者のマージトを作ったうえでソートを行う際に考慮すべき点として，負荷ソートによる実装に着目し，新たな GPU 上のソートアルゴリ

分散がある．本研究ではこの点を考慮して，すべてのバケットズムを提案する．

内を GPU 上で協調的にマージしていく方法を提案する．さら

図 1 は Baxter のマージソートの性能調査の結果を示している．ここでは，一様分布に基づき 32-bit と 64-bit の整数の配に，評価実験によって，提案アルゴリズムが Baxter のマージソートより最大で 39%高速となることを示す．

列をそれぞれ生成し，各配列のソートを行った．この図から，以降の本稿の構成は以下の通りである．2. 節で予備知識につ

要素数が 220 付近で性能がピークを迎え，入力を大きくするにいて説明する．具体的には，GPU コンピューティングと GPU

つれて，性能が低下していくことが分かる．これは，メモリバ上のソートアルゴリズムについて述べる．3. 節で提案アルゴリ

ンド幅がボトルネックとなるためと考えられる．マージソート

ズムについて説明し，4. 節でその評価を行う．最後に 5. 節では計算量が O(n log n) であり，これにしたがってメモリアクセ

まとめと今後の課題について述べる．

2. 予 備 知 識

Algorithm 1: サンプルソート
入力: ソート対象の配列 A

本節では提案アルゴリズムの関連事項について説明する．

1

if |A| > M then

2. 1 節で GPU コンピューティングにおいて利用される CUDA

2S ← 入力 A から αk 個の要素をランダムに選択について，2. 2 節で GPU を用いたソートの既存研究について

3

S をソート

述べる．

4

2. 1 CUDA
GPU コ ン ピ ュ ーティ ン グ とは，元来グラフィックス処理用(s0 , s1 , ..., sk ) ← (−∞, Sα , S2α , ..., S(k−1)α , ∞) ⊲ si は Sの iα 番目の要素．s0 と sk は最小と最大の値を表す番兵．

5
6

for x ∈ A dosi−1 <= x < si となる i を探索

プロセッサであった GPU を，汎用的なプロセッサとして利用

7

し計算の高速化を図る技術を指す．GPU を用いたアプリケー

8

end

ションを作成する際には，NVIDIA 社によって提供されている

9

各バケット Bi に対して再帰的にサンプルソートを適用

CUDA (Compute Unified Device Architecture) （注 1） が主に

10

バケット Bi に x を格納

A ← B1 , ..., Bk を結合したもの用いられる．CUDA が対象とする GPU は以下のような構成に11成っている．GPU は複数の SM (Streaming Multiprocessor)

12を含んでおり，各 SM は 数十から数百の SP (Streaming Pro-

13

end

cessor) を搭載している．CUDA では，一つのカ ーネ ル (GPU14

return A

else
他のアルゴリズムで A をソート上で実行される関数) に対して，数万から数十万の軽量のスレッドを並行に走らせることで，高い並列度を実現している．

2. 2 GPU 上のソ ート ア ルゴリ ズム

CUDA のスレッドは以下のような階層構造を成している．ま既 存 の GPU 上 の ソ ー ト ア ル ゴ リ ズ ム は 大 き く 分 け て ，ず，一つのカーネルの実行にともなって，一つのグリッドが作

基数ソートによるもの [11, 15] と，比較ソートベースのも成される．グリッドは複数のブ ロ ッ ク から構成されており，各の [2, 4, 5, 9, 14, 15, 17] の二種類がある. GPU 上の基数ソート

ブロックは数百のスレッドを持つ．ブロックは 一つの SM にに関する研究はいくつか存在するが，Merrill ら [11] によるも割り当てられ，その SM 上で処理が実行される．ブロック内ののが最速と考えられる．この研究をもとに，更に最適化を施しスレッドは， SM が持つ高速にアクセス可能な共有メ モリ を利た基数ソートの実装が，ライブラリの一部として公開されて用することで，同じブロック内のスレッドとデータのやり取り

いる [10]．基数ソートは，32-bit のキーを対象としたソートのができる．一方，SM は複数のブロックを並行に走らせること

場合，他のアルゴリズムよりも高速である．一方，基数ソートが可能である．このとき，同じ SM に同時に割り当てられていはキーのビット表現を利用しているため，ソート可能なデータるブロックは，SM が持つ資源を共有しながら実行することにが限定されるという欠点がある．また，計算量がソートに用いなる．基本的にはブロック間の同期はできないため，ブロック

るキーのビット長に依存しているため，64-bit 整数キーなど，

間ではデータの依存関係が無いように問題空間を分割する必要

キー長が長い場合には著しく性能が落ちるという問題もある．

がある．ブロック間同期が必要な場合には，複数のカーネルを

マージソートなどの比較ソートベースの研究も数多く存在す

用いることで対処する．

る．ここで比較ソートとは，基数ソートとは異なり，キー同士GPU 向けの最適化で重要な点の一つに，種々のメモリの活の比較のみによってソートするアルゴリズムを指す．GPU 上の用がある．CUDA で主に利用されるメモリには三種類のものが比較ソートの代表的なものとしては，クイックソート [4]，マーある．グ ロ ーバルメ モ リ は GPU 上で最大のメモリで，数 GB

ジソート [2, 5, 15, 17]，バイトニックソート [1, 14]，サンプルの容量を持つ．GPU で処理するデータは基本的にはこのメモ

ソート [9] などがある．これらのアルゴリズムは，キー同士のリ上に置かれる．データの再利用や計算の途中経過を保持する

比較関数さえ用意できれば，ソート対象のデータには制限がなには，高速にアクセス可能なオンチップの共有メモリやレジス

い．ただし，理論上，比較ソートは基数ソートよりも計算量自タが利用可能である．

体は大きくなってしまう [8]．一方，計算量がビット長に依存しグローバルメモリの利用において留意すべき点として，コ ア
レ ス ド ア ク セス がある．通常，グローバルメモリへのアクセスないため，64-bit のキーのソートの場合では基数ソートより高速となる．

は複数のトランザクションにわけられ処理される．しかし，一本研究では，高速な GPU 上の比較ソートアルゴリズムを提つの warp (32 スレッドの集まり) 内のスレッド全てが，アラ案する．以降では，本研究がベースとしている，サンプルソー

インされた 128 バイト領域にアクセスしている場合，これは一ト [3, 9] とマージソートの一実装 [2] について説明する．

つのトランザクションで実行される．これをコアレスドアクセ

2. 2. 1 サンプルソート

スと呼ぶ．コアレスドアクセスが実現されているか否かで大きLeischner ら [9] はサンプルソート [3] の GPU 上の実装を提く性能が変化するため，アルゴリズム設計時に十分に考慮する案した．サンプルソートは，簡単に言うと，クイックソートを必要がある．

拡張したアルゴリズムである．クイックソートが入力を再帰的（注 1）：https://developer.nvidia.com/cuda-toolkitていく．

に二分割していくのに対し，サンプルソートでは k 個に分割しバケットへ分割

タイル内ソート

パディング

B1
2 ブロック/1 マージ

B2

B3

B4

タイル内ソート

4 ブロック/1 マージ

2 ブロック/1 マージ

8 ブロック/1 マージ

4 ブロック/1 マージ

図 2 マージソートの処理の例．

図3

提案アルゴリズムの処理の例．

サンプルソートの擬似コードはアルゴリズム 1 のように書力の配列は複数のタイルに分割され，1 ブロックが 1 タイルをける．入力データがあるしきい値 M より大きければ，以下のソートする．ここでタイルとは，一つのブロックによって独立ように k 個のバケットに分割し再帰的にサンプルソートを適用に処理されるデータ領域を指す．タイルの大きさは，利用するする．

GPU やデータの性質によって適切な値を設定する必要がある．

（ 1 ） 入力 A から k − 1 個の分割点 si を抽出 (2–4 行)

一段階のマージは以下の二つのステップで行っていく．

（ 2 ） 入力の各要素を k 個のバケット Bi に分割 (5–8 行)（ 1 ） Merge Path によって各ブロックが担当する範囲を計算（ 3 ） 各バケット Bi に対しサンプルソートを適用 (9–10 行)（ 2 ） ステップ 1 の計算結果をもとに，各ブロックが並列に入力が十分に小さくなったら他のアルゴリズムによってソート

マージをして出力を完了する．

Leischner らは，以上のアルゴリズムをもとに GPU 上でのサンプルソートを提案した．彼らの実装では，すべてのバケッ

図 2 から分かるようにブロックの数は常に一定で，各ブロックが処理する要素数も一定である．これによって，負荷の偏りが発生しないようにしている．

トがしきい値以下の大きさになるまで，並列に分割していく．

本研究では，上記のマージアルゴリズムをもとに，複数バ

その後，一つのバケットを一つのブロックでソートしていく．

ケット内のマージを一つのカーネルで協調的に行うアルゴリズ

この実装の問題点の一つに，バケットの大きさの違いによる負ムを提案する．多くのバケットに分割することで必要なマージ

荷の偏りがある．さらに，バケットが大き過ぎたり小さ過ぎるの回数を減らすことができ，全体のメモリアクセス回数が少な場合，一つのブロックでソートを行うと，十分に GPU の性能

く済む．そのため，メモリバンド幅がボトルネックとなるようを引き出せないと考えられる．

な場合に高速化が期待できる．

本研究で提案するアルゴリズムでは，サンプルソートのデー
タ分割方法によって，一度だけ分割を行う．そして，すべての

3. 提案ア ルゴリ ズム

バケット内を一つのカーネルで徐々にマージしていく．分割を

本節では，サンプルソートのデータ分割方法と，Merge Path

どのように GPU 上で実現するかについては 3. 1 節で説明する．

を利用した協調的マージに基づく，提案アルゴリズムについて

基本的には Leischner らの研究に基づいているが，いくつかの説明する．提案アルゴリズムの処理の流れは図 3 のようになる．

改良点が存在する．

まず，入力の配列を受け取ったら，この配列を k 個のバケット

2. 2. 2 マージソート

へと分割する．その後，各バケット内をマージソートによって上記のように，GPU 上の比較ソートの研究は多数存在する

ソートしていく．データを k = 2l 個のバケットに分割する場

が，現在最速のものと考えられるのは， Baxter によるマージ合，最大でマージの回数を l 回削減可能である．ただし，バケッ

ソート [2] である．このアルゴリズムは，Merge Path [7, 12]

トサイズにはある程度ずれが生じるため，l 回未満の削減にとを利用した高速なマージに基づいている．Merge Path によっ

どまることもある．以降では，3. 1 節でデータ分割について述て，二つのソートされた配列を並列にマージする際の，各プロ

べ，3. 2 節で協調的マージの説明をする．

セッサへのデータの分割が，高速かつ負荷の偏りが無いように実現できる．データを分割した後は，各プロセッサはそれぞれ

3. 1 データ 分割
ここでは， Leischner らの研究 [9] をもとに改良を加えた，独立に担当範囲のデータをマージし，結果を出力する．

GPU 上のデータ分割アルゴリズムについて説明していく．こBaxter のマージソートは以下の 2 フェイズから成る．

のアルゴリズムは以下の四つのステップからなる．(ここでは k（ 1 ） タイル内ソート：各ブロックがそれぞれ担当する範囲

個のバケットへ分割することを考える．)のデータをソートする（ 1 ） 分割点の選択（ 2 ） マージ：各ブロックがソートしたデータをもとに，協（ 2 ） バケットサイズの計算

調的にマージをしていく（ 3 ） Prefix sum の計算

図 2 はマージソートの処理の流れを示した図である．まず，入（ 4 ） データの再配置

タイル数

バケット数

バケット数

c1,1

c1,2

···

c1,T

s1,1

s1,2

···

s1,T

s1,1

s2,1

···

sk,1

c2,1
..
.

c2,2
..
.

···
..
.

c2,T
..
.

s2,1
..
.

s2,2
..
.

···
..
.

s2,T
..
.

s1,2
..
.

s2,2
..
.

···
..
.

sk,2
..
.

ck,1

ck,2

···

ck,T

sk,1

sk,2

···

sk,T

s1,T

s2,T

···

sk,T

Prefix Sum

タイル数

転置

図 4 Prefix sum の計算前後の処理の流れ．
ここで ci,j はタイル j 内のバケット i に含まれるキーの数を表す．

3. 1. 1 分割点の選択
第一ステップでは，k − 1 個の分割点を入力データから選択

順序で保持される．

3. 1. 3 Prefix Sum の計算する．まず，入力からランダムに αk 個の要素を取得する．こ第三ステップでは，前ステップで計算したバケットサイズ

こで α は oversampling factor であり，大きいほど分割後のバの配列に対して， prefix sum の計算を行う．ここで，pre-ケットサイズが均等になる [3]．実際の分割点 si には，取得しfix sum とは，ある配列 [a1 , a2 , ..., an ] に対する，新たな配列

P0, a1 , a1 + a2 , ..., n−1i=1 ai を指す．本研究では，Yan ら [16]たサンプルをソートしたときの，iα 番目の要素を用いる．
このステップは一つのカーネルとして実現され，1 ブロックののアルゴリズムを実装し利用した．

みで処理を行う．入力から得られたサンプルは，共有メモリに本研究では，このステップで新たに二つの最適化を導入す置きバイトニックソート [1] によりソートする．Leischner らはる．一つめは，prefix sum の行列の転置である．これは，次ス

oversampling factor α を固定の値としていたが，本研究では，テップでコアレスドアクセスを最大限発生させるための最適化

バケットサイズの平滑化のため，共有メモリを可能な限り使用すとなる．次ステップでは，j 番目のタイルを処理するブロック

るように設定する．例えば，共有メモリが 48 KB，k = 128, 入力が，prefix sum 結果の si,j (1 <=i<= k) を読む必要がある．すデータが 32-bit の整数であるとき，α = 48 · 1024 / (4 · k) = 96

なわち，図 4 中央の行列における一列が必要となる．しかし，となる．抽出した分割点は Leischner らが提案しているように，

この行列は行優先で保持されているため，このままでは一つの二分探索木として保持する [9]．これによって，後の，ある要素

warp 中のスレッドがばらばらのアドレスにアクセスしてしまがどのバケットに属するかの計算が高速に行える．

い，コアレスドアクセスが全く発生しない．一方，行列の転置

3. 1. 2 バケットサイズの計算を行うことで，今度は図 4 右の行列における一行を読めばよく第二ステップでは，入力配列を複数のタイルに分割して，入なる．この場合，一つの warp 中のスレッドが連続した領域を力データがどのバケットにいくつ含まれるかのヒストグラム

読むので，コアレスドアクセスの条件を満たす．

を各タイルごとに計算する．この結果に対して，次ステップで二つめの最適化は，バケットのアラインメントである．このprefix sum (後述) を計算することで，データ分割後のバケット最適化もコアレスドアクセスに関するもので，次ステップとの位置が決定できる．

マージに影響する．コアレスドアクセスは，2. 1 節で述べたよ

一つのブロックが一つのタイルを担当し，各ブロックが並列

うに，一つの warp 内のスレッドがア ラ イ ン さ れた 128 バイに計算を行う．タイル内のデータに対するバケットサイズを計

ト 領域にアクセスしている場合に発生する．しかし，各バケッ

算するために，共有メモリ上に，現在のタイル内におけるバトの開始アドレスはデータやパラメータによって決まり，通常

ケットのデータ数を保持しておく．ブロック中のスレッドは，

128 バイト境界にアラインされていない．すなわち，あるバケッ

タイル内の異なる要素を並列に調べる．前ステップで作成したトに対する一つの warp の読み書きは，少なくとも二つ以上の分割点の二分探索木を用いることで，各要素がどのバケットに

トランザクションとして処理される．これを解消するためには，属するかを決定する．そして，この結果をもとに，共有メモリ

各バケットの開始アドレスを調べ，アラインメントが達成され上の対応するバケットサイズを 1 増加させる．この際，同じバるように調整をする必要がある．これは，prefix sum 結果のケットのサイズを同時に変更する可能性があるため，atomic 関

si,1 (1 <=i<= k) を見ることで容易に実現可能である．

数を用いて加算を行う．

3. 1. 4 データの再配置

タイル内の計算が終了したら，そのタイルのヒストグラムを

ここまでに計算した結果をもとにデータの再配置を行う．第グローバルメモリに出力する．このとき，異なるタイルの同じ

二ステップの場合と同様に，入力配列を複数のタイルに分割しバケット ID のバケットサイズが連続するようにメモリ上に配て，各タイルを並列に処理する．ブロックはタイル内のデータ

置する．こうすることで，次ステップの prefix sum の計算で，に対して，どのバケットに属するかの判定をもう一度行う．こデータ分割後における各タイルの各バケットが，どの位置から

れは，再計算した方が，計算結果を保持してそれを読むよりも，はじまるかが分かるようになる．最終的に，出力結果は図 4 左

高速に処理できるためである [9]．判定を行った後，対応するの行列のようになる．この行列はグローバルメモリ上に行優先

バケットサイズを atomic 関数により増加させる．この atomic

t1

t3

t2

t4

t5

t7

t6

···
アドレス

t1

···
112

116

120

124

128

132

136

140

t3

t2

t4

t1

t3

t2

···
アドレス

···
112

116

120

124

128

132

バケット開始アドレス

バケット開始アドレス(a) 最適化前(b) 最適化後

136

140

図 5 バケット位置を元に戻す際の最適化．
この図は 32-bit 要素の配列に，複数のスレッドが並列にアクセスする様子を表している．
ここで，ti はあるブロック中の i 番目のスレッドを意味する．

関数の返り値は，増加させる前の値であり，この値と前ステップで計算した prefix sum の値を調べることで，ある要素をど

表 1 実験に用いた環境

OSの位置に出力すべきかが決定できる．前ステップでアラインメ

Ubuntu 12.04.3

ントのために調整を行ったため，データの再配置後は，図 3 に

LTS 64 bit

ToolkitGPU

CUDA 5.5

Tesla C2050(14 SM, 32 SP@1.15 GHz,
メモリ 3 GB, ECC on)

示されているように，バケット間に隙間がある状態となる．
ここで，グローバルメモリへの出力を単純に行うと，一つの最後のマージでは，データ分割時に行ったバケットのアライ

warp 中のスレッドがばらばらな位置に出力してしまう．すな

ンメントを元に戻す必要がある．そのためには，マージ結果

わち，2. 1 節で述べたコアレスドアクセスの条件を満たさず，を元のバケット開始アドレスに出力すればよい．しかし，単純性能が劣化してしまう．そこで，一度共有メモリ上でデータにこの処理を行うと，3. 1. 4 節で述べたように完全なコアレス

を並び替えてから，グローバルメモリに出力する．これによっ

ドアクセスにはならない．これは，図 5(a) から分かるように，て，可能な限りコアレスドアクセスが発生するようにできる．

1 warp 中のスレッドのアクセスが 128 バイト境界を跨いでしLeischner ら [9] は単純にグローバルメモリに出力したほうがまっているためである．この問題を解決するために， 128 バイ

高速であると述べていたが，本研究の予備実験では，一度共有

ト境界前後の書き込みを別々に扱う．すなわち，128 バイト境

メモリ上で並び替えたほうが高速であったため，こちらを採用

界前の出力は必要なスレッドだけが行い，境界後は 1 番目のス

した．

レッドから順に出力を行うように割り当てる（図 5(b)）．これ

3. 2 協調的マ ージ

によって，境界後の出力は全てコアレスドアクセスの条件を満

前節のアルゴリズムによってデータを k 個のバケットに分たすようにできる．

割した後，バケット内をソートしていく．本研究ではこれを2. 2. 2 節で説明したマージソートのアルゴリズムを基に行う．
マージをはじめる前に，あるブロックがどのバケットを担当

4. 評 価 実 験
本節では，3. 節で提案したアルゴリズムの性能を実験によりするかを決定する必要がある．これは，バケットごとにサイズ評価する．実験に用いた環境は表 1 に示すとおりである．まず，が異なり，静的には割り当てが決定できないためである．まず，

4. 1 節で，提案アルゴリズムと既存アルゴリズムとの比較を行

データ分割中に計算した結果をもとに，バケットサイズを保持

う．4. 2 節以降では，提案アルゴリズムの詳細な分析をしていする配列を作る．この配列から，各バケットが何個のタイルにく．具体的には，4. 2 節でバケット数を変化させた場合の性能分割されるかを計算し，その結果の prefix sum を計算する．この変化を示す．また，バケットサイズとマージの回数について

れによって，各ブロックに関する以下の二つの情報を取得する．

4. 3 節で述べ，4. 4 節で最適化の効果について検証する．

•

どのバケットを処理するか

4. 1 既存ア ルゴリ ズム と の比較•

同じバケットを処理するブロックのうち何番目のブロッ

ここでは，提案アルゴリズム，Baxter のマージソート [2]，ブロックの割り当てが終了したら，マージソートの場合と同様分布による以下の四種類のデータを用いて比較し，のちに他

クかMerrill の基数ソート [10] の三つの性能を比較する．まず，一様に，タイル内ソートを行う．ここでの変更としては，バケッの分布による比較を行う．

トの境界は跨がないようにソートをするという点がある．通常，（ 1 ） 32-bit 整数のキーバケットはタイルによって綺麗に分割はできないので，図 3 に（ 2 ） 64-bit 整数のキー

見られるように，マージソートの場合と比較して，必要なブ（ 3 ） 32-bit 整数のキーと値のペアロック数は多少増加する．その後，バケット内のタイルを対象（ 4 ） 64-bit 整数のキーと値のペアにマージを進めていく．このとき，図 3 からわかるように，対(3) と (4) で，値とはキーに付随するデータであり，ソートは応するタイルが存在しない場合には，何もせずそのまま出力す

キーをもとに行なう．キーと値のソートの応用例としては，関る．マージを繰り返していき，全てのバケット内がソートされ係データベースのリレーションのソートがある．例えば，キーたら終了となる．図 2 と図 3 を比較すると，マージの回数が削としてソートに用いる属性を，値としてタプル ID を利用する

減できていることがわかる．

ことが考えられる．

500
400
300
200
100
proposed

merge

radix

0 19 20 21 22 23 24 25 26 27
2
2
2
2
2
2
2
2
2

400
350
300
250
200
150
100
50

proposed

merge

radix

0 19 20 21 22 23 24 25 26 27
2
2
2
2
2
2
2
2
2

Number of keys

400
350
300
250
200
150
100
50

proposed

Number of keys(a) 32-bit 整数

merge

radix

0 19 20 21 22 23 24 25 26 27
2
2
2
2
2
2
2
2
2

Sorting rate (Mkeys/sec)

600

Sorting rate (Mkeys/sec)

Sorting rate (Mkeys/sec)

Sorting rate (Mkeys/sec)

700

250
200
150
100
50
proposed(b) 64-bit 整数

merge

radix

0 19 20 21 22 23 24 25 26 27
2
2
2
2
2
2
2
2
2

Number of keys

Number of keys(c) 32-bit 整数のペア(d) 64-bit 整数のペア

図 6 一様分布によって生成したデータをソートしたときの性能．
縦軸は性能（一秒間にどれだけソートできたか）を示している．
表2

データサイズを変化させ，ソートの性能を計測した結果を計算時間の内訳とカーネル呼び出し回数（32-bit 整数，キー数 225 ）

図 6 に示す．ここで性能は，同じ分布に基づき 10 種類のデータを作成し，それらをソートしたときの性能の平均とした．分

カーネル

割数 k は，予備実験をもとに，最適な値を設定した．
（詳細は提案アルゴリズムマージソート [2]

時間 (ms)

時間 (ms) 回数

回数

4. 2 節で示す．
）なお，メモリ不足などで正常に動作しなかった

タイル内ソート

16.5

1

13.7

1

箇所は，図では省いている．また，整数のペアを用いた実験結

マージ

22.4

6

56.5

14

果は，整数のみのソートの結果とほぼ同じ傾向となったため，

分割点の選択0.8

1以降では基本的には説明しないこととする．

バケットサイズの計算

4.4

1

Prefix sum

0.3

1

転置

0.3

1

データの再配置

8.5

1合には，提案アルゴリズムが良い性能を示している．例えば

32-bit 整数の場合（図 6(a)），データサイズが 222 以上のとき

合計に，提案アルゴリズムがマージソートより高速となる．特に，倍の性能となった．また，64-bit 整数が対象の場合（図 6(b)），
最大でマージソートの 1.39 倍となり，32-bit 整数の場合より多少有効性が高い結果となった．これは，64-bit 整数の方が必
要メモリバンド幅が大きいため，データ分割によるマージ回数の削減がより有効に働くからと考えられる．一方，データが小さい場合には提案アルゴリズムよりマージソートが高速となった．これは，逆に分割をしすぎて GPU の性能を十分に引き出せない状態となってしまうためと考えられる．このような場合

Sorting rate (Mkeys/sec)

キー数が 225 のとき，提案アルゴリズムはマージソートの 1.32

53.2

700
600
500
400
300
200
100
proposed

merge

radix

0 19 20 21 22 23 24 25 26 27
2
2
2
2
2
2
2
2
2

Sorting rate (Mkeys/sec)

マージソートと比較して，データサイズが十分に大きい場

70.2

900
800
700
600
500
400
300
200
100

Number of keys(a) 正規分布（平均 230 ，標準偏差 226 ）proposed

merge

radix

0 19 20 21 22 23 24 25 26 27
2
2
2
2
2
2
2
2
2

Number of keys(b) 指数分布（λ = 10−3 ）図 7 一様分布以外の分布を利用したときのソート性能比較．

には，マージソートをそのまま利用したほうが良い．
基 数 ソ ー ト と 比 較 す る と ，32-bit 整 数 の ソ ー ト の 場 合
2225ているが，提案アルゴリズムとマージソートのマージ回数の差の間では，提案アルゴリが 8 と，7 より多くなっている．これは，それぞれのタイルサ

ズムが基数ソートと同等かそれ以上の性能を示している．一方，イズの設定が異なるためである．また，提案アルゴリズムのタ

キー数が増加するにつれて，基数ソートが優位となっていく．

イル内ソートが少し遅くなっているが，これは 3. 2 節で述べた

基数ソートは計算量がデータサイズに比例するため，データがとおり，必要なタイル数が多少増加してしまうためと考えられ大きくなっても性能は劣化せずに済む．しかし，マージソート

る．データ分割は 14.3 ms で実現できており，全体の処理時間の計算量は O(n log n) であるため，データサイズが 2 倍になる

は提案アルゴリズムが 17 ms 速い結果となった．

（図 6(a)），キー数が 2から 2と処理時間は 2 倍超となるため，性能は徐々に低下していくこ

最後に，一様分布以外の分布を用いて生成した 32-bit 整数とになる．64-bit 整数の場合（図 6(b)）では，2. 2 節で触れた

キーの配列を用いて性能比較を行う．図 7 に，正規分布と指数

とおり基数ソートの性能は大幅に劣化しており，提案アルゴリ

分布を利用した結果を示す．正規分布の場合（図 7(a)）は一様

ズムが高速となる．特に，キー数が 222 のときに，基数ソート

分布のときとほぼ同じ結果となった．一方，指数分布を用いたと比較して提案アルゴリズムが 2.33 倍の性能を達成している．

とき（図 7(b)）は，提案アルゴリズムとマージソートには大き次に，提案アルゴリズムとマージソートのより詳細な比較をな変化がなかったが，基数ソートの性能が大幅に向上した．こ行う．表 2 に，それぞれの計算時間の内訳と主なカーネルの呼れは，指数分布では多くのキーが小さな値になるためと考えら

び出し回数を示す．この表から，提案アルゴリズムのマージ回れる．この場合，上位ビットのほとんどが 0 となり，基数ソー

数は 6 回であり，マージソートの 14 回と比較して 8 回削減でトではソートが早い段階で完了する．そのため，一様分布の場

きていることがわかる．ここでは分割数を k = 128 = 27 とし合と比較して，性能が向上した．

400
300
200
100
k = 26

0

2

19

2

20

2

k = 27

21

2

22

2

23

2

24

k = 28

2

25

2

26

2

350
300
250
200
150
100
50
k = 26

0

27

2

19

2

Number of keys

20

2

21

k = 27

2

22

2

23

2

24

k = 28

2

25

2

26

2

27

400
350
300
250
200
150
100
50
k = 26

0

2

19

2

20

2

Number of keys(a) 32-bit 整数

Sorting rate (Mkeys/sec)

500

400

Sorting rate (Mkeys/sec)

600

Sorting rate (Mkeys/sec)

Sorting rate (Mkeys/sec)

700

21

k = 27

2

22

2

23

2

24

k = 28

2

25

2

26

2

27

250
200
150
100
50
k = 260

Number of keys(b) 64-bit 整数

k = 27k = 28219 220 221 222 223 224 225 226 227

Number of keys(c) 32-bit 整数のペア(d) 64-bit 整数のペア

図 8 異なるバケット数を用いた場合の性能比較．
表 3 最大バケットサイズ（32-bit 整数，バケット数 k = 27 , タイルサイズ 5888）
データサイズ

219

220

221

222

223

224

225

226

227

最大バケットサイズ

5,643 11,252 23,752 49,854 88,789 183,971 340,139 736,388 1,415,086

しきい値

5,888 11,776

23,552

47,104 94,208 188,416 376,832 753,664 1,507,328

4. 2 バケ ッ ト 数の違いによ る 性能変化

1.30

本節では，バケット数 k の違いによる性能の変化について
7

8

ケット数 k を 2 , 2 , 2 としたときの性能を図 8 に示す．この図から，32-bit 整数と 32-bit 整数のペアでは k = 27 が最適
78な値となることがわかる．64-bit 整数では，k = 2 と k = 2の場合が良い性能を示しており，データサイズが 222 –225 では

k = 28 が，それより大きい場合には k = 27 が優位となる．本研究で行った実験では，大きなデータに対して有効な k = 27を採用した．一方，64-bit 整数のペアでは k = 28 が最適な値

1.25
1.20

Speedup6

1.20

Speedup

調査する．前節と同様に一様分布を用いてデータを生成し，バ

1.30Transposition
Alignment
Transposition+
Alignment

1.25

1.15
1.10
1.05

1.15
1.10
1.05

1.00
0.95

Transposition
Alignment
Transposition+
Alignment

1.00
2

19

2

20

2

21

2

22

2

23

2

24

2

25

2

26

2

27

Number of keys(a) 32-bit 整数

0.95

219 220 221 222 223 224 225 226 227

Number of keys(b) 64-bit 整数

図 9 転置（Transposition）とアラインメント（Alignment）を適用したときの性能向上率．

となった．これは，要求されるメモリバンド幅が四種類のデー
タの中で最も大きく，より多くマージ回数を削減できる方が良の違いによるものと考えられる．アラインメントは，データ分

いためと考えられる．

割以降全体に影響を与えるが，転置はデータの再配置にしか影

4. 3 バケ ッ ト サイ ズと マ ージ 回数

響しない．データ分割以降は全体の処理時間の半分以上を占め

ここでは，データ分割によって作られるバケットの大きさとる部分であるため，アラインメントの効果が大きい結果となっマージ回数について調べる．特に，ある一様分布によって生成た．ただし，図 9(a) におけるキー数 227 のように，転置のみをされた 10 種類のデータを分割したときの，最大のバケットサイ

適用した場合より，同時に最適化を施した方が転置がより効果ズに注目する．これは，最大のバケットサイズによって，マー的となることがある．これは以下のように考えられる．最適化

ジカーネルの呼び出し回数が決まるためである．表 3 に結果ををしていない場合，再配置のカーネルは，入力と出力の両方で示す．表中のしきい値は，バケット数が k = 2 の場合に，最コアレスドアクセスの条件を満たさない状態である．そのため，

大バケットサイズがしきい値以下ならば l 回マージが削減可能

転置によって入力はコアレスドアクセスになるが，出力は変わとなるような値である．表 3 から，データサイズが 221 と 222

らないため，こちらがボトルネックとなる．そこで，アライン

以外の場合には，最大バケットサイズはしきい値以下であり，lメントを行うことで，出力もコアレスドアクセスになり，問題

回マージが削減できていることがわかる．

が解消される．

l

4. 4 最適化の効果

また，キー数が増加するにしたがって，性能向上率も上昇し最後に，3. 1. 3 節で導入した二つの最適化の効果を検証する．

ていくことが見てとれる．さらに，32-bit 整数より 64-bit 整

具体的には，prefix sum の計算後の行列の転置と，バケット開数の場合に最適化が有効であることもわかる．これは，データ

始位置のアラインメントである．後者の最適化には，3. 2 節でが大きくなるにつれて，要求メモリバンド幅が増加していくた

述べた，最後のマージにおける最適化も含めている．図 9 に，めである．二つの最適化はどちらも，コアレスドアクセスに関

最適化を適用した際の性能向上率を示す．この図では，二つのするものであり，これらが有効に機能していることがわかる．

最適化を適用していない場合の性能を 1 としている．

一方，データが小さい場合には，逆に最適化によって遅くなっ

図から，二つの最適化のうち，アラインメントがより効果的てしまうこともある．これは，最適化の導入によって発生するであることがわかる．これは二つの最適化の影響範囲の大きさ

オーバーヘッドによるものと考えられる．例えば図 9(a) における 219 では，マージの回数が 1 回のみであり，データも小さ
いため，最適化によって性能が劣化する結果となった．

5. ま

と

め

本研究では，GPU 上の高効率なソートアルゴリズムを提案した．提案アルゴリズムはサンプルに基づくデータ分割と協調的マージによってソートを行う．はじめにデータ分割を行うことで，マージソートで行われるマージの回数を減らし，必要メ
モリバンド幅の削減を図った．さらに，コアレスドアクセスが最大限発生するように， prefix sum の行列の転置とバケットのアラインメントの二つの最適化を導入した．その後，一つのカーネルで全てのバケット内のマージを行う，協調的マージによってバケット内のソートを完了する．実験により，データサ
イズが十分大きい場合に，マージソートより高速化できている
ことを示した．また，一様分布に基づく 32-bit 整数の配列を用いた場合には，基数ソートと同等の性能となることもわかった．
さらに，64-bit 整数の場合には，基数ソートの 2 倍以上の性能となった．
今後の課題としては，データ分割方法の改良がある．今回はデータ分割の回数を 1 回に限定していたが，これを 2 回以上行うことで必要なメモリアクセスをさらに削減できると考えられる．また，現在の分割方法では atomic 関数を用いているため，ソート結果の安定性は保証されていない．この問題を改善することも検討する予定である．
謝辞

本研究の一部は，平成 25 年度共同研究（富士通研究所 CPE25149），JSPS 科研費（24240015），「エクサスケール
計算技術開拓による先端学際計算科学教育研究拠点の充実」事
業によるものである．

文

献[1] K. E. Batcher, “Sorting networks and their applications,”
Proc. AFIPS Spring Joint Computer Conference, pp. 307–
314, 1968.
[2] S. Baxter, “Modern GPU,” http://nvlabs.github.io/
moderngpu/, (last accessed Feb. 2014).
[3] G. E. Blelloch, C. E. Leiserson, B. M. Maggs, C. G. Plaxton, S. J. Smith, and M. Zagha, “A Comparison of Sorting
Algorithms for the Connection Machine CM-2,” Proc. 3rd
Annual ACM Symp. Parallel Algorithms and Architectures(SPAA), pp. 3–16, 1991.
[4] D. Cederman and P. Tsigas, “GPU-Quicksort: A Practical Quicksort Algorithm for Graphics Processors,” J. Exp.
Algorithmics, vol. 14, no. 4, pp. 1–24, Jul. 2009.
[5] A. Davidson, D. Tarjan, M. Garland, and J. D. Owens, “Efficient Parallel Merge Sort for Fixed and Variable LengthKeys,” Proc. Innovative Parallel Computing (InPar), pp.
131–139, May 2012.
[6] G. Graefe, “Implementing Sorting in Database Systems,”
ACM Comput. Surv., vol. 38, no. 3, article 10, Sep. 2006.
[7] O. Green, R. McColl, and D. A. Bader, “GPU Merge Path:
A GPU Merging Algorithm,” Proc. 26th ACM Int’l Conf.
Supercomputing (ICS), pp. 331–340, 2012.
[8] D. E. Knuth, “The Art of Computer Programming, Volume 3: Sorting and Searching,” Second Edition, AddisonWesley, 1998.

[9] N. Leischner, V. Osipov, and P. Sanders, “GPU Sample
Sort,” Proc. 24th IEEE Int’l Parallel & Distributed Processing Symp. (IPDPS), pp. 1–10, 2010.
[10] D. Merrill,“CUB,”
http://nvlabs.github.io/cub/,ver. 1.1.1, 2013.
[11] D. Merrill and A. Grimshaw. “High Performance and Scalable Radix Sorting: A case study of implementing dynamic
parallelism for GPU computing,” Parallel Process. Lett.,
vol. 21, no. 2, pp. 245–272, Jun. 2011.
[12] S. Odeh, O. Green, Z. Mwassi, O. Shmueli, and Y. Birk,“Merge Path - Parallel Merging Made Simple,” Proc. 26th
IEEE Int’l Parallel & Distributed Processing Symp. Workshops & PhD Forum (IPDPSW), pp. 1611–1618, 2012.
[13] J. D. Owens, M. Houston, D. Luebke, S. Green, J. E. Stone,
and J. C. Phillips, “GPU Computing,” Proc. IEEE, vol. 96,
no. 5, pp. 879–899, May 2008.
[14] H. Peters, O. Schulz-Hildebrandt, and N. Luttenberger, “Anovel sorting algorithm for many-core architectures basedon adaptive bitonic sort,” Proc. 26th IEEE Int’l Parallel &Distributed Processing Symp. (IPDPS), pp. 227–237, 2012.
[15] N. Satish, C. Kim, J. Chhugani, A. D. Nguyen, V. W. Lee,D. Kim, and P. Dubey, “Fast Sort on CPUs and GPUs:
A Case for Bandwidth Oblivious SIMD Sort,” Proc. ACM
SIGMOD Int’l Conf. Management of Data (SIGMOD), pp.
351–362, 2010.
[16] S. Yan, G. Long, and Y. Zhang, “StreamScan: Fast Scan
Algorithm for GPUs without Global Barrier Synchronization,” Proc. 18th ACM SIGPLAN Symp. Principles andPractice of Parallel Programming (PPoPP), pp. 229–238,
2013.
[17] X. Ye, D. Fan, W. Lin, N. Yuan, and P. Ienne, “High Performance Comparison-Based Sorting Algorithm on Many-Core
GPUs,” Proc. 24th IEEE Int’l Parallel & Distributed Processing Symp. (IPDPS), pp. 1–10, 2010.