GPGPU Advent Calendar 4日目です。
論文紹介ネタでもいいとのことなのでCPU比100倍とか景気の良いことを言ってる拙著論文の紹介を…、というのは誰も得しない感しかないので、ここ1年ぐらいのPyOpenCLのアップデート状況についての話をしましょう。

PyOpenCLとは

PyOpenCLはPythonからOpenCLのAPIにアクセスするライブラリで、APIのラッパーを提供するだけでなく、スクリプト言語らしい便利な機能が多数実装されています。

同じ作者が開発したPyCUDAというCUDA用の同様のライブラリもあるのですが、私は最近ではずっとPyOpenCLの方を利用しています。その主な理由としては、NVIDIA GPU以外の環境でも実行ができることと、カーネル関数の実行時コンパイルの仕組みがOpenCLのAPIレベルで提供されているのでライブラリのソースが読みやすいことが挙げられます。また、OpenCL単体では機能的に貧弱な部分をPythonで補えるのは魅力的です。

PyCUDAとPyOpenCLのどちらを利用するべきかという話が時々ありますが、性能的にはあまり違いはないので、CUDAかOpenCLの慣れている方で選ぶか、PyCUDA / PyOpenCLでテスト実装したカーネル関数を将来的にC/C++に移植するなどであれば、ターゲットにするプラットフォームにあらかじめ合わせておいた方がいいでしょう。

PyOpenCLのアップデート状況

PyOpenCLはバージョン0.92以降、2011.1、2011.2、2012.1というナンバリングでリリースが行われていて、現在はバージョン2012.2リリースに向けて開発が行われています。主要な追加機能の実装は落ち着いてるようなので、今年中かはわかりませんが遠くないうちにリリースされるんではないでしょうか。

さて、本題のここ1年でPyOpenCLに追加された機能についてお話しましょう。
大きい変更点は以下の2点です。

複素数サポート及びベッセル関数の追加 (2012.1)
Parallel Scanのインタフェース変更とそれに伴う各種アルゴリズムの追加 (2012.2)

ベッセル関数という言葉はPyOpenCLのリリースノートで初めて見たぐらい詳しくないので、今回は刷新されたParallel Scanとそれをベースにしたアルゴリズムについて紹介します。

Parallel Scan

Scanとは、Pythonコードで大雑把に説明すると以下のような処理になります。

from functools import reduce


def add(a, b): 
    return a + b 


def inclusive_scan(f, seq):
    return [reduce(f, seq[0:i + 1]) for i in range(len(seq))]


def exclusive_scan(f, seq, init):
    return [reduce(f, seq[0:i], init) for i in range(len(seq))]


seq = [1, 3, 5, 7, 9]

print(inclusive_scan(add, seq))
# [1, 4, 9, 16, 25]

print(exclusive_scan(add, seq, 0)) 
# [0, 1, 4, 9, 16]

以上のように、リストを与えると、各要素が自分より手前の要素までの総和を持ったリストが得られます。自分自身を含むかどうかでinclusive_scan、exclusive_scanの二種類があり、exclusive_scanの場合は単位元が必要です。

Scanは和以外にも、ある条件を満たす二項演算に対して適用可能で、例えば自分より前の要素までの積や最小値、最大値を求めることなどができます。並列的なScanが可能な演算の条件については以下のブログ記事で議論されています。

並列環境におけるreduceとscanアルゴリズム

Scanをベースにしたアルゴリズム

並列Scanアルゴリズムは様々な実装があり、よくチューニングされたものがCUDAやOpenCLのサンプルコードとして添付されていたりします。それでは、Scanが効率よくできるとどう嬉しいのでしょうか？

それは、以下に挙げるようなstream compaction系のアルゴリズムやソートまでもがScanをもとに実装することができるからです。新しく導入されたアルゴリズムについて実際のPyOpenCLを使ったコードで見てみましょう。

from __future__ import print_function
import numpy as np
import pyopencl as cl
from pyopencl import array as clarray
from pyopencl import algorithm as clalgorithm

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

size = 100
host_a = np.random.randint(0, 10, size).astype(np.int32)
print(host_a)
# [1 1 1 6 6 9 8 9 8 1 2 7 0 6 1 6 5 4 4 5 3 5 1 0 3 6 9 4 8 5 1 4 8 8 1 3 8
#  3 7 5 1 7 1 7 6 1 9 9 7 3 4 1 8 2 8 3 0 6 4 9 0 2 7 4 4 7 7 1 3 4 7 7 9 6
#  6 7 9 1 2 3 1 2 3 9 1 7 3 0 9 5 8 1 5 4 5 2 8 8 1 9]
a = clarray.to_device(queue, host_a)


# 1. 条件に一致する要素のみを残す（arrayの前側に寄せる)
b, b_size = clalgorithm.copy_if(a, 'ary[i] % 2 == 0')
host_b = b.get()[:b_size.get()]
print(host_b)
# [6 6 8 8 2 0 6 6 4 4 0 6 4 8 4 8 8 8 6 4 8 2 8 0 6 4 0 2 4 4 4 6 6 2 2 0 8
#  4 2 8 8]


# 2. 条件に一致する要素を取り除く(一致しない要素をarrayの前側に寄せる)
c, c_size = clalgorithm.remove_if(a, 'ary[i] % 2 == 0')
host_c = c.get()[:c_size.get()]
print(host_c)
# [1 1 1 9 9 1 7 1 5 5 3 5 1 3 9 5 1 1 3 3 7 5 1 7 1 7 1 9 9 7 3 1 3 9 7 7 7
#  1 3 7 7 9 7 9 1 3 1 3 9 1 7 3 9 5 1 5 5 1 9]


# 3. 条件に一致する要素を前側に、そうでない要素を後側に寄せる
d, e, d_size = clalgorithm.partition(a, 'ary[i] % 2 == 0')
host_d = d.get()[:d_size.get()]
host_e = e.get()[:size - d_size.get()]
print(host_d)
# [6 6 8 8 2 0 6 6 4 4 0 6 4 8 4 8 8 8 6 4 8 2 8 0 6 4 0 2 4 4 4 6 6 2 2 0 8
#  4 2 8 8]
print(host_e)
# [1 1 1 9 9 1 7 1 5 5 3 5 1 3 9 5 1 1 3 3 7 5 1 7 1 7 1 9 9 7 3 1 3 9 7 7 7
#  1 3 7 7 9 7 9 1 3 1 3 9 1 7 3 9 5 1 5 5 1 9]
# partitionの各結果はcopy_if, remove_if個別の結果と等しい
print(np.all(host_b == host_d))
# True
print(np.all(host_c == host_e))
# True


# 4. ソート
sort_kernel = clalgorithm.RadixSort(ctx, 'int* ary', 'ary[i]', ['ary'])
(f,) = sort_kernel(a, key_bits=32)
host_f = f.get()
print(host_f)
# [0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 3 3
#  3 3 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 6 6 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7
#  7 7 7 7 8 8 8 8 8 8 8 8 8 8 8 9 9 9 9 9 9 9 9 9 9 9]


# 5. 連続する重複要素を取り除く
g, g_size = clalgorithm.unique(f)
host_g = g.get()[:g_size.get()]
print(host_g)
# [0 1 2 3 4 5 6 7 8 9]

以上のようなアルゴリズムを使いこなすことができればGPU上で効率よく実行できる処理の幅は広がるでしょう。

最後に

今回紹介したScanをベースにした新しいアルゴリズムは、現在のリリースバージョン(2012.1)には含まれておらず、リポジトリの最新版を入手する必要があります。

PyOpenCLは如何せん、様々なOpenCLプラットフォームやPythonのバージョンに対応しなければならないため、バグには出会いやすいプロダクトでもあるかと思います。機能追加だけでなく、バグ修正も頻繁に行われているので、githubから最新版を試してみるのをお勧めします。

第2弾があるかは不明。