PL/CUDAとmatrix型
PG-Stromには↓のような利点があるが、
- SQLから自動的にGPU用バイナリ命令列を生成するため、GPUプログラミングを意識する必要がない、
- PostgreSQLの行指向データを用いるため、既存DBからデータの移行を必要としない。
その裏返しとして、同時に↓のような特徴も持っている。
- NULL値チェックや四則演算ごとのオーバーフローチェックなどSQLに由来する動作のため、専用に設計されたGPUプログラムほど高速には動作できない。
- 行指向であるためデータの密度が低く、目的のデータにアクセスするために何度もRAMアクセスを繰り返す必要がある。
比較対象がSQLであれば、これらの欠点には目をつぶっても、数百～数千コアの計算能力で殴りつけて十分なアクセラレーション性能を得る事ができる。（というか、GPUにデータを供給する足回りの方が先に悲鳴を上げる）
が、例えばRなどのように、専用の統計解析パッケージに対して計算能力で優位に立つ事を考えると、計算の中核部分などはSQL由来のアレやコレやは省く事で、GPU本来のパフォーマンスを引き出してやりたいところである。
そこで実装してみたのが、PL/CUDA言語とmatrix型。
PL/CUDAはSQL関数を実装するための言語で、CREATE FUNCTION構文で定義したSQL関数の定義部分をそのままCUDAのkernel関数としてビルド。関数の引数を自動的にGPUへ送出してユーザ定義の処理ロジックを並列実行し、実行結果をまたCPU側へ書き戻す。
引数にはfloatやtextなどが利用できるが、新たにmatrix型というデータ型を定義した。これの実体は単精度浮動小数点型の2次元配列であるが、全ての値が非NULLである事を保証しているので、行列の(i,j)要素をダイレクトにアクセスできる。要は、計算密度の高い部分には密度の高いデータを供給するようにしようという話である。
試しにCREATE FUNCTIONでPL/CUDA関数を実装してみる。簡単な例という事で、行列同士の乗算を行うGPU関数を実装してみた。
CREATE OR REPLACE FUNCTION matrix_gpu_mul(matrix, matrix)
RETURNS matrix
AS $$
#plcuda_prep
#plcuda_num_threads 1
if (!pg_matrix_sanitychecks(kcxt, arg1) ||
!pg_matrix_sanitychecks(kcxt, arg2))
PLCUDA_ERROR_RETURN(StromError_InvalidValue);
MatrixType *X = (MatrixType *)arg1.value;
MatrixType *Y = (MatrixType *)arg2.value;
MatrixType *R;
size_t len;
if (X->width != Y->height)
PLCUDA_ERROR_RETURN(StromError_InvalidValue);
len = sizeof(cl_float) * (size_t)X->height * (size_t)Y->width;
if (kplcuda->results_bufsz < len)
PLCUDA_ERROR_RETURN(StromError_DataStoreNoSpace);
retval->isnull = false;
retval->value = (varlena *)results;
R = (MatrixType *) results;
pg_matrix_init_fields(R, X->height, Y->width);
#plcuda_begin
#plcuda_num_threads matrix_gpu_mul_num_threads
MatrixType *X = (MatrixType *)arg1.value;
MatrixType *Y = (MatrixType *)arg2.value;
MatrixType *R = (MatrixType *)retval->value;
cl_uint x_height = X->height;
cl_uint y_height = Y->height;
cl_float *xval;
cl_float *yval;
cl_uint i, j, k, nloops;
assert(X->width == Y->height);
nloops = X->width;
i = get_global_id() % x_height;
j = get_global_id() / x_height;
xval = X->values + i;
yval = Y->values + j * y_height;
if (get_global_id() < x_height * y_height)
{
cl_float sum = 0.0;
for (k=0; k < nloops; k++)
{
sum += (*xval) * (*yval);
xval += x_height;
yval ++;
}
R->values[get_global_id()] = sum;
}
#plcuda_end
#plcuda_results_bufsz matrix_gpu_mul_results_bufsz
$$ LANGUAGE plcuda;何ヶ所か#plcuda_で始まる見慣れないディレクティブがあるが、これがPL/CUDAのディレクティブで、#plcuda_prepから始まるコードブロック、#plcuda_beginから始まるコードブロックがそれぞれ、GPU kernel関数の実体として差し込まれる事になる。
基本的に、#plcuda_prepブロックはバッファの初期化を、#plcuda_beginはアルゴリズムの中核部分の処理を、そしてここでは使っていないが#plcuda_postブロックは結果バッファへの書き戻しと、3ステップでGPU kernel関数を実行する事を想定している。
ただ、Dynamic Parallelismを使って他のGPU kernel関数をバンバン起動する事もできるので、別にこの3ステップに拘る必要はない。
コードブロックの中で指定されている#plcuda_num_threadsは、GPU kernel関数を起動する時のスレッド数を指定する。バッファ初期化を行う#plcuda_prepブロックは1スレッドで起動され、行列積の計算を行う#plcuda_begin...#plcuda_endブロックは、別のSQL関数であるmatrix_gpu_mul_num_threadsの返り値と同じ数のスレッドで起動される。
行列積の大きさは入力された行列の大きさによって変わるので、予め特定の値で決め打ちする事ができないためである。
末尾で指定している#plcuda_results_bufszは、結果バッファの大きさを指定する。これは定数でもよいが、行列積のように出力結果が入力行列の大きさによって変わる場合には、別のSQL関数を呼びだして結果バッファを動的に決められるようにすべきである。
で、これを実行した結果、以下のような計算結果が得られる。
パッと見た感じ、単純な3x3行列同士の演算であるが、裏方ではGPU用のコードをビルドし、引数の行列を2つGPUへ転送して計算を行い、結果をCPU側へ書き戻している。
postgres=# select matrix_gpu_mul('{{2,0,0},{0,2,0},{0,0,2}}'::matrix,
'{{1,2,3},{4,5,6},{7,8,9}}'::matrix);
matrix_gpu_mul
--------------------------------
{{2,4,6},{8,10,12},{14,16,18}}
(1 row)plcuda_function_source関数を使えば、ここで定義したPL/CUDA関数がどのようなGPU kernel関数に置き換えられ、ビルドされているのかを確認する事ができる。
以下のコード中、/* ---- code by pl/cuda function ---- */で挟まれた部分が、PL/CUDA関数の定義から差し込まれた部分である。
postgres=# select plcuda_function_source('matrix_gpu_mul'::regproc);
plcuda_function_source
---------------------------------------------------------------------
#include "cuda_common.h" +
#include "cuda_matrix.h" +
#include "cuda_plcuda.h" +
+
STATIC_INLINE(void) +
__plcuda_matrix_gpu_mul_prep(kern_plcuda *kplcuda, +
void *workbuf, +
void *results, +
kern_context *kcxt) +
{ +
pg_matrix_t *retval __attribute__ ((unused)); +
pg_matrix_t arg1 __attribute__((unused)); +
pg_matrix_t arg2 __attribute__((unused)); +
assert(sizeof(*retval) <= sizeof(kplcuda->__retval)); +
retval = (pg_matrix_t *)kplcuda->__retval; +
assert(retval->isnull || (void *)retval->value == results); +
arg1 = pg_matrix_param(kcxt,0); +
arg2 = pg_matrix_param(kcxt,1); +
+
/* ---- code by pl/cuda function ---- */ +
if (!pg_matrix_sanitychecks(kcxt, arg1) || +
!pg_matrix_sanitychecks(kcxt, arg2)) +
PLCUDA_ERROR_RETURN(StromError_InvalidValue); +
MatrixType *X = (MatrixType *)arg1.value; +
MatrixType *Y = (MatrixType *)arg2.value; +
MatrixType *R; +
size_t len; +
if (X->width != Y->height) +
PLCUDA_ERROR_RETURN(StromError_InvalidValue); +
len = sizeof(cl_float) * (size_t)X->height * (size_t)Y->width; +
if (kplcuda->results_bufsz < len) +
PLCUDA_ERROR_RETURN(StromError_DataStoreNoSpace); +
retval->isnull = false; +
retval->value = (varlena *)results; +
R = (MatrixType *) results; +
pg_matrix_init_fields(R, X->height, Y->width); +
/* ---- code by pl/cuda function ---- */ +
} +
+
KERNEL_FUNCTION(void) +
plcuda_matrix_gpu_mul_prep(kern_plcuda *kplcuda, +
void *workbuf, +
void *results) +
{ +
kern_parambuf *kparams = KERN_PLCUDA_PARAMBUF(kplcuda); +
kern_context kcxt; +
+
assert(kplcuda->nargs == kparams->nparams); +
INIT_KERNEL_CONTEXT(&kcxt,plcuda_prep_kernel,kparams); +
__plcuda_matrix_gpu_mul_prep(kplcuda, workbuf, results, &kcxt); +
kern_writeback_error_status(&kplcuda->kerror_prep, kcxt.e); +
} +
+
STATIC_INLINE(void) +
__plcuda_matrix_gpu_mul_main(kern_plcuda *kplcuda, +
void *workbuf, +
void *results, +
kern_context *kcxt) +
{ +
pg_matrix_t *retval __attribute__ ((unused)); +
pg_matrix_t arg1 __attribute__((unused)); +
pg_matrix_t arg2 __attribute__((unused)); +
assert(sizeof(*retval) <= sizeof(kplcuda->__retval)); +
retval = (pg_matrix_t *)kplcuda->__retval; +
assert(retval->isnull || (void *)retval->value == results); +
arg1 = pg_matrix_param(kcxt,0); +
arg2 = pg_matrix_param(kcxt,1); +
+
/* ---- code by pl/cuda function ---- */ +
MatrixType *X = (MatrixType *)arg1.value; +
MatrixType *Y = (MatrixType *)arg2.value; +
MatrixType *R = (MatrixType *)retval->value; +
cl_uint x_height = X->height; +
cl_uint y_height = Y->height; +
cl_float *xval; +
cl_float *yval; +
cl_uint i, j, k, nloops; +
assert(X->width == Y->height); +
nloops = X->width; +
i = get_global_id() % x_height; +
j = get_global_id() / x_height; +
xval = X->values + i; +
yval = Y->values + j * y_height; +
if (get_global_id() < x_height * y_height) +
{ +
cl_float sum = 0.0; +
for (k=0; k < nloops; k++) +
{ +
sum += (*xval) * (*yval); +
xval += x_height; +
yval ++; +
} +
R->values[get_global_id()] = sum; +
} +
/* ---- code by pl/cuda function ---- */ +
} +
+
KERNEL_FUNCTION(void) +
plcuda_matrix_gpu_mul_main(kern_plcuda *kplcuda, +
void *workbuf, +
void *results) +
{ +
kern_parambuf *kparams = KERN_PLCUDA_PARAMBUF(kplcuda); +
kern_context kcxt; +
+
assert(kplcuda->nargs == kparams->nparams); +
if (kplcuda->kerror_prep.errcode != StromError_Success) +
kcxt.e = kplcuda->kerror_prep; +
else +
{ +
INIT_KERNEL_CONTEXT(&kcxt,plcuda_main_kernel,kparams); +
__plcuda_matrix_gpu_mul_main(kplcuda, workbuf, results, &kcxt);+
} +
kern_writeback_error_status(&kplcuda->kerror_main, kcxt.e); +
} +
+
(1 row)
今回は、まずPL/CUDA関数を定義し、実行できるようになったところまで。
次はもう少し実際的なケースで、計算量の多いものをPL/CUDA関数で実装する例をご紹介したい。
