Numpy(np.exp)とPyCUDAで指数関数計算の速度比較

今日はこのサイトを参考にして、numpy.exp, cumath.exp, ElementwiseKernel, SourceModuleの速度比較をする。numpy.expというのは、”Calculate the exponential of all elements in the input array”「入力配列の全要素のネイピア数のべき乗(指数関数)を算出する。」関数のことらしい。

スポンサーリンク

速度比較用コードを作成

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray, cumath
from pycuda.elementwise import ElementwiseKernel
from pycuda.compiler import SourceModule
import numpy as np

start = drv.Event()
end = drv.Event()

kernel = ElementwiseKernel(
   "double *a,double *b",
   "b[i] = exp(a[i]);")

mod = SourceModule("""
  __global__ void gexp(double *a,double *b,int n)
  {
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    while (i < n) {
    b[i] = exp(a[i]);
    i += blockDim.x * gridDim.x;
    }
  }
  """)
knl = mod.get_function("gexp")

results = []
for N in [10**4, 10**5, 10**6, 10**7, 10**8]:
    a = 2*np.ones(N,dtype=np.float64)
    a_gpu = gpuarray.to_gpu(a)
    b_gpu = gpuarray.zeros_like(a_gpu)
    c_gpu = gpuarray.to_gpu(a)
    d_gpu = gpuarray.zeros_like(c_gpu)
    
    start.record()
    np.exp(a)
    end.record()
    end.synchronize()
    sec1 = start.time_till(end)*1e-3
    print ("Numpy",sec1)
    
    start.record() # start timing
    kernel(a_gpu,b_gpu)
    end.record() # end timing
    end.synchronize()
    sec2 = start.time_till(end)*1e-3
    print ("Kernel",sec2)
    print (np.allclose(np.exp(a),b_gpu.get()))
    
    start.record() # start timing
    knl(c_gpu,d_gpu,block=(1024,1,1),grid=(N//1024,1,1))
    end.record() # end timing
    end.synchronize()
    sec3 = start.time_till(end)*1e-3
    print ("knl",sec3)
    print (np.allclose(np.exp(a),d_gpu.get()))
    
    start.record()
    cumath.exp(a_gpu)
    end.record()
    end.synchronize()
    sec4 = start.time_till(end)*1e-3
    print ("Cumath", sec4)
    #print (np.allclose(np.exp(a),cumath.exp(a_gpu))
    results.append([N,sec1,sec2,sec3,sec4])
Numpy 2.0479999948292972e-06
Kernel 0.058564609527587894
True
knl 5.9519998729228974e-05
True
Cumath 0.00011289600282907486
Numpy 0.0014398399591445924
Kernel 6.447999924421311e-05
True
knl 7.606399804353714e-05
True
Cumath 0.00441871976852417
Numpy 0.014872480392456054
Kernel 0.000333407998085022
True
knl 0.0003580479919910431
True
Cumath 0.0006108800172805786
Numpy 0.15476223754882812
Kernel 0.002812959909439087
True
knl 0.0030813119411468508
True
Cumath 0.003414207935333252
Numpy 1.515461669921875
Kernel 0.02290995216369629
True
knl 0.024895999908447267
True
Cumath 0.025443199157714844

結果をグラフ化するコードを作成

import matplotlib.pyplot as plt
results = np.array(results)
legends = []
nH = results[:5, 0:1]
rows = results[:5,1:6]
plt.semilogx(nH,rows, 'o-')
legends += ['' + s for s in ['numpy','wise','source','cumath']]
plt.rcParams['figure.figsize'] = 18, 10
plt.rcParams["font.size"] = "20"
plt.ylabel('Seconds')
plt.xlabel('Value of N')
plt.legend(legends);

cumath.exp, ElementwiseKernelバージョン, SourceModuleバージョンの速度は、N=10000を除いて(ElementwiseKernelバージョンがメチャクチャ遅い)ほとんど変わらなかった。