pyopencl.arrayの使い方

PyOpenCLはPythonからOpenCLAPIにアクセスするためのライブラリですが、それだけではなく素のOpenCLにはない機能もいくつか追加されています。pyopencl.arrayは、numpy.ndarrayライクのインタフェースでOpenCLのデバイスを意識することなく基本的な行列演算を行うことができます。しかしながら、PyOpenCLのドキュメントからは使い方がわかりにくい面もありますので、今回はその基本的な使い方について紹介します。numpyの基本的な知識を要求する気がします。

Arrayオブジェクトの作成

Arrayオブジェクトを作る方法は大きくわけて3つあります。

  1. hostのデータをもとに作成する
  2. device上で生成する
  3. 既にあるArrayオブジェクトから作成する

それぞれの方法についてソースコードを例示します。

import pyopencl as cl
from pyopencl import array as clarray
from pyopencl import clmath
from pyopencl import clrandom
import numpy

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

# 1. hostのデータをもとに作成する

host_a = numpy.random.rand(100)
a = clarray.to_device(queue, host_a)

# 2. device上で生成する

# 空の浮動小数点数型Arrayオブジェクト
b = clarray.empty(queue, 1000, numpy.float32)
# 0で初期化された2次元整数型Arrayオブジェクト
c = clarray.zeros(queue, (10, 10), numpy.int32)
# range
d = clarray.arange(queue, 0, 100, 10, dtype=numpy.int32)
# 乱数生成
e = clrandom.rand(queue, 100, numpy.float32)

# 3. 既にあるArrayオブジェクトをもとに作成する

f = clarray.empty_like(a)
g = clarray.zeros_like(a)

なお、emptyやempty_likeで作ったArrayオブジェクトはfillメソッドやsetメソッド、clrandom.fill_randを使って初期化することができます。

b.fill(2.5)

clrandom.fill_rand(f, a=-10, b=10)

h = clarray.empty_like(a)
h.set(host_a)

上記のように、randとfill_rand関数に引数a, bを与えると乱数の生成範囲が[a, b)になります。デフォルト値は[0, 1)です。

Arrayオブジェクトの演算

Arrayオブジェクトはnumpy.ndarrayと同様に演算をすることができます。

スカラー値との演算

a = clrandom.rand(queue, 5, numpy.int32, a=-10, b=10)
print a.get()

print a + 1
print a - 1.5
print a * 3
print a / 2.0
print a ** 2.0

getメソッドはdevice上にあるArrayオブジェクトの内容をhost側のndarrayオブジェクトにコピーします。Arrayオブジェクトの__str__メソッド、__repr__メソッドの中ではgetメソッドが使われており、printなどで表示するときには自動的にhost側に転送されるため、本記事では以降省略するものとします。

実行結果

[ 4 -8 -6 -8  2]  
[ 5 -7 -5 -7  3]  
[ 2.5 -9.5 -7.5 -9.5  0.5] 
[ 12 -24 -18 -24   6] 
[ 2. -4. -3. -4.  1.]
[ 16.  64.  36.  64.   4.]

Arrayオブジェクト同士の演算

b = clrandom.rand(queue, 5, numpy.int32, a=-10, b=10)
print b

print a + b
print a - b
print a * b
print a / b
print a ** b.astype(numpy.float32)

実行結果

[ 9  7 -6  2 -8] 
[ 13  -1 -12  -6  -6]  
[ -5 -15   0 -10  10]
[ 36 -56  36 -16 -16]
[ 0 -1  1 -4  0] 
[  2.62144000e+05  -2.09715200e+06   2.14334705e-05   6.40000000e+01
   3.90625000e-03]

単項演算

print -a
print abs(a)
print len(a)

実行結果

[-4  8  6  8 -2] 
[4 8 6 8 2]
5

numpy.ndarrayでは他にもビット演算などがサポートされていますが、PyOpenCLでは今のところサポートされていません。

数学関数

Arrayオブジェクトの各要素に数学関数などOpenCLでサポートされた関数を適用するためにはpyopencl.clmath以下の関数を使います。OpenCLの仕様上Arrayオブジェクトのdtypeは浮動小数点数型である必要があります。利用可能な関数についてはPyOpenCLのドキュメントOpenCLの仕様書もご覧ください。

a = clrandom.rand(queue, 5, numpy.float32)
print a
     
print clmath.sqrt(a)
print clmath.sin(a)
print clmath.exp(a)

実行結果

[ 0.45427519  0.19955009  0.48750299  0.80684304  0.14642352]
[ 0.67399943  0.44671029  0.69821417  0.89824444  0.38265327]
[ 0.43881115  0.19822837  0.46842125  0.72210687  0.14590086]
[ 1.5750314   1.22085333  1.62824535  2.24082255  1.15768635]

Reduction

配列の総和を求めたり、最小値、最大値を求めるような処理はリダクションと呼ばれます。GPUでこのようなリダクション処理を効率よく行わせる場合には数々のチューニングテクニックが必要で、結構な量のコードを書く必要がありますが、PyOpenCLではこのようなショートカットが用意されているので詳細を気にせずに高速な処理を行わせることができます。

a = clrandom.rand(queue, 1000, numpy.float32)
b = clrandom.rand(queue, 1000, numpy.float32)

print clarray.sum(a)
print clarray.min(a)
print clarray.max(a)
print clarray.dot(a, b)

実行結果

510.130096436          
0.00353157520294
0.998835802078         
251.52911377

Subset

Arrayオブジェクトから指定したインデックスの要素を取り出した新しいArrayオブジェクトを作ったり、そのような部分配列からリダクションを行う機能が用意されています。

a = clrandom.rand(queue, 1000, numpy.float32)
b = clrandom.rand(queue, 1000, numpy.float32)
indices = clarray.arange(queue, 0, 1000, 2, dtype=numpy.int32)

c = clarray.take(a, indices)
print clarray.subset_dot(indices, a, b)
print clarray.subset_min(indices, a)
print clarray.subset_max(indices, a)

実行結果

120.567108154
0.00146424770355
0.999368607998

Conditionals

以下の関数は2つのArrayオブジェクトの各インデックスにおいて条件にあう要素を取り出します。if_positiveでは第1引数の要素が0以下の場合第3引数の要素が、そうでない場合第2引数の要素が選ばれます。

a = clrandom.rand(queue, 5, numpy.float32)
b = clrandom.rand(queue, 5, numpy.float32)
criterion = clrandom.rand(queue, 5, numpy.int32, b=2)
print a
print b
print criterion

print clarray.if_positive(criterion, a, b)
print clarray.maximum(a, b)
print clarray.minimum(a, b)

実行結果

[ 0.45082718  0.26803267  0.16043031  0.3467738   0.0948478 ]
[ 0.35381329  0.95943987  0.29650748  0.53794491  0.68245429]
[0 1 0 0 0]
[ 0.35381329  0.26803267  0.29650748  0.53794491  0.68245429]
[ 0.45082718  0.95943987  0.29650748  0.53794491  0.68245429]
[ 0.35381329  0.26803267  0.16043031  0.3467738   0.0948478 ]

カスタムKernel

先述のリダクションのように効率よく並列できるパターンというのがいくつかあります。PyOpenCLではそういったパターンの内Elementwise(transform)、Reduction、Scanについてメタプログラミングの仕組みが提供されています。

from pyopencl.elementwise import ElementwiseKernel
from pyopencl.reduction import ReductionKernel
from pyopencl.scan import ExclusiveScanKernel
from pyopencl.scan import InclusiveScanKernel
    
a = clrandom.rand(queue, 5, numpy.float32)
print a

k1 = ElementwiseKernel(ctx, 'float* a', 'a[i] = sin(a[i]);')
b = a * 1
k1(b)                        
print b
                             
k2 = ReductionKernel(ctx, numpy.float32, neutral='0', reduce_expr='a + b', map_expr='x[i]', arguments='__global float* x')
print k2(a)

k3 = ExclusiveScanKernel(ctx, numpy.float32, scan_expr='a + b', neutral='0')
print k3(a * 1)
           
k4 = InclusiveScanKernel(ctx, numpy.float32, scan_expr='a + b')
c = clarray.empty_like(a)
print k4(a, c)

実行結果

[ 0.48551953  0.27677017  0.63273305  0.30237836  0.28056026]
[ 0.46666792  0.27325016  0.59135091  0.29779151  0.27689403]
1.9779613018
[ 0.          0.48551953  0.7622897   1.39502275  1.69740105]
[ 0.48551953  0.7622897   1.39502275  1.69740105  1.9779613 ]

ここではPyOpenCLについて解説しましたが、PyCUDAでもgpu_arrayという同等の機能が提供されていて、ほぼ同じように利用することが可能です。