AWSのP2.*インスタンスで PG-Strom を試す

従前、AWSの提供するGPUインスタンス g2.* に搭載されているGPUはGRID K520というちょっと古いモデルで、PG-Stromは非対応だった。
理由は、一年ほど前にComputing Capability 3.5以降で対応のDynamic Parallelism機能を使うように全面的に作り直したからで、詳細は以下のエントリを参照。

kaigai.hatenablog.com

その後、昨年の10月にAWSは新世代*1のGPUインスタンスを新たにリリースした。

japan.zdnet.com

これでPG-Stromの動作要件を満たすようになった上に、特にメモリ搭載量で相応の強化が行われたため、例えばPGconf.ASIAで発表を行った創薬領域の類似度サーチのような、I/Oが支配的でないようなワークロードであれば相応の効果が見込める、ハズである。

発表から少し間が空いてしまったが、p2.*インスタンスで動作検証を行い、当該環境をAMIイメージとして公開してみた。

環境構築

基本的にはAMIを探してポチポチするだけである。

最初に、リージョンを p2.* インスタンス提供済みのリージョンに切り替える。
現状、東海岸の「バージニア北部」、西海岸の「オレゴン」、欧州の「アイルランド」での提供が開始されている模様。東京はまだである*2。
とりあえず、今回は「オレゴン」リージョンを使用。

また、AWSアカウントを作成直後は p2.* インスタンスの起動制限数が 0 に設定されているため、p2.xlargeインスタンスを起動できるよう、EC2ダッシュボードの『制限』から『制限緩和のリクエスト』を行う必要がある。豆知識。

次に、コミュニティAMIから『PGStrom』で検索すれば、夜なべして作ったAMIイメージがヒットする。
これは PG-Strom v1.0 + PostgreSQL v9.5.5 + CUDA 7.5環境を CentOS7(x86_64) 上に構築したもので、テスト用のデータも既にセットアップ済みのものである。

次に、起動すべきインスタンスを選択する。
p2.*インスタンスは以下の3種類提供されており、今回は『p2.xlarge』タイプを使用した。
なお、Amazon様も『高性能データベース』用途ときちんと書いていてくれておりますｗ

一応、PG-StromはマルチGPUにも対応しているが、CPU並列を使えないPostgreSQL v9.5系列ではGPUへのデータ供給の方が先にボトルネックになってしまうので、あまり意味はない。
この辺は、現在開発中の PostgreSQL v9.6 ベースの実装でうまくハンドリングできるようになる。

インスタンスの構成を確認して起動。ストレージはroot区画に24GBをアタッチするだけのシンプルなものがデフォルト構成。

しばらく待っているとインスタンスが立ち上がってくる。

『ec2-user』でログインする事ができる。

動作確認

早速動作確認を行う。ec2-userはlocalhost経由でPostgreSQLサーバに接続できるよう設定されている。

[ec2-user@ip-172-31-24-196 ~]$ psql postgres
psql (9.5.5)
Type "help" for help.

postgres=#

Tesla K80 GPUが搭載されている事が分かる。
物理的には一枚のカード上にGPU 2個を搭載するデバイスなので、一枚のカードを他の誰かとシェアしている事になる。

postgres=# SELECT * FROM pgstrom_device_info();
 id |                       property                        |    value
----+-------------------------------------------------------+--------------
  0 | Device name                                           | Tesla K80
  0 | Total global memory size                              | 11519 MBytes
  0 | max threads per block                                 | 1024
  0 | Maximum block dimension X                             | 1024
  0 | Maximum block dimension Y                             | 1024
  0 | Maximum block dimension Z                             | 64
  0 | Maximum grid dimension X                              | 2147483647
  0 | Maximum grid dimension Y                              | 65535
  0 | Maximum grid dimension Z                              | 65535
  0 | Maximum shared memory available per block             | 48 KBytes
  0 | Memory available on device for __constant__           | 64 KBytes
  0 | Warp size in threads                                  | 32
  0 | Maximum pitch in bytes allowed by memory copies       | 2147483647
  0 | Maximum number of 32bit registers available per block | 65536
  0 | Typical clock frequency in kilohertz                  | 823 MHz
  0 | Alignment requirement for textures                    | 512
  0 | Number of multiprocessors on device                   | 13
  0 | Has kernel execution timeout                          | False
  0 | Integrated with host memory                           | False
  0 | Host memory can be mapped to CUDA address space       | True
  0 | Compute mode                                          | Default
  0 | Alignment requirement for surfaces                    | 512
  0 | Multiple concurrent kernel support                    | True
  0 | Device has ECC support enabled                        | True
  0 | PCI bus ID of the device                              | 0
  0 | PCI device ID of the device                           | 30
  0 | Device is using TCC driver model                      | False
  0 | Peak memory clock frequency                           | 2505 MHz
  0 | Global memory bus width                               | 384 bits
  0 | Size of L2 cache in bytes                             | 1536 KBytes
  0 | Maximum threads per multiprocessor                    | 2048
  0 | Number of asynchronous engines                        | 2
  0 | Device shares unified address space                   | True
  0 | PCI domain ID of the device                           | 0
  0 | Major compute capability version number               | 3
  0 | Minor compute capability version number               | 7
  0 | Device supports stream priorities                     | True
  0 | Device supports caching globals in L1                 | True
  0 | Device supports caching locals in L1                  | True
  0 | Maximum shared memory per multiprocessor              | 112 KBytes
  0 | Maximum number of 32bit registers per multiprocessor  | 131072
  0 | Device can allocate managed memory on this system     | True
  0 | Device is on a multi-GPU board                        | False
  0 | Unique id if device is on a multi-GPU board           | 0
(44 rows)

早速、シンプルな集計のクエリを叩いてみるが、その前にメインのテーブルに対してpg_prewarmを実行しておくが吉。
gp2ストレージタイプの場合、バッファに載っていないデータをロードするのに割と時間がかかるようなので、I/Oネックになってしまうと、CPUとかGPUとかいう次元の遥か手前で引っかかってしまう。

postgres=# SELECT pg_prewarm('t0'::regclass);
 pg_prewarm
------------
     833334
(1 row)

実行計画を確認する。JOIN+GROUP BYをGPUで実行するようにプランが選択された事がわかる。

postgres=# EXPLAIN SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
                                         QUERY PLAN
---------------------------------------------------------------------------------------------
 HashAggregate  (cost=3578311.16..3578311.48 rows=26 width=12)
   Group Key: t0.cat
   ->  Custom Scan (GpuPreAgg)  (cost=14139.24..2873626.84 rows=234 width=48)
         Reduction: Local + Global
         GPU Projection: cat, ax
         ->  Custom Scan (GpuJoin) on t0  (cost=10139.24..2838812.04 rows=98599882 width=12)
               GPU Projection: t0.cat, t1.ax
               Depth 1: GpuHashJoin, HashKeys: (t0.aid)
                        JoinQuals: (t0.aid = t1.aid)
                        Nrows (in/out: 98.60%), KDS-Hash (size: 13.47MB, nbatches: 1)
               ->  Seq Scan on t1  (cost=0.00..1935.00 rows=100000 width=12)
(11 rows)

1億行のテーブルのJOINとGROUP BYで11.3秒要している。
セッションに接続後の初回実行だと、これに加えて、GPUコンテキストの初期化やGPUコードのJITコンパイルの時間も余分にかかるが、オンプレに比べるとAWSの環境ではこの辺が多少もっさりする印象がある。

postgres=# SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
 cat |  count  |       avg
-----+---------+------------------
 nnn | 3845736 | 49.9122204644341
 ccc | 3842975 | 49.9253161146813
 ddd | 3841209 | 49.9346989307005
 aaa | 3848221 | 49.9265219996308
 kkk | 3843481 | 49.9232348058601
 fff | 3845484 | 49.9434949581969
 iii | 3846743 | 49.9227199719054
 jjj | 3846076 | 49.9292863471083
 qqq | 3845646 |  49.945213365697
 hhh | 3842519 | 49.9266693322143
 ttt | 3846725 | 49.9232478784252
 ooo | 3847927 | 49.9346102129999
 zzz | 3847116 | 49.9320751724648
 lll | 3848447 |  49.927865906661
 www | 3846691 | 49.9537192102014
 bbb | 3845315 | 49.9281298849625
 ppp | 3850842 | 49.9149069792416
 eee | 3846285 |  49.931458202446
 xxx | 3845570 | 49.9577920754281
 ggg | 3845044 | 49.9409383291351
 rrr | 3847816 | 49.9341189910578
 uuu | 3845813 | 49.9295202543591
 vvv | 3849157 | 49.9253944163053
 yyy | 3843414 | 49.9364087656463
 mmm | 3848758 | 49.9033622681507
 sss | 3846990 | 49.9213589517191
(26 rows)

Time: 11326.834 ms

一方、PG-Stromを無効にした場合は以下のように、同じクエリの実行に51.0秒を要している。ﾄｩｯﾄｩﾙｰ。

postgres=# EXPLAIN SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
                                 QUERY PLAN
-----------------------------------------------------------------------------
 HashAggregate  (cost=3937016.94..3937017.26 rows=26 width=12)
   Group Key: t0.cat
   ->  Hash Join  (cost=3185.00..3197517.82 rows=98599882 width=12)
         Hash Cond: (t0.aid = t1.aid)
         ->  Seq Scan on t0  (cost=0.00..1833334.00 rows=100000000 width=8)
         ->  Hash  (cost=1935.00..1935.00 rows=100000 width=12)
               ->  Seq Scan on t1  (cost=0.00..1935.00 rows=100000 width=12)
(7 rows)

postgres=# SELECT cat,count(*),avg(ax) FROM t0 NATURAL JOIN t1 GROUP BY cat;
 cat |  count  |       avg
-----+---------+------------------
 nnn | 3845736 | 49.9122204644368
 ccc | 3842975 |  49.925316114681
 ddd | 3841209 |  49.934698930695
 aaa | 3848221 | 49.9265219996321
 kkk | 3843481 | 49.9232348058579
 fff | 3845484 | 49.9434949581975
 iii | 3846743 |    49.9227199719
 jjj | 3846076 | 49.9292863471066
 qqq | 3845646 | 49.9452133657041
 ttt | 3846725 | 49.9232478784215
 hhh | 3842519 |  49.926669332217
 zzz | 3847116 |  49.932075172461
 ooo | 3847927 | 49.9346102129966
 lll | 3848447 | 49.9278659066662
 www | 3846691 | 49.9537192102029
 bbb | 3845315 | 49.9281298849615
 ppp | 3850842 | 49.9149069792369
 eee | 3846285 |  49.931458202444
 xxx | 3845570 | 49.9577920754269
 ggg | 3845044 | 49.9409383291352
 uuu | 3845813 | 49.9295202543602
 rrr | 3847816 | 49.9341189910555
 vvv | 3849157 | 49.9253944163048
 yyy | 3843414 | 49.9364087656454
 mmm | 3848758 | 49.9033622681445
 sss | 3846990 | 49.9213589517182
(26 rows)

Time: 51031.084 ms

p2.* インスタンスの登場でPG-Stromを使用するための環境の準備はずいぶんお手軽になった。デプロイ即利用できるようAMIイメージも用意してみたので、まだPG-Stromを触った事がないという方は、ぜひ試してみて頂ければと思う。

*1:と言っても、Tesla K80なんですが…。

*2:予定なしとか悲しい事は言わないでほしいけど

2017年の開発ロードマップについて考える

あけましておめでとうございました。（やや出遅れ感）

新年という事で、この一年、どういった技術開発に取り組んでいきたいかをざーっと書き出してみる事にする。
これらのうち、いくつかはPostgreSQL本体機能の強化を伴うものであったりするので、ある程度計画的にモノゴトを進めないといけないワケで…。

PG-Strom v2.0

先ず最優先で取り組むのが、PostgreSQL v9.6への対応。
CPUパラレル実行と、新しいオプティマイザへの対応でかなり大きなアーキテクチャ上の変更を伴ったものの、全体としてはよりシンプルな設計に落とし込む事ができている。

ちなみに、現状だとこの程度までは動くようになっている。
集約演算がGroupAggregateとGpuPreAggの二段階に分解されており、GpuPreAggはGatherの配下で並列に動作している事に注目。

postgres=# EXPLAIN (ANALYZE, VERBOSE)
           SELECT cat, count(*), avg(aid), max(bid)
             FROM t0
            WHERE aid < 50000 and cid > 50000
            GROUP BY cat;

                                           QUERY PLAN
------------------------------------------------------------------------------------------------
 GroupAggregate  (cost=91050.50..91070.84 rows=26 width=48)
                 (actual time=1754.432..1755.056 rows=26 loops=1)
   Output: cat, pgstrom.sum((pgstrom.nrows())),
                pgstrom.favg((pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                             (pgstrom.psum((aid)::bigint))))),
                max((pgstrom.pmax(bid)))
   Group Key: t0.cat
   ->  Sort  (cost=91050.50..91052.32 rows=728 width=48)
             (actual time=1754.381..1754.524 rows=910 loops=1)
         Output: cat, (pgstrom.nrows()), (pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                                                       (pgstrom.psum((aid)::bigint)))),
                 (pgstrom.pmax(bid))
         Sort Key: t0.cat
         Sort Method: quicksort  Memory: 160kB
         ->  Gather  (cost=90938.54..91015.89 rows=728 width=48)
                     (actual time=1749.313..1753.732 rows=910 loops=1)
               Output: cat, (pgstrom.nrows()), (pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                                                             (pgstrom.psum((aid)::bigint)))),
                       (pgstrom.pmax(bid))
               Workers Planned: 4
               Workers Launched: 4
               ->  Parallel Custom Scan (GpuPreAgg) on public.t0  (cost=89938.54..89943.09 rows=182 width=48)
                                                            (actual time=1670.139..1673.636 rows=182 loops=5)
                     Output: cat, (pgstrom.nrows()),
                             pgstrom.pavg((pgstrom.nrows((aid IS NOT NULL))),
                                          (pgstrom.psum((aid)::bigint))),
                             (pgstrom.pmax(bid))
                     Reduction: Local
                     GPU Projection: t0.cat, pgstrom.nrows(), pgstrom.nrows((t0.aid IS NOT NULL)),
                                     pgstrom.psum((t0.aid)::bigint), pgstrom.pmax(t0.bid), t0.aid, t0.cid
                     Outer Scan: public.t0  (cost=4000.00..87793.21 rows=2496387 width=12)
                                            (actualtime=14.663..274.534 rows=500187 loops=5)
                     Outer Scan Filter: ((t0.aid < 50000) AND (t0.cid > 50000))
                     Rows Removed by Outer Scan Filter: 1499813
                     Extra: slot-format
                     Worker 0: actual time=1724.195..1728.058 rows=182 loops=1
                     Worker 1: actual time=1570.952..1573.837 rows=182 loops=1
                     Worker 2: actual time=1738.205..1742.053 rows=182 loops=1
                     Worker 3: actual time=1569.055..1571.961 rows=182 loops=1
 Planning time: 0.907 ms
 Execution time: 1759.557 ms
(25 rows)

また、PG-Strom v2.0では、PostgreSQL v9.6へのキャッチアップだけではなく、いくつか目玉となる機能を準備中である。
一つは、これまで何度か紹介している SSD-to-GPU P2P DMA 機構。そしてもう一つは、BRINインデックスへの対応である。

SSD-to-GPU P2P DMA

SSD-to-GPU P2P DMA (NVMe-Strom) は、NVIDIA社製GPUのGPUDirect RDMA機構を利用したもので、PostgreSQLのデータブロックが格納されているNVMe-SSDのデータブロックからGPUへとダイレクトにデータ転送を行う。ファイルシステムを介する事によるオーバーヘッドや、RAMへの無駄なコピーが発生しないため、スループットを稼げるという特長がある。
現状では、GpuScanワークロード下においてNVMe-SSD 1個から成る区画からのデータ転送に対応しており、シングルプロセス性能で1.4GB/sのスキャン性能を出している。
PostgreSQL v9.6対応の過程で、GpuJoinやGpuPreAggの直下にテーブルスキャンが入る場合、これらのロジックはGpuScanがテーブルをスキャンするための関数を直接呼ぶように改良されているので、特別な事は何もしなくても『ストレージからデータを読んだ時点で既にJOIN/GROUP BYが完了している』という状態を作り出す事はできるはず。

PG-Strom v2.0に向けた課題はSoft-RAID0/1への対応。Linuxの場合、基本的には128KB単位で順番にストライピングがかかっているだけなので、技術的にはそう難しい話ではないと考えている。
DC用途向けに、PCIe x8スロット接続で5～6GB/s程度のSeqRead性能を持つNVMe-SSD製品が各社から出てきているので、計算上は、SSD二枚から全力でGPUにデータを流し込む事ができれば、GPUの持つPCIe x16スロットの帯域を飽和させられる事になる。

BRINインデックス対応

BRINインデックス自体はPostgreSQL v9.5から搭載されている機能で、特に時系列データのように

1. ある一定範囲の値を持つデータが物理的な近傍に集まっている
2. データの更新頻度が小さい
3. データサイズが大きい

といった特徴を備えたデータセットに向いており、例えば、センサデータをPostgreSQLに収集して解析するといったワークロードに有効な機能。

BRINインデックス自体は、永安さんのこちらの記事が詳しいです。
pgsqldeepdive.blogspot.jp

PG-Stromとしては、搭載RAMが比較的小さなGPUを使うという事もあり、B-treeのようなランダムメモリアクセスを前提としたインデックスへの対応は厳しい。
ただ、条件句の評価はGPUの数千コアを使って並列処理が可能であるものの、インデックスの選択率が高くなると分の悪い勝負なので『このブロックは明らかに該当行なし』という事が分かっているなら、それを読み飛ばしたい。

GpuScanがBRINインデックスを理解し、必要のないブロックを読み飛ばす事ができるようになれば、例えばIoTのキーワードに絡めてセンサデータの集積・解析用途に、という使い方もできるハズである。
特に、PG-Stromはカラムナーを前提としたDWHではないので、生データをそのまま処理させても高速化できるという点は強みになるだろう。

PL/CUDA

PL/CUDAに関しては、言語バインディングに関してする事は（できる事は）多くないので、その周辺領域を拡充していきたい。

一つは、PostgreSQLにおける可変長データの1GB越え。
現状、全ての可変長データの基盤となっている varlena 構造は、最大でも1GBまでのデータしか持てないため、PL/CUDA関数の引数としてarray-matrixを渡す時には、例えば問題領域をうまく分割してデータサイズを1GB未満に抑えてやらないといけない。

昨年10月のCBI学会で発表した研究でも、1000万件の化合物データ（1.3GB）をロードするにはサイズが大きすぎたので、安全マージンも見て4分割した上でGPU側へロードしている。

しかし、昨今のGPUでは10GB近く、あるいはそれ以上のメモリを搭載するのが常識的になりつつあり、問題領域を1GB以内に抑えねばならない、、、というのはユーザにいかにも不都合である。

先にpgsql-hackersにデザインプロポーザルを投げたところ、Robert Haasから『varlenaとは別の体系で可変長データを保持するフォーマットを作成すべき』とサジェスチョンがあり、自分もこの方針には同意。3月のcommit-festまでにはパッチを投稿し、2018年リリース予定のPostgreSQL v11でのマージを目指したい。

もう一つは、現状、複雑な計算ロジックを個々のユーザ毎に書かねばならないという点である。
2017年1月時点でPL/CUDAを実証できたワークロードとしては、k-NN法類似度サーチや、k-meansクラスタリングがあるが、例えば MADLib のような統計解析パッケージで提供されているアルゴリズムの、全部とは言わないまでも、使用頻度が高く計算負荷の高いものをGPUで計算するようパッケージ化できれば、ユーザの裾野はより広がるだろうし、仮にカスタマイズが必要となっても骨格となるアルゴリズムのGPU実装が既に存在する事で省力化が可能となるハズである。

PG-Strom v3.0へ向けた種まき

In-memory Columnar Cache

NVMe-SSDとの密連携の他にI/O系処理を高速化する方策として、通常のPostgreSQLテーブル（行形式）の脇に、予めDMAバッファ上に列形式にデータを再編したキャッシュを持たせる機構を考えている。
このキャッシュはBackground-workerにより非同期で作成され、そのため、スキャンする区間のうち一部領域しか列形式のキャッシュが構築されていないかもしれない。しかし、その場合でもPG-StromはGPUで行形式データを扱えるので（多少のパフォーマンス差に目をつぶれば）大きな問題とはならない。

データ本体とキャッシュを別に持つ場合、必ず一貫性制御が複雑で頭の痛い問題として立ちはだかる。
アイデアとしては、これを避けるために ALL_VISIBLE フラグが 1 であるブロックのみを列形式キャッシュに持たせる。
ALL_VISIBLE=1であるブロックは、MVCC制御に関わらず、全てのレコードが全てのトランザクションから可視である事が保証されている。そのため、複雑な同時実行制御に頭を悩ませる必要はなく、全てのレコードの内容を単純に列データとして展開すればよい。

問題は、PostgreSQL側でテーブルが更新され、ALL_VISIBLE フラグが 0 にクリアされた時のinvalidation処理である。
現状、ここにフックを挟む事はできないので、PostgreSQL側の機能強化を行う必要がある。
デザインプロポーザルを出し、この処理を行うにふさわしい場所とフックの仕様を固めていきたい。

GpuSort(+LIMIT)

PG-Strom v2.0では、実はGpuSort機能の廃止を予定している。これは、ソートという問題の性質上、ある程度問題規模が大きくならないとGPUによる処理時間メリットが出てこない一方で、一度にGPUでソートを行う件数が多くなればなるほど、初期データのローディングに時間がかかるようになり、非同期・多重処理のメリットを得にくいからである。
そのため、GPUでソートを行うセグメントサイズを小さくして単位ローディング時間を短くする一方で、CPUでのMergeSort処理の割合が大きくなるか、それとも、セグメントサイズを大きくするかというジレンマに悩まされてきた。
（で、最終的にはそれほど速くならない事が多かったり・・・。）

根本的な原因は、GPUで処理を行ってもソートではデータ件数を減らす事が不可能な点にある。なので、GpuSortを有効に活かすには、何がしか『データ件数を減らせる』パターンでのソートに限った方が利口だ。
例えば、LIMIT句で『上位xx件を取り出す』という事が明らかな場合に限り、GpuSortを使用するというパターンであれば十分に効果を発揮する事ができるだろう。
これにはPostgreSQL本体側で、『LIMIT句でxx件のデータが要求されますよ』という事を下位のノードに伝えてやる仕組みが必要だが、PostgreSQL v10 向けにパッチを出しており、ある程度レビューも進んでいるので間に合うだろうとは踏んでいる。

こんな感じで2017年の開発ロードマップについて考えてみたが、さてさて、大晦日に振り返ってどの程度きちんとやれているでしょうか。

Beyond the 1GB limitation of varlena

This article is a part of the PostgreSQL Advent Calendar 2016.

According to the request by Joe Conway (@josepheconway), I wrote this article in English.

I like to share the discussion we had at the PostgreSQL developer unconference on the day before PGconf.ASIA 2016, and the related idea of mine.

PostgreSQL supports variable length data types; like text, bytea, json, numeric, array and so on.
These data types individually have their own characteristics and own internal format, however, all of them are built on a common structure to represent variable length field; that is called varlena.

It has very simple internal format. Contents of the variable length fields are followed by either 1-byte or 4-bytes header.

We can identify the header type by the least bit of the first 1 byte. If least bit of the first byte, it means 1-byte header, elsewhere, 4-byte header.
In case of 4-byte header, the second bit is also used to show whether it is compressed or not. Therefore, rest of 30-bits can be available to represent the contents length, so, maximum length of the variable length field is 1GB.
A special case exists if the first byte is 00000001. It is an external TOAST reference which consists of OID of TOAST table and unique ID within TOAST table.

Below is the source code comment at include/postgres.h.

    :
 * Bit layouts for varlena headers on little-endian machines:
 *
 * xxxxxx00 4-byte length word, aligned, uncompressed data (up to 1G)
 * xxxxxx10 4-byte length word, aligned, *compressed* data (up to 1G)
 * 00000001 1-byte length word, unaligned, TOAST pointer
 * xxxxxxx1 1-byte length word, unaligned, uncompressed data (up to 126b)
   :

My concern is we have no way to represent a larger variable length datum more than 1GB.

It is quite natural for users who want to process the maximum available data size as long as system capability allows. It is not uncommon for the recent GPU models to have 10GB class device memory or more.
On the other hands, due to the restriction of variable length datum of PostgreSQL, we cannot have a datum larger than 1GB. It also means we cannot provide a large matrix (= 2D-array) onto PL/CUDA function at once. It is an unignorable restriction if a problem user want to solve by PL/CUDA is unavailable or expensive to split into multiple portions.

According to the background, we discussed a few options to support variable length datum larger than 1GB.

64bits varlena header

The first idea was straightforward but got less interest because 99% of existing variable length data types are satisfied with 1GB limitation.
If we have one more varlena header, VARDATA() and VARSIZE() which are widely used in the PostgreSQL core and extensions need to have branch operation inside the macro, and it is not easy to justify the penalty for this niche usage.

Use of large-object and its ID instead

The second idea was a suggestion from audience. Now PostgreSQL has a feature of large object which allows to store up to 4TB data chunk with a unique identifier. If PL/CUDA function supports ID of large object, instead of varlena datum, PL handler can extract larger data chunk by itself.
However, here is a problematic scenario. It requires users to construct large objects preliminary. It is inconvenient when user wants to deal with a large matrix constructed on the fly. For example, the query below constructs a matrix based on the underlying table scan with qualifier. The qualifier can be changed for each execution.

SELECT * FROM matrix_unnest(
  (SELECT my_plcuda_function(Q.matrix, D.matrix)
     FROM (SELECT array_matrix(a,b,c) matrix
             FROM table_q
            WHERE tag LIKE '%abc%') Q,
          (SELECT array_matrix(x,y,z) matrix
             FROM table_d
            WHERE tag LIKE '%xyz%') D
  )
)

In addition, creation of a large object makes unnecessary i/o; which leads performance slow-down.

Special type that includes indirect reference

A overall consensus was to define a special data type which support indirect references to data chunks less than 1GB. If only narrow use-case wants to have a datum larger than 1GB, it is quite natural to define a special data type for the purpose.

My interest is representation of matrix in database.

In case of matrix, it has less necessity to have all the items in a flat memory chunk.

If we have 8GB of matrix, we can split it into 9 portions to keep individual chunk size less than 1GB.

Then, once we define a special matrix type that consists of small metadata and indirect references to the chunks less than 1GB, it is available to represent a big matrix larger than 1GB as a usual SQL data type.

A remaining problem is serialization/deserialization when the matrix data type is saved/loaded.
Right now, PostgreSQL saves contents portion of the variable length datum onto the storage as is, thus, pointers to reference sub-matrix has to be serialized appropriately, however, we have no infrastructure to manipulate type specific data structure on toast_insert_or_update() which set up physical tuple for INSERT/UPDATE.
Likely, pg_type needs to have an optional callback functions for serialization/deserialization. It shall be a mandatory requirement if data structure has indirect reference to other memory chunks.

I expect the serialization callback will use TOAST relation to save a large but less than 1GB chunks, then put its unique ID instead of the pointers. We will be able to have another advantage in this approach, because all the sub-matrix we have to update are the portion actually updated. If some of sub-matrix were not updated, we don't need to delete old toast pages and insert new ones. It will make a performance benefit than existing flat varlena structure.

The timing for deserialization needs a bit more consideration. Because heap_tuple_fetch_attr handles deserialization of the existing flat varlena, but no type OID is supplied to the function. It is not a good option to change the function prototype because many existing code already uses this function without type OID.
We have two another options here. The first one packs type OID within the serialized structure. It needs to define a new VARTAG_* label to distinct from the existing flat varlena. The second one is delayed load because indirectly referenced data chunks will not be used without functions/operators which support the data type. It enables not to load unreferenced chunk, however, it is uncertain whether functions/operators can manipulate a value supplied as an argument. *1

Last, even if we can have variable length datum larger then 1GB from the viewpoint of data format, it is never small data chunks. It involves not a small amount of i/o (or memory copy) stuff.
Therefore, it is a significant to have special optimization based on the knowledge for usual use case of the types.

For example, some workloads takes sparse matrix which have small amount of non-zero values, but mostly zero. In this case, type may be able to assume empty sub-matrix are all-zero instead of data size reduction.

Diagonal matrix is also often used, in case when valid values are located around the diagonal axis only.

I hope making a proof of the concept patch near future, then have a discussion at pgsql-hackers.

*1:If not acceptable, we may need to load sub-matrix multiple times when a particular matrix object is referenced by multiple functions/operators.

PGconf.SV 2016 and PL/CUDA

I've participated PGconf Silicon Valley 2016 held at the South San Francisco Conference Center.

I could have a presentation here. Its title is PL/CUDA - Fusion of HPC Grade Power with In-Database Analytics.
Its slides are below:
https://www.slideshare.net/kaigai/pgconfsv2016-plcuda/

What is PL/CUDA?

PL is an abbreviation of Procedural Language. In PostgreSQL, we have some PL options like PL/pyhton, PL/R and so on.
This functionality allows to implement SQL functions with different program languages individually.
For example, this SQL function definition is described with "SQL" language.

CREATE OR REPLACE FUNCTION sample(int, int)
RETURNS int
$$
  SELECT $1 + $2;
$$ LANGUAGE 'sql'

We can describe the code block between the $$ marker in SQL-language.

It means we can use other programming language for SQL function definition if another LANGUAGE is given.
So, it is never strange to have CUDA language here.

PL/CUDA is a procedural language which supports CUDA code block for user defined SQL function.
It allows to embed a heavy computing portion within a SQL query, and invokes GPU kernels transparently.

As literal, DBMS is database management system. It implies heavy computing is not a primary workload for DBMS, and it is almost correct. So, people who want to run heavy analytic algorithm usually exports the data from the database for more computing performance. However, here is two downsides. The first is network traffic to export the data set. If DBMS would have enough computing capability, all we have to export is results of the analytic algorithm. The other is loss of SQL flexibility on the post-process of the analytic algorithm. If we could finish the algorithm core in database, it is quite easy to join the result with other tables, to make a summary using GROUP BY or window functions.

PL/CUDA is a solution to solve these problems. The code block in the function declaration shall be put into GPU kernel, then PG-Strom's infrastructure build the GPU kernel on the fly and execute on GPU devices with thousands threads. Its function arguments and results are automatically move to/from GPU devices.

Actually, PL/CUDA changes the concept of PG-Strom a bit. It originally intended to accelerate SQL by GPU transparently, thus existing SQL shall be accelerated with no query changes.
On the other hands, PL/CUDA assumes the core logic of algorithm is implemented manually for each use case.
Why we changed this concept? It is based on user's feedback.
Even if we don't need to describe CUDA code, right now, nobody implements complicated analytic algorithms in SQL, thus, it means user has to write up a new code to process them either SQL or CUDA. I'm not certain whether SQL is suitable language for this purpose, and most of known algorithm assumes procedural languages.
So, we thought writing a CUDA code may be a valuable option in some use cases, especially, processing of advanced algorithms. We call this type of GPU utilization manual optimization mode, in contradistinction to the existing automatic code generation from SQL; automatic optimization mode.

Array-Matrix

One other thing we need to pay attention is data-format and memory bus utilization rate when GPU kernel handles the data-format.
PostgreSQL adopts row-oriented data format. It is an I/O optimal data format for transactional workloads.
On the other hands, it has a downside when we reference individual columns within a record. Because the location of column is not predictable unless reading the record actually, we have to walk on records from the head, thus, it takes larger number of memory access.
It is unignorable penalty when we run a logic on GPU which has relatively smaller memory cache. In case of column-oriented data format, location of the data can be determined uniquely by pointer of the array and thread-id.

In addition, column-oriented data format has an advantage for GPU devices from the standpoint of memory bus utilization rate.
Unit size of GPU's memory transaction is usually 32bytes/256bits. If multiple processor cores requires coalesced memory location simultaneously, these multiple memory accesses are consolidated to a single memory transaction, then, a single memory transaction will load 8 x 32bits variables at once. This scenario pulls up the maximum utilization rate from the memory bus of GPUs; that is much higher than CPU/RAM bandwidth.
When we run a particular calculation by multiple processor cores, it tends to reference same columns simultaneously. Columns are scattered in row-oriented data structure, so it is hard to pull out the maximum utilization rate of memory bus.

However, columnar-format is not always winner. Transformation from row-format to column-format needs extra cost prior to GPU invocation, and it is hard to justify the data transformation cost for usual SQL operations like JOIN, GROUP BY, etc...
One candidate is highly computing intensive workload like advanced analytic algorithm which we will implement using PL/CUDA.

Fortunately, PostgreSQL already has a data-type which is almost equivalent to the column oriented data format. It is array data type, and we call the array data type with no NULLs and consists of fixed-length data type array-matrix.
An aggregate function array_matrix(variadic datatype[]) consumes multiple rows then produce a 2-dimensional array data. Any elements are re-ordered per column, and data location is predictable by the offset and width of the element data type.

We usually call PL/CUDA function with array-matrix arguments.

Usecase.1 - Similarity Search on Drug Discovery

The first usecase I introduced is a similarity search workloads on drug-discovery *1.

First of all, I like to introduce the background of the problem we tackled.
When researcher tried to find out a candidate of new drug, it has to be "active" to the target protein which is related to the target disease. Some chemical compounds are inactive, some other chemical compounds are active but toxicity. Academic papers, daily published by somewhere, report relationship between a target protein and chemical compounds. So, researcher can gather the public data of chemical compounds relevant to a particular protein.

On the other hands, people said the total number of chemical compounds we can obtain as a test reagent and investigate are about 10,000,000. It is too expensive to purchase and test them in laboratory, with brute-force approach. So, researcher is motivated to pick up chemical compounds which have higher probability of active to the target protein.
It is an approach to search "similar" chemical compounds to the known ones; which are already reported in academic papers and etc...

Similarity is to define a distance between two chemical compounds. Euclid distance is one approach usually used to define geometrical distance between two points.
We used a combination of Jaccard index *2 and k-NN (nearest neighbor) method in our trial.

The k-NN method is used to define distance between a set of items and an item. Once we calculate all the distance between items, than picks up the largest k distances and considers its average the representative distance between them.
Jaccard index can be used when characteristics of two items are described in a bitmap. We consider the ratio of bitmap-AND and bitmap-OR as similarity of two items.

According to the logic, if number of query compounds set is Q, we once need to calculate Q similarities for each database item; that is about 10M. Than, these similarity must be sorted for each database item. Thus, total amount of computing order is almost O(DxQ) + O(DxQlogQ).
It is not small. If Q is 1000, we need to handle 10billion combination to process the logic.

The PL/CUDA function we implemented (knn_gpu_similarity()) did nothing special. It once computes the similarity between two chemical compounds by a processor core for each combination, then run the bitonic-sorting to pick up the largest k-similarities.
This PL/CUDA function takes two array-matrix that are constructed from the relation scan on the fly, or statically pre-built.
Once similarity search gets completed, it unnest the matrix to a generic stream of records by matrix_unnest, then processes the results set using usual SQL features (tables JOIN and window function here).

Its performance is quite impressive. We compared to the same logic based on CPU.
In case of the largest query set (1000), it reduces the response time more than x150 times faster.
It may have a power to change the workstyle of researcher. 3,000sec (=50min) is a timescale for batch job. It is too long to wait in from of the workstation. On the other hands, 20sec is scale for interactive operations; which allows try&error on research process.

Usecase.2 - Data Clustring (k-means)

The second usecase is clustering analysis.
k-means is a representative method for non-hierarchical cluster analysis, used for unsupervized learning.

This algorithm tries to get the final clustering with iteration of simple operations.
First, it assigns a cluster for each item randomly as initial state.

Second, makes a centroid for each cluster according to the initial state.

Next, choose the cluster of every items according to the centroid updated on the last step.

Next, update the centroid according to the latest cluster; updated on the last step.

Then, break the iteration if cluster assignment gets converged or reached to the threshold of the iteration.

The data-set we tried to apply k-means() clustering is a collection of city traffic data at Aarhus, Denmark, by the CityPlus project.
iot.ee.surrey.ac.uk

This dataset contains traffic data (car count, avg speed, ...) in 449 observation point all around the city. The total amount of sensor data is 13.5M records.
In this trial, we applied k-means() algorithm on the raw dataset, then make a summary for each observation point. We assumed each observation points belong to the most frequent cluster.
Then, mapped it on the google map using static map APIs.

The result is below:

It looks to me the blue observation points are bypass highway; many car count and high speed.
On the other hands, green line seems to me the way to the downtown, thus many car but slow speed.
I'm not sure the characteristics of the red line, but appears on the beltway frequently.

Because all the clustering calculation is executed inside of the database, it is quite easy to control the input data.
Below is the same diagram when weekdays and weekend.

• weekdays

• weekend

Here is a little difference between them, however, the diagram of weekend categorize the road to the downtown with same color of the highway. So, it has less traffic on the weekend, thus, car speed is higher than the case of weekday.

For the arrangement, all I did is just adding a WHERE clause on the query to kick k-means clustering. Utilization of the SQL flexibility is one of the biggest advantage of PL/CUDA, in addition to the performance benefit.

SELECT report_id, k, c
  FROM (SELECT report_id, k, c,
               row_number() OVER (PARTITION BY report_id
                                  ORDER BY c DESC) rank
          FROM (SELECT report_id, k, count(*) c
                  FROM matrix_unnest(
                        (SELECT gpu_kmeans ( array_matrix(
                                               int4_as_float4(report_id),
                                               avg_measured_time,
                                               avg_speed,
                                               vehicle_count),
                                             5)
                           FROM tr_rawdata
                          WHERE extract('dow' from timestamp) in (0,6)
                        )
                       ) R(report_id int, k int)
                 GROUP BY report_id, k
               ) __summary_1
       ) __summary_2
   WHERE rank = 1;

*1:KaiGai Kohei, Yoshimori Atsushi, Efficient Similarity Search Using Multiple Reference Molecules on PG-Strom architecture, CBI Annual Meeting, 2016

*2:also called Tanimoto index

PL/CUDAでk-means法を実装する

前回のエントリでは、CBI学会で発表を行った、PL/CUDAによる類似化合物の検索について説明した。

今回は、コレとはまた別のワークロードに対する応用という事で、クラスタリング処理のIn-Database実装に挑戦してみた。
トライしてみたのは k-means法 によるクラスタリング。非階層クラスタリングの領域では最も頻繁に使用される（と、言われている）アルゴリズムで、計算量もそこそこ大きい。

k-meansクラスタリングとは

教師なし学習方式の一つで、所与のデータ群を一定数(=k個)のグループに分類するためのアルゴリズムである。
以下のステップを一定回数、またはクラスタに属するデータ群に変化がなくなるまで繰り返す。

１．初期クラスタをランダムに設定する。

２．各クラスタの中心点を計算する。

３．各データ要素の属するクラスタを更新する。各データ要素はクラスタ中心点が最も近傍であるクラスタに割り当てられる。

４．新しいクラスタ定義に基づいて、各クラスタ中心点を再計算する。

５．以下繰り返し・・・。

これをPL/CUDAで実装するとどうなるか？

基本的には上記のアルゴリズムを愚直にCUDAで実装する事になるので、特別な事はしていない。
一回のPL/CUDA関数呼び出しの中で何度も繰り返し処理を行う形になるので、メインのGPU Kernelは1スレッドで起動し、Dynamic Parallelismを使って別のGPU Kernelを起動する形にする方が合理的である。

関数定義の中身に興味がある方は、こちらをどうぞ。
toybox/gpu_kmeans.sql at master · kaigai/toybox · GitHub

In-databaseでクラスタリングする

以下の図はお試し的に10000件のランダムデータを5つのクラスタに分割してみたもの。綺麗に分類されている。

ただ、ランダムデータだけだと面白くないので、来週のPGconf.SVでの発表に向け、そこそこ件数が多くてクラスタリングに適したデータは無いかと探してみたら・・・あった。

iot.ee.surrey.ac.uk

FP7ファンドの支援を受けたCityPulseというプロジェクトが収集しているデータで、デンマークのオーフス(Aarhus)市の自動車通行状況を調査したデータが公開されている。
これは、ある二地点間を通過した自動車の台数や平均速度がタイムスタンプと共に収集されているもので、2014年2月～6月の四ヵ月間のデータが全部で1350万レコード。
これをクラスタリングにかけると、どこの道路とどこの道路の通行状況が似ているのか分かるはず。

と、いう事でやってみた。

まず、PL/CUDAを用いて実装したgpu_kmeans()関数は以下のように定義されている。
第一引数にクラスタ化すべきデータとそのID値をreal[]行列として与え、第二引数にはクラスタ数を与える。

CREATE OR REPLACE FUNCTION
gpu_kmeans(real[],    -- ID + Data Matrix
           int,       -- k-value (number of clusters)
           int = 10,  -- max number of iteration
           int = 1)   -- seed of initial randomness
RETURNS int[]

これをクエリの中で使用すると以下のようになる。
少しGoogle Map Static APIを勉強して、SQLから直接URLを生成するようにしてみた。
ただ、URL長の制限が8KBという事で、線を何百本も引くような地図は作れないため、泣く泣く125件分だけに限定して地図を生成している。

WITH
summary AS (
SELECT report_id, k, c
  FROM (SELECT report_id, k, c,
               row_number() OVER (PARTITION BY report_id
                                  ORDER BY c DESC) rank
          FROM (SELECT report_id, k, count(*) c
                  FROM matrix_unnest(
                          (SELECT gpu_kmeans(array_matrix(
                                               int4_as_float4(report_id),
                                               avg_measured_time,
                                               avg_speed,
                                               vehicle_count),
                                             5)
                             FROM tr_rawdata
                          )
                       ) R(report_id int, k int)
                 GROUP BY report_id, k
               ) __summary_1
       ) __summary_2
   WHERE rank = 1
),
location AS (
SELECT point_1_lat, point_1_lng,
       point_2_lat, point_2_lng,
       CASE k WHEN 1 THEN 'red'
              WHEN 2 THEN 'blue'
              WHEN 3 THEN 'green'
              WHEN 4 THEN 'purple'
              ELSE 'orange'
       END col
  FROM summary s, tr_metadata m
 WHERE s.report_id = m.report_id
),
path_definition AS (
SELECT 'path=color:' || col || '|weight:3|' ||
       point_1_lat::text || ',' || point_1_lng::text || '|' ||
       point_2_lat::text || ',' || point_2_lng::text path_entry
  FROM location
 LIMIT 125 -- becuase of Goole Map API restriction
)
SELECT 'http://maps.google.com/maps/api/staticmap?' ||
       'zoom=11&' ||
       'size=640x480&' ||
       'scale=2&' ||
       string_agg(path_entry, '&') ||
       '&sensor=false'
  FROM path_definition;

$ wget -O map.png "`psql traffic -At -f ~/traffic.sql`"

その結果がこれ。

どうやら、青の区間は高速道路かバイパス道路で、比較的自動車の台数が多くスピードが出ているクラスタに見える。
一方、緑の区間はダウンタウンに向かう道路なので、慢性的な渋滞に悩まされているのだろうか？
赤は少し判断に迷うが、そのどちらでもない、といった所のようである。

また、これは中核となるkmeansクラスタリングをSQLの中に埋め込んでいるので、抽出条件を変えるだけで簡単に母集団を切替える事ができる。
以下の例は、個々のデータのタイムスタンプを用いて、ウィークディ（月曜～金曜）とウィークエンド（土曜、日曜）でそれぞれクラスタリング、図を作ってみたケース。

■平日

■週末

微妙な差ではあるが、週末のケースでは、ダウンタウンに向かう道路がハイウェイと同じ色に分類されている。
平日に比べると交通量が減る分、自動車がスイスイ進むという事だろうか。

閑話休題。今週のPGconf.SVでは、この辺のネタも交えてPL/CUDAメインで話をしてくる事になります。

PL/CUDAによるIn-database Analytics ～創薬におけるワークロードを例として～

やや場違い感が否めないが、今週、ＣＢＩ学会（計算情報化学生物学会）の2016年大会でポスター発表を行ってきた。

発表タイトルは『Efficient Similarity Search Using Multiple Reference Molecules on PG-Strom architecture』というもので、要は、創薬分野で新薬の候補となる化合物を類似度サーチする際に、PG-Stromを用いて高速かつIn-Databaseで行ってしまおうという試みである。

前提となる利用シーンは、例えば次のようなシチュエーションである。
ある疾患には特定のたんぱく質が作用する事が分かっており、また、世界各地で発表される研究論文の中には『このたんぱく質には〇〇〇という化学物質が活性を持っている』というような事が記載されている。こういった組合せは別に１：１という訳ではなく、別の研究論文には化合物△△△が、別の論文には□□□が、程度の多寡はあれ活性を持っていたという事が書かれている訳である。

研究者は、このような研究論文からデータを引っ張ってきて、『タンパク質Ｘに対して活性を持つ化学物質の組』というのを作る。これを便宜上、Query Chemical Compounds (Q)と称する事にする。

一方、薬剤の候補となる化合物の探索空間としては概ね1000万件程度が知られている。これを便宜上、Database Chemical Compounds (D)と称する事にする。これらがどのような化学的特徴を持つのかという事は分かっているものの、目的となるたんぱく質に活性を示すのか否かは個別に調べてみる必要がある。但し、調べるといってもそれなりに費用・時間を要するため、適切なものを効率的に絞り込まねばならない。

ここで使用するのが類似度探索（Similarity Search）である。既に学術論文等で目的のたんぱく質に活性を示す事が報告されている化学物質に"近い"化合物を、1000万件の探索空間から絞り込む事ができれば、関連しそうな化合物に一気に網をかける事ができる。

選択の基準となるのが、化合物同士の"近さ"であるが、幾つかの方法が考えられる。中心点距離や平均値を用いる方法もあるが、今回はk-NN(k-Nearest Neighbor)法を用いた。
これは、Q化合物群とあるD化合物間の距離（類似度スコア）を全て計算した上で、そのうち最も近傍にあるk個の距離の平均値を代表距離とするものである。

この時、類似度スコアの計算量としては概ねO(QxD)のオーダーとなる。
つまり、Q化合物群が1000個、D化合物群が1000万個である場合には、100億通りの組合せに対して類似度スコアを計算する必要があり、そこそこしんどい計算になる。

今回の研究で使用したアプローチは、Q化合物群、D化合物群を共にPostgreSQLのテーブルに格納し、ここから読み出したデータをPL/CUDA関数の引数として与え、この関数の内部で①類似度スコアの計算、②並び替え（ソート）、③上位k件の平均値導出、までを並列に実行してやろうというものである。

PL/CUDAというのはPG-Stromが提供するPostgreSQLのProcedural Languageの一つで、SQL関数の定義部分にCUDAコードを直接記述する事ができる。実行時には、このコードブロックがGPU Kernel関数のテンプレートの中にスポっと挿入され、PG-Stromのインフラを使ってコードをビルド、GPUでの並列実行を行う。
関数の引数や返り値は、本来はRAM<-->GPU間で転送する必要があるが、この辺はSQL関数の定義から機械的に生成できる部分であるので、PG-Stromが全て面倒をみてくれる。今回は、Matrixに相当するデータ構造としてPostgreSQLの2D配列データ構造を使用してデータの受け渡しを行った。

そして、ひとたびk-NN法類似度サーチのロジックをPL/CUDA関数に組み込めば、あとは見慣れたSQLのパズルである。
例えば、以下のようなクエリを用いてfinger_print_queryテーブルから読み込んだQ化合物群、finger_print_10m_matrixテーブルから読み込んだD化合物群を、PG-Stromの提供するarray_matrix()関数を用いて2D配列データに変換し、PL/CUDA関数の引数として与えている。

PL/CUDA関数の返り値もまた2D配列であるが、これも同様にPG-Stromの提供するmatrix_unnest()関数により、一般的なN行のレコードへと戻している。

PREPARE knn_sim_rand_10m_gpu_v2(int)	-- arg1:@k-value
AS
SELECT float4_as_int4(key_id) key_id, similarity
  FROM matrix_unnest((SELECT rbind(knn_gpu_similarity($1,Q.matrix,D.matrix))
                        FROM (SELECT cbind(array_matrix(id),
                                           array_matrix(bitmap)) matrix
                                FROM finger_print_query
                               LIMIT 99999) Q,
                             (SELECT matrix
                              FROM finger_print_10m_matrix) D))
    AS sim(key_id real, similarity real)
ORDER BY similarity DESC
LIMIT 1000;

ではこれで性能はどれ位出るのか？手元の環境で計測してみたところ、同じロジックをCPUで実装((knn_gpu_similarity()と等価なknn_cpu_similarity()関数をCで実装))して実行した場合と比べ、問題サイズが大きい場合では130倍を超える応答時間の改善を確認できた。

Q-size	CPU(E5-2670v3)	GPU(Tesla K20)	GPU(Tesla M40)	高速化率(CPU vs M40)
10	30.25s	13.23s	13.41s	x2.3
50	145.29s	14.49s	13.54s	x10.7
100	295.95s	15.99s	14.04s	x21.1
500	1503.31s	27.11s	17.47s	x86.1
1000	3034.94s	43.88s	22.58s	x134.4

In-Databaseでのデータ解析処理でここまで応答時間が変わってくると、計算の運用に関わる従来の前提も色々と変わってくることになる。

一つ考えられるシナリオとしては、計算性能を担保するために外部のアプリケーションを使用してデータ解析を行っていたような利用シーンで、In-Database処理が十分な計算性能を出すようになった結果、データをエクスポートする必要がなくなるという事。これは、Hadoop方面の人が言うように『データのある場所で計算を行う』アプローチであり、最近の技術トレンドと合致している。

そしてもう一つの副次的な効果として、PL/CUDA関数の結果がSQLのデータ構造として返却されるという事は、データ解析ロジックの中で最もヘビーな計算処理を終えた後のデータに対して、JOIN、Window関数、ORDER BY、GROUP BYといったSQLの構文を使って非常にフレキシブルなデータの加工・後処理が可能になるという事である。

これらは、研究者がデータを多面的・多角的に眺めて新たな知見を得る事に役立つが、従来50分を要していた処理が22秒まで短縮された事と併せ、研究に不可欠なTRY&ERRORに適した環境を提供してくれることだろう。

同期DMAと非同期DMA

おっとっと、やらかしてしまった（但し、良い方に）。

PG-Strom + NVMe-Stromでのパフォーマンス計測の際に、SSDからロードしたデータ以外に、例えばテーブル定義情報や定数パラメータといったSQLの実行に必要な情報は一般的なRAM-to-GPU DMAで転送していたのだけども、ココがうっかり同期DMAになっていたために、本来の性能を発揮できないでいた。

そこで、きちんと非同期DMAを実行できるようにコードを修正し、改めてPG-Strom + NVMe-Stromの実行速度を測り直した数字が以下の通り。じゃん。

ワークロードは変わらず、以下の三種類のクエリを64GB/7億件のテーブルに対して実行した。

• Q1: 比較的シンプルな検索条件を持つスキャン
• Q2: 比較的複雑な検索条件を持つスキャン
• Q3: 文字列マッチ（LIKE句）を持つスキャン

応答時間が概ね42～43secの範囲に収まっているのがお分かりいただけると思う。
これをスループットに直すと、64GB / 43sec = 1524MB/sec のレートでデータを処理できており、Intel SSD 750のカタログスペック 2200MB/s からすると概ね70%程度となる。

応答時間に殆ど差が見られないという事で、これは、GPUで実行するクエリの評価はI/Oよりも短時間で完了するために、非同期DMA転送の裏側に隠ぺいされてしまったと考える事ができる。

CUDAでは非同期で実行する個々のタスク（例えば、RAM=>GPUへのデータ転送、GPU Kernelの実行、など）の順序関係を制御するために、ストリーム（CUstream）という制御構造を持っている。
ある種当然ではあるが、ホスト側から送出されてくるデータを用いて計算しようというGPU Kernelは、少なくともデータ転送が終わらなければ実行できないし、計算結果をデバイス側⇒ホスト側に転送する時も、GPU Kernelの実行が終わっていないと、計算途中のGPU RAMのイメージを送りつけられても困惑する限りである。

CUDAの持つRAM=>GPUのデータ転送用APIには二種類あり、一つは cuMemcpyHtoD() などの同期API。これは、関数が呼び出された時点で呼び出し元をブロックし、RAM-to-GPU DMAを用いてデータを転送する。関数から戻った時点で、既にGPU側のバッファにデータの転送は終わっており、ストリームとは無関係に使用できる。
もう一つは cuMemcpyHtoDAsync() などの非同期API。これは、ストリームにDMA要求が突っ込まれた順番に、非同期にデータ転送を行う。GPU Kernelの実行開始なども、ストリームに要求が突っ込まれた順序を元にした依存関係を壊さないように、ただし、DMAチャネルなどのリソースが空けば待っているタスクはどんどん実行される。

ただし、cuMemcpyHtoDAsync()を用いても必ずしも非同期DMAになるかと言えば、少々落とし穴があり、CPU側のRAMを『ここはDMAバッファだからスワップアウトしないでね』とCUDAランタイムに登録した領域以外を転送元/転送先に指定した場合、黙って同期DMAモードで動作するのである。
今回の場合がそれで、本来はDMAバッファを cuMemHostRegister()で登録しておかねばならないところを忘れていた。Orz

結果、本来であれば次々とデータ転送を行えるところが、ストリームに突っ込んだ cuMemcpyHtoDAsync() が実行可能になるまでブロックされた挙句、同期DMAを行ったものだから、トータルの処理時間が随分と間延びした形になっていた。あーあ。

まぁ、スコアとしては前に計測した時から50%ほど伸びて、対PostgreSQLで見てみると、３～４倍程度のスループットを発揮しているので善しとしよう。