cudaを使ったparallel programming(並列計算)

cuda programmingには、pycudaの他にnumba.cuda.jitも利用できる。numbaの方がチュートリアルやサンプルコードがネット上に多ゴロゴロしているので、むしろnumbaの方が学習しやすいというメリットがある。ただ、スピード的にはpycudaの方が圧倒的に速いので、そこのところはデメリットと言えるが、cudaプログラミングの学習には参考資料が多いので都合が良い。

import numba.cuda
import numpy as np
import math

my_gpu = numba.cuda.get_current_device()
print("Running on GPU:", my_gpu.name)
cores_per_capability = {
    1: 8,
    2: 32,
    3: 192,
    5: 128,
    6: 128,
    7: 64
}
cc = my_gpu.compute_capability
print("Compute capability: ", "%d.%d" % cc, "(Numba requires >= 2.0)")
majorcc = cc[0]
print("Number of streaming multiprocessor:", my_gpu.MULTIPROCESSOR_COUNT)
cores_per_multiprocessor = cores_per_capability[majorcc]
print("Number of cores per mutliprocessor:", cores_per_multiprocessor)
total_cores = cores_per_multiprocessor * my_gpu.MULTIPROCESSOR_COUNT
print("Number of cores on GPU:", total_cores)
Running on GPU: b'GeForce GTX 1060 with Max-Q Design'
Compute capability:  6.1 (Numba requires >= 2.0)
Number of streaming multiprocessor: 10
Number of cores per mutliprocessor: 128
Number of cores on GPU: 1280

コア数がたった1280しかないのは知らなかった。Quadro GV100は、cudaコアが5120にtensorコアが640、HBM2が32Gという構成で価格が150万円らしい。高いんだか安いんだか分からん。

CUDA Programmingを始める前に

先ずは恒例のmodule import等の事前準備

from __future__ import division
import os
import sys
import glob
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from IPython.display import HTML, Image
from timeit import default_timer as timer
%matplotlib inline
%precision 4
plt.style.use('ggplot')
from numba import cuda, vectorize, guvectorize
from numba import void, uint8 , uint32, uint64, int32, int64, float32, float64, f8
import numpy as np

Cuda Programmingのステップ

カーネル関数によるGPU上のベクトルの加算から始める。段階は以下の通り。

  1. カーネル関数(GPU上で並列計算を実行するためのコード)を定義する。最も単純なモデルは一つのカーネルが同時に命令を処理してCPUへのリターンをコントロールする。
  2. 多くのスレッドが一つのカーネルを実行する。
  3. CPU上に加算されるベクトルと加算されたベクトルの場所を割り当てる。
  4. ベクトルをGPUにコピーする。
  5. カーネルを与えられた変数(grid and block dimension)で実行する。
  6. 加算したベクトルをコピーしてCPUに返す。

カーネルバージョン1を作成

このバージョンはCUDA Cで実行される必要がある事を基本的にきっちりとこなせる。

@cuda.jit('void(float32[:], float32[:], float32[:])')
def cu_add1(a, b, c):
    """This kernel function will be executed by a thread."""
    bx = cuda.blockIdx.x # which block in the grid?
    bw = cuda.blockDim.x # what is the size of a block?
    tx = cuda.threadIdx.x # unique thread ID within a blcok
    i = tx + bx * bw

    if i > c.size:
        return

    c[i] = a[i] + b[i]

カーネル Ver.1を実行する

device = cuda.get_current_device()

n = 100

# Host memory
a = np.arange(n, dtype=np.float32)
b = np.arange(n, dtype=np.float32)

# Assign equivalent storage on device
da = cuda.to_device(a)
db = cuda.to_device(b)

# Assign storage on device for output
dc = cuda.device_array_like(a)

# Set up enough threads for kernel
tpb = device.WARP_SIZE
bpg = int(np.ceil(float(n)/tpb))
print ('Blocks per grid:', bpg)
print ('Threads per block', tpb)

# Launch kernel
cu_add1[bpg, tpb](da, db, dc)

# Transfer output from device to host
c = dc.copy_to_host()

print (c)
Blocks per grid: 4
Threads per block 32
[  0.   2.   4.   6.   8.  10.  12.  14.  16.  18.  20.  22.  24.  26.
  28.  30.  32.  34.  36.  38.  40.  42.  44.  46.  48.  50.  52.  54.
  56.  58.  60.  62.  64.  66.  68.  70.  72.  74.  76.  78.  80.  82.
  84.  86.  88.  90.  92.  94.  96.  98. 100. 102. 104. 106. 108. 110.
 112. 114. 116. 118. 120. 122. 124. 126. 128. 130. 132. 134. 136. 138.
 140. 142. 144. 146. 148. 150. 152. 154. 156. 158. 160. 162. 164. 166.
 168. 170. 172. 174. 176. 178. 180. 182. 184. 186. 188. 190. 192. 194.
 196. 198.]

カーネルバージョン2を作成

@cuda.jit('void(float32[:], float32[:], float32[:])')
def cu_add2(a, b, c):
    """This kernel function will be executed by a thread."""
    i  = cuda.grid(1)

    if i > c.shape[0]:
        return

    c[i] = a[i] + b[i]

カーネルVer. 2を走らせる

device = cuda.get_current_device()

n = 100
a = np.arange(n, dtype=np.float32)
b = np.arange(n, dtype=np.float32)
c = np.empty_like(a)

tpb = device.WARP_SIZE
bpg = int(np.ceil(float(n)/tpb))
print ('Blocks per grid:', bpg)
print ('Threads per block', tpb)

cu_add2[bpg, tpb](a, b, c)
print (c)
Blocks per grid: 4
Threads per block 32
[  0.   2.   4.   6.   8.  10.  12.  14.  16.  18.  20.  22.  24.  26.
  28.  30.  32.  34.  36.  38.  40.  42.  44.  46.  48.  50.  52.  54.
  56.  58.  60.  62.  64.  66.  68.  70.  72.  74.  76.  78.  80.  82.
  84.  86.  88.  90.  92.  94.  96.  98. 100. 102. 104. 106. 108. 110.
 112. 114. 116. 118. 120. 122. 124. 126. 128. 130. 132. 134. 136. 138.
 140. 142. 144. 146. 148. 150. 152. 154. 156. 158. 160. 162. 164. 166.
 168. 170. 172. 174. 176. 178. 180. 182. 184. 186. 188. 190. 192. 194.
 196. 198.]

ベクトル化ベクトル加算

@vectorize(['int64(int64, int64)',
            'float32(float32, float32)',
            'float64(float64, float64)'],
           target='cuda')
def cu_add(a, b):
    return a + b
n = 100
a = np.arange(n, dtype=np.float32)
b = np.arange(n, dtype=np.float32)
c = cu_add(a, b)
print (c)
[  0.   2.   4.   6.   8.  10.  12.  14.  16.  18.  20.  22.  24.  26.
  28.  30.  32.  34.  36.  38.  40.  42.  44.  46.  48.  50.  52.  54.
  56.  58.  60.  62.  64.  66.  68.  70.  72.  74.  76.  78.  80.  82.
  84.  86.  88.  90.  92.  94.  96.  98. 100. 102. 104. 106. 108. 110.
 112. 114. 116. 118. 120. 122. 124. 126. 128. 130. 132. 134. 136. 138.
 140. 142. 144. 146. 148. 150. 152. 154. 156. 158. 160. 162. 164. 166.
 168. 170. 172. 174. 176. 178. 180. 182. 184. 186. 188. 190. 192. 194.
 196. 198.]

2D(2次元)バージョン

@cuda.jit('void(float32[:,:], float32[:,:], float32[:,:])')
def cu_add_2d(a, b, c):
    """This kernel function will be executed by a thread."""
    i, j  = cuda.grid(2)

    if (i < c.shape[0]) and (j < c.shape[1]):
        c[i, j] = a[i, j] + b[i, j]
    cuda.syncthreads()

Low level cuda.jitは、以下の例のように、ブロック毎のスレッド数とグリッド毎のブロック数によるカーネルの正確なインスタンス化が要求される。

device = cuda.get_current_device()

n = 480
p = 320
a = np.random.random((n, p)).astype(np.float32)
b = np.ones((n, p)).astype(np.float32)
c = np.empty_like(a)

threadsperblock = (16, 16)
blockspergrid_x = (n + threadsperblock[0]) // threadsperblock[0]
blockspergrid_y = (p + threadsperblock[1]) // threadsperblock[1]
blockspergrid = (blockspergrid_x, blockspergrid_y)

print (blockspergrid, threadsperblock)

cu_add_2d[blockspergrid, threadsperblock](a, b, c)
print (a[-5:, -5:])
print (b[-5:, -5:])
print (c[-5:, -5:])
(31, 21) (16, 16)
[[0.3248 0.9058 0.2799 0.242  0.4865]
 [0.0489 0.8854 0.3703 0.037  0.6805]
 [0.067  0.8357 0.9608 0.5922 0.5895]
 [0.3291 0.6385 0.1935 0.7923 0.1367]
 [0.5513 0.4608 0.1582 0.9282 0.6274]]
[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]
[[1.3248 1.9058 1.2799 1.242  1.4865]
 [1.0489 1.8854 1.3703 1.037  1.6805]
 [1.067  1.8357 1.9608 1.5922 1.5895]
 [1.3291 1.6385 1.1935 1.7923 1.1367]
 [1.5513 1.4608 1.1582 1.9282 1.6274]]

cudaプログラミングは後々の機械学習にかなり役立つらしいので、機械学習と並行して学習するのがいいのかもしれない。

参考サイトMassively parallel programming with GPUs