Numba:高性能計算的高生產率

在這篇文章中，筆者將向你介紹一個來自Anaconda的Python編譯器Numba，它可以在CUDA-capable GPU或多核cpu上編譯Python代碼。Python通常不是一種編譯語言，你可能想知道為什么要使用Python編譯器。答案當然是：運行本地編譯的代碼要比運行動態的、解譯的代碼快很多倍。Numba允許你為Python函數指定類型簽名，從而在運行時啟用編譯(這就是“Just-in-Time”，即時，也可以說JIT編譯)。Numba動態編譯代碼的能力意味著你不會因此而拋棄Python的靈活性。這是向提供高生產率編程和高性能計算的完美結合邁出的一大步。

使用Numba可以編寫標準的Python函數，并在CUDA-capable GPU上運行它們。Numba是為面向數組的計算任務而設計的，很像大家常用的NumPy庫。在面向數組的計算任務中，數據并行性對于像GPU這樣的加速器是很自然的。Numba了解NumPy數組類型，并使用它們生成高效的編譯代碼，用于在GPU或多核CPU上執行。所需的編程工作可以很簡單，就像添加一個函數修飾器來指示Numba為GPU編譯一樣。例如，在下面的代碼中，@ vectorize decorator會生成一個編譯的、矢量化的標量函數在運行時添加的版本，這樣它就可以用于在GPU上并行處理數據數組。

要在CPU上編譯和運行相同的函數，我們只需將目標更改為“CPU”，它將在編譯水平上帶來性能，在CPU上向量化C代碼。這種靈活性可以幫助你生成更可重用的代碼，并允許你在沒有GPU的機器上開發。

import numpy as np

from numba import vectorize

@vectorize(['float32(float32, float32)'], target='cuda')

def Add(a, b):

return a + b

# Initialize arrays

N = 100000

A = np.ones(N, dtype=np.float32)

B = np.ones(A.shape, dtype=A.dtype)

C = np.empty_like(A, dtype=A.dtype)

# Add arrays on GPU

C = Add(A, B)

關于Python 的GPU-Accelerated庫

CUDA并行計算平臺的優勢之一是其可用的GPU加速庫的闊度。Numba團隊的另一個項目叫做pyculib，它提供了一個Python接口，用于CUDA cuBLAS(dense linear algebra，稠密線性代數)，cuFFT(Fast Fourier Transform，快速傅里葉變換)，和cuRAND(random number generation，隨機數生成)庫。許多應用程序都能夠通過使用這些庫獲得顯著的加速效果，而不需要編寫任何特定于GPU的代碼。例如，下面的代碼使用“XORWOW”偽隨機數生成器在GPU上生成100萬個均勻分布的隨機數。

import numpy as np

from pyculib import rand as curand

prng = curand.PRNG(rndtype=curand.PRNG.XORWOW)

rand = np.empty(100000)

prng.uniform(rand)

print rand[:10]

CUDA Python的高并行性

Anaconda(原名Continuum Analytics)認識到，在某些計算上實現大的速度需要一個更具表現力的編程接口，它比庫和自動循環矢量化更詳細地控制并行性。因此，Numba有另一組重要的特性，構成了其非正式名稱“CUDA Python”。Numba公開了CUDA編程模型，正如CUDA C/ C++，但是使用純python語法，這樣程序員就可以創建自定義、調優的并行內核，而不會放棄python帶來的便捷和優勢。Numba的CUDA JIT(通過decorator或函數調用可用)在運行時編譯CUDA Python函數，專門針對你所使用的類型，它的CUDA Python API提供了對數據傳輸和CUDA流的顯式控制，以及其他特性。

下面的代碼示例演示了一個簡單的Mandelbrot設置內核。請注意，mandel_kernel函數使用Numba提供的cuda.threadIdx，cuda.blockIdx，cuda.blockDim和cuda.gridDim架構來計算當前線程的全局X和Y像素索引。與其他CUDA語言一樣，我們通過插入在括號內一個“執行配置”(CUDA-speak用于線程數和線程塊啟動內核)，在函數名和參數列表之間中: mandel_kernel[griddim, blockdim](-2.0, 1.0, -1.0, 1.0, d_image, 20)。你還可以看到使用to_host和to_device API函數來從GPU中復制數據。

Mandelbrot的例子將在Github上持續更新。

@cuda.jit(device=True)

def mandel(x, y, max_iters):

"""

Given the real and imaginary parts of a complex number,

determine if it is a candidate for membership in the Mandelbrot

set given a fixed number of iterations.

"""

c = complex(x, y)

z = 0.0j

for i in range(max_iters):

z = z*z + c

if (z.real*z.real + z.imag*z.imag) >= 4:

return i

return max_iters

@cuda.jit

def mandel_kernel(min_x, max_x, min_y, max_y, image, iters):

height = image.shape[0]

width = image.shape[1]

pixel_size_x = (max_x - min_x) / width

pixel_size_y = (max_y - min_y) / height

startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x

startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y

gridX = cuda.gridDim.x * cuda.blockDim.x;

gridY = cuda.gridDim.y * cuda.blockDim.y;

for x in range(startX, width, gridX):

real = min_x + x * pixel_size_x

for y in range(startY, height, gridY):

imag = min_y + y * pixel_size_y

image[y, x] = mandel(real, imag, iters)

gimage = np.zeros((1024, 1536), dtype = np.uint8)

blockdim = (32, 8)

griddim = (32,16)

start = timer()

d_image = cuda.to_device(gimage)

mandel_kernel[griddim, blockdim](-2.0, 1.0, -1.0, 1.0, d_image, 20)

d_image.to_host()

dt = timer() - start

print "Mandelbrot created on GPU in %f s" % dt

imshow(gimage)

在一臺帶有NVIDIA Tesla P100 GPU和Intel Xeon E5-2698 v3 CPU的服務器上，這個CUDA Python Mandelbrot代碼運行的速度比純Python版本快1700倍。1700倍似乎有些不切實際，但請記住，我們正在比較編譯的、并行的、GPU加速的Python代碼來解釋CPU上的單線程Python代碼。

今天開始使用Numba吧

Numba為Python開發人員提供了一種簡單的進入GPU加速計算的方法，并提供了一種使用越來越復雜的CUDA代碼的方法，其中至少有新語法和術語。你可以從簡單的函數decorator開始實現自動編譯函數，或者使用pyculib的強大的CUDA庫。當你提高對并行編程概念的理解時，當你需要對于并行線程的有表現力且靈活的控制時，CUDA可以在不需要你第一天就完全了解的情況下使用。

Numba是一個BSD認證的開源項目，它本身嚴重依賴于LLVM編譯器的功能。Numba的GPU后端使用了基于LLVM的NVIDIA編譯器SDK。CUDA庫周圍的pyculib包裝器也是開源且經過BSD認證的。

要開始使用Numba，第一步是下載并安裝Anaconda Python發行版，這是一個“完全免費的、用于大規模數據處理、預測分析和科學計算的Python發行版”，其中包括許多流行的軟件包(Numpy、Scipy、Matplotlib、iPython等)和“conda”，這是一個強大的包管理器。一旦您安裝了Anaconda，通過鍵入conda安裝numba cudatoolkit pyculib，安裝所需的CUDA包。然后在ContinuumIO github存儲庫中查看CUDA的Numba教程。筆者建議你在Anaconda的博客上查看Numba的帖子。

Nvidia的CUDA 架構為我們提供了一種便捷的方式來直接操縱GPU 并進行編程，但是基于

C語言的CUDA實現較為復雜，開發周期較長。而python 作為一門廣泛使用的語言，具有

簡單易學、語法簡單、開發迅速等優點。作為第四種CUDA支持語言，相信python一定會

在高性能計算上有杰出的貢獻–pyCUDA。

pyCUDA特點

pyCUDA工作流程

調用的基本例子

包含內容

pyCUDA特點

CUDA完全的python實現

編碼更為靈活、迅速、自適應調節代碼

更好的魯棒性，自動管理目標生命周期和錯誤檢測

包含易用的工具包，包括基于GPU的線性代數庫、reduction和scan，添加了快速傅里葉變換包和線性代數包LAPACK

完整的幫助文檔Wiki

pyCUDA的工作流程

具體的調用流程如下：開始

編寫python程序

python程序檢查？

調用pyCUDA編譯CUDA 源碼并上傳GPU

編譯正確？

PyCUDA’s numpy進行數據讀入處理

數據讀入處理成功？

輸出GPU 加速處理結果

結束

調用基本例子

import pycuda.autoinit

import pycuda.driver as drv

import numpy

from pycuda.compiler import SourceModule

mod = SourceModule("""

__global__ void multiply_them(float *dest, float *a, float *b)

{

const int i = threadIdx.x;

dest[i] = a[i] * b[i];

}

""")

multiply_them = mod.get_function("multiply_them")

a = numpy.random.randn(400).astype(numpy.float32)

b = numpy.random.randn(400).astype(numpy.float32)

dest = numpy.zeros_like(a)

multiply_them(

drv.Out(dest), drv.In(a), drv.In(b),

block=(400,1,1), grid=(1,1))

print dest-a*b

補充內容：對于GPU 加速python還有功能包，例如處理圖像的pythonGPU加速包—— pyGPU

以及專門的GPU 加速python機器學習包—— scikitCUDA

GPU入門

現在GPU編程正變得越來越流行，由于CPU串行執行的局限性，程序如果在GPU上運行，則可真正做到多線程，并發執行，極大減少運行時間，這對于分秒必爭的科學計算，以及新興的人工智能領域都帶來了極大的便利。 目前，GPU編程以NVIDIA的CUDA平臺為主，支持四種語言C、C++、Fortran(PGI)以及Python。目前CUDA的最新版本已經達到7.5，具體配置可以看官方指導和其它教程，這里不做具體介紹。 下面我們具體來看Python的GPU編程。 我的顯卡是GeForce GT 740M,安裝CUDA7.5,使用Python2.7搭配相關庫。 首先我們要引入一些必要的包?from numbapro import cuda?是cuda包是必須導入的，否則不能使用GPU。 引入之后就可以調用cuda對象了。例如，創建一個一維網格

tx = cuda.threadIdx.x

bx = cuda.blockIdx.x

bw = cuda.blockDim.x

i = tx + bx * bw

array[i] = something(i)

也可以簡化成

i = cuda.grid(1)

array[i] = something(i)

上面兩段代碼實現的功能是一樣的。

接下來我們了解一下CUDA Stream的操作 CUDA流是對CUDA設備的命令隊列，通過特殊的流，CUDA API可以變為異步，這也意味著請求肯在隊列結束前返回。存儲器傳送指令和內核調用都可以使用CUDA流。 下面我們來看示例：

stream = cuda.stream()

devary = cuda.to_device(an_array, stream=stream)

a_cuda_kernel[griddim, blockdim, stream](devary)

cuda.copy_to_host(an_array, stream=stream)

# 在an_array中的數據可能尚未就緒

stream.synchronize()

# an_array中的數據已經就緒

另一種語法是使用Python環境

stream = cuda.stream()

with stream.auto_synchronize():

devary = cuda.to_device(an_array, stream=stream)

a_cuda_kernel[griddim, blockdim, stream](devary)

devary.copy_to_host(an_array, stream=stream)

# an_array中的數據已經就緒

接下來是關于共享內存： 為了達到最大性能，CUDA內核需要使用共享內存用于緩存數據，CUDA編譯器支持使用cuda.shared.array(shape, dtype)方法用來指定使用內核中的對象。 下面看一個例子

bpg = 50

tpb = 32

n = bpg * tpb

@jit(argtypes=[float32[:,:], float32[:,:], float32[:,:]], target='gpu')

def cu_square_matrix_mul(A, B, C):

sA = cuda.shared.array(shape=(tpb, tpb), dtype=float32)

sB = cuda.shared.array(shape=(tpb, tpb), dtype=float32)

tx = cuda.threadIdx.x

ty = cuda.threadIdx.y

bx = cuda.blockIdx.x

by = cuda.blockIdx.y

bw = cuda.blockDim.x

bh = cuda.blockDim.y

x = tx + bx * bw

y = ty + by * bh

acc = 0.

for i in range(bpg):

if x < n and y < n:

sA[ty, tx] = A[y, tx + i * tpb]

sB[ty, tx] = B[ty + i * tpb, x]

cuda.syncthreads()

if x < n and y < n:

for j in range(tpb):

acc += sA[ty, j] * sB[j, tx]

cuda.syncthreads()

if x < n and y < n:

C[y, x] = acc

以上就是python GPU編程的入門介紹。

我們以FFT 為例，

看看究竟如何利用GPU進行加速。 先看示例代碼，然后進行講解。

import sys

import numpy as np

from scipy.signal import fftconvolve

from scipy import misc, ndimage

from matplotlib import pyplot as plt

from numbapro.cudalib import cufft

from numbapro import cuda, vectorize

from timeit import default_timer as timer

@vectorize(['complex64(complex64, complex64)'], target='gpu')

#目標平臺是64位機器且擁有GPU

def vmult(a, b):

return a * b

def best_grid_size(size, tpb):

bpg = np.ceil(np.array(size, dtype=np.float) / tpb).astype(np.int).tolist()

return tuple(bpg)

def main():

# 構建過濾器

laplacian_pts = '''

-4 -1 0 -1 -4

-1 2 3 2 -1

0 3 4 3 0

-1 2 3 2 -1

-4 -1 0 -1 -4

'''.split()

laplacian = np.array(laplacian_pts, dtype=np.float32).reshape(5, 5)

# 構建圖像

try:

filename = sys.argv[1]

image = ndimage.imread(filename, flatten=True).astype(np.float32)

except IndexError:

image = misc.lena().astype(np.float32)

print("Image size: %s" % (image.shape,))

response = np.zeros_like(image)

response[:5, :5] = laplacian

# CPU

ts = timer()

cvimage_cpu = fftconvolve(image, laplacian, mode='same')

te = timer()

print('CPU: %.2fs' % (te - ts))

# GPU

threadperblock = 32, 8

blockpergrid = best_grid_size(tuple(reversed(image.shape)), threadperblock)

print('kernel config: %s x %s' % (blockpergrid, threadperblock))

# cuFFT系統，觸發器初始化.

# 對于小數據集來說，效果更明顯.

# 不應該計算這里浪費的時間

cufft.FFTPlan(shape=image.shape, itype=np.complex64, otype=np.complex64)

# 開始GPU運行計時

ts = timer()

image_complex = image.astype(np.complex64)

response_complex = response.astype(np.complex64)

d_image_complex = cuda.to_device(image_complex)

d_response_complex = cuda.to_device(response_complex)

cufft.fft_inplace(d_image_complex)

cufft.fft_inplace(d_response_complex)

vmult(d_image_complex, d_response_complex, out=d_image_complex)

cufft.ifft_inplace(d_image_complex)

cvimage_gpu = d_image_complex.copy_to_host().real / np.prod(image.shape)

te = timer()

print('GPU: %.2fs' % (te - ts))

# 繪制結果

plt.subplot(1, 2, 1)

plt.title('CPU')

plt.imshow(cvimage_cpu, cmap=plt.cm.gray)

plt.axis('off')

plt.subplot(1, 2, 2)

plt.title('GPU')

plt.imshow(cvimage_gpu, cmap=plt.cm.gray)

plt.axis('off')

plt.show()

if __name__ == '__main__':

main()

看看運行結果：

時間對比：

Image size: (512L, 512L)

CPU: 0.66s

kernel config: (16, 64) x (32, 8)

GPU: 0.09s

[Finished in 61.4s]

1.導入的包

Numpy提供常見的數學函數，包含許多有用的數學庫

Scipy是python下的數值計算庫，和Numpy一樣是科學計算不可或缺的庫

Matplotlib是用以繪制二維圖形的Python模塊

Numbapro是CUDA提供的專用庫

timeit是計時工具

2.首先程序提供了一個測試數據集并進行轉化，對于CPU部分，直接調用Scipy中的fftconvolve函數計算出結果，而GPU則主要調用了numbapro中cufft庫，具體使用參考官方文檔。

3.計算出結果后繪圖展示，主要使用matplotlib的方法，只需設置plt對象并傳入參數即可。

以上在代碼中有詳細注釋，不一樣展開。總的來說，使用GPU還是得參考官方的文檔，并結合Python編程，才能解決復雜的問題并達到良好的效果。

第三篇

相信如果你使用過Python Numpy包，一定了解NumPy(Numeric Python)提供了許多高級的數值編程工具，如：矩陣數據類型、矢量處理，以及精密的運算庫。它專為進行嚴格的數字處理而產生。多為很多大型金融公司使用，以及核心的科學計算組織如：Lawrence Livermore，NASA用其處理一些本來使用C++，Fortran或Matlab等所做的任務。 　　但是由于復雜的計算，Numpy的計算效率難免受到影響，因此我們對它進行了許多優化，用于優化的包有PyPy、Numba 與 Cython，而NumbaPro就是建立在Numba和cuda基礎上的高級優化方法。 　　下面我們一起來看。 　　使用NumbaPro,我們可以對Numpy中的方法進行優化，使Python代碼可以動態編譯為機器碼，并在運行中加載，使得GPU充分發揮多線程的優勢。針對GPU，Numbapro也可以自動完成工作，并優化GPU體系結構的代碼。另外，基于CUDA API編寫的Python代碼也可以有效地利用硬件。 　　說了這么多，下面就讓我們從簡單的示例開始學習。

from numbapro import vectorize

@vectorize(['float32(float32, float32)'], target='cpu')

def sum(a, b):

return a + b

如果需要使用GPU來運行，只需要將第二行改成@vectorize(['float32(float32, float32)'], target='gpu')

對于更復雜的操作，可以使用Just-In-Time (JIT)來編譯。

from numbapro import cuda

@cuda.jit('void(float32[:], float32[:], float32[:])')

def sum(a, b, result):

i = cuda.grid(1) # 等價于threadIdx.x + blockIdx.x * blockDim.x

result[i] = a[i] + b[i]

# 調用: sum[grid_dim, block_dim](big_input_1, big_input_2, result_array)

下面繼續看一個具體的應用：

import numpy as np

import math

import time

from numba import *

from numbapro import cuda

from blackscholes_numba import black_scholes, black_scholes_numba

#import logging; logging.getLogger().setLevel(0)

RISKFREE = 0.02

VOLATILITY = 0.30

A1 = 0.31938153

A2 = -0.356563782

A3 = 1.781477937

A4 = -1.821255978

A5 = 1.330274429

RSQRT2PI = 0.39894228040143267793994605993438

@cuda.jit(argtypes=(double,), restype=double, device=True, inline=True)

def cnd_cuda(d):

K = 1.0 / (1.0 + 0.2316419 * math.fabs(d))

ret_val = (RSQRT2PI * math.exp(-0.5 * d * d) *

(K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5))))))

if d > 0:

ret_val = 1.0 - ret_val

return ret_val

@cuda.jit(argtypes=(double[:], double[:], double[:], double[:], double[:],

double, double))

def black_scholes_cuda(callResult, putResult, S, X,

T, R, V):

# S = stockPrice

# X = optionStrike

# T = optionYears

# R = Riskfree

# V = Volatility

i = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x

if i >= S.shape[0]:

return

sqrtT = math.sqrt(T[i])

d1 = (math.log(S[i] / X[i]) + (R + 0.5 * V * V) * T[i]) / (V * sqrtT)

d2 = d1 - V * sqrtT

cndd1 = cnd_cuda(d1)

cndd2 = cnd_cuda(d2)

expRT = math.exp((-1. * R) * T[i])

callResult[i] = (S[i] * cndd1 - X[i] * expRT * cndd2)

putResult[i] = (X[i] * expRT * (1.0 - cndd2) - S[i] * (1.0 - cndd1))

def randfloat(rand_var, low, high):

return (1.0 - rand_var) * low + rand_var * high

def main (*args):

OPT_N = 4000000

iterations = 10

if len(args) >= 2:

iterations = int(args[0])

callResultNumpy = np.zeros(OPT_N)

putResultNumpy = -np.ones(OPT_N)

stockPrice = randfloat(np.random.random(OPT_N), 5.0, 30.0)

optionStrike = randfloat(np.random.random(OPT_N), 1.0, 100.0)

optionYears = randfloat(np.random.random(OPT_N), 0.25, 10.0)

callResultNumba = np.zeros(OPT_N)

putResultNumba = -np.ones(OPT_N)

callResultNumbapro = np.zeros(OPT_N)

putResultNumbapro = -np.ones(OPT_N)

time0 = time.time()

for i in range(iterations):

black_scholes(callResultNumpy, putResultNumpy, stockPrice,

optionStrike, optionYears, RISKFREE, VOLATILITY)

time1 = time.time()

print("Numpy Time: %f msec" %

((1000 * (time1 - time0)) / iterations))

time0 = time.time()

for i in range(iterations):

black_scholes_numba(callResultNumba, putResultNumba, stockPrice,

optionStrike, optionYears, RISKFREE, VOLATILITY)

time1 = time.time()

print("Numba Time: %f msec" %

((1000 * (time1 - time0)) / iterations))

time0 = time.time()

blockdim = 1024, 1

griddim = int(math.ceil(float(OPT_N)/blockdim[0])), 1

stream = cuda.stream()

d_callResult = cuda.to_device(callResultNumbapro, stream)

d_putResult = cuda.to_device(putResultNumbapro, stream)

d_stockPrice = cuda.to_device(stockPrice, stream)

d_optionStrike = cuda.to_device(optionStrike, stream)

d_optionYears = cuda.to_device(optionYears, stream)

time1 = time.time()

for i in range(iterations):

black_scholes_cuda[griddim, blockdim, stream](

d_callResult, d_putResult, d_stockPrice, d_optionStrike,

d_optionYears, RISKFREE, VOLATILITY)

d_callResult.to_host(stream)

d_putResult.to_host(stream)

stream.synchronize()

time2 = time.time()

dt = (time1 - time0) * 10 + (time2 - time1)

print("numbapro.cuda time: %f msec" % ((1000 * dt) / iterations))

delta = np.abs(callResultNumpy - callResultNumba)

L1norm = delta.sum() / np.abs(callResultNumpy).sum()

print("L1 norm: %E" % L1norm)

print("Max absolute error: %E" % delta.max())

delta = np.abs(callResultNumpy - callResultNumbapro)

L1norm = delta.sum() / np.abs(callResultNumpy).sum()

print("L1 norm (Numbapro): %E" % L1norm)

print("Max absolute error (Numbapro): %E" % delta.max())

if __name__ == "__main__":

import sys

main(*sys.argv[1:])

運行結果是：

Numpy Time: 1178.500009 msec

Numba Time: 424.500012 msec

numbapro.cuda time: 138.099957 msec

可以看出，該程序是通過引入cuda對象并使用即時編譯方法進行加速。 　　比較發現運行時間發現，運用numbapro方式加速效果明顯。

總結：

1.可以通過GPU編寫數據并行處理的程序來加快處理速度。

2.可以使用CUDA庫，例如cuRAND, cuBLAS, cuFFT 3.Python CUDA程序也可以最大化的使用硬件資源。

其他鏈接：

https://dask.org/

https://anaconda.org/talley/repo?type=conda&label=main

---------------------

作者：WeisongZhao

來源：CSDN

原文：https://blog.csdn.net/weixin_41923961/article/details/83687809

版權聲明：本文為博主原創文章，轉載請附上博文鏈接！

總結

以上是生活随笔為你收集整理的python gpu编程_Python笔记_第四篇_高阶编程_进程、线程、协程_5.GPU加速的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。