CUDA&GPU in Python

CPU

既然讲到CPU就来复习一下计算机组成原理的一点关于CPU的知识吧~冯诺依曼机有五大组成部分，包括控制器，运算器，存储器，输入设备，输出设备，早期的冯诺依曼机是以运算器为核心的，现代计算机依然沿用冯诺依曼体系，只是不再以运算器为核心，而是以存储器为核心了。
现代计算机的系统包括两大部分，一是硬件系统，二是软件系统，详情参考下面我做的思维导图

那么CPU是什么作用呢？它的功能主要是解释计算机指令以及处理计算机软件中的数据。是一台计算机的运算核心和控制核心。

从cpu的结构中可以看出，cpu需要大量的空间去放置存储单元和控制单元，相比之下计算单元只占据了很小的一部分，所以它在大规模并行计算能力上极受限制，而更擅长于逻辑控制。因为遵循冯诺依曼架构（存储程序，顺序执行），CPU总是按照指令一步一步进行计算，随着人们对更大规模与更快处理速度的需求的增加，渐渐变得有些力不从心。

GPU

并行计算(Parallel Computing)是指同时使用多种计算资源解决计算问题的过程，是提高计算机系统计算速度和处理能力的一种有效手段。它的基本思想是用多个处理器来共同求解同一问题，即将被求解的问题分解成若干个部分，各部分均由一个独立的处理机来并行计算。
GPU的构成相对简单，有数量众多的计算单元和超长的流水线，特别适合处理大量的类型统一的数据。但GPU无法单独工作，必须由CPU进行控制调用才能工作。CPU可单独作用，处理复杂的逻辑运算和不同的数据类型，但当需要大量的处理类型统一的数据时，则可调用GPU进行并行计算。

CPU VS GPU

CUDA in Python

编程工具：numba

import matplotlib.pyplot as plt
import numpy as np
import math
import time
from numba import cuda

winter = plt.imread('./winter.jpg')
print(winter.shape)

plt.imshow(winter)
plt.show()

原图展示：

不用GPU,尝试CPU进行伽马变换：

#伽马变换
def cpu_process(src,dst,gamma=2.5):
    h,w,c = src.shape
    for h_id in range(h):
        for w_id in range(w):
            for c_id in range(c):
                color = (src[h_id,w_id,c_id]/255)**gamma*255
                dst[h_id,w_id,c_id] = np.uint8(color)

cpu_img = np.zeros(winter.shape,dtype=np.uint8)
tick = time.time()
cpu_process(src=winter,dst=cpu_img)
tock = time.time()
print(f'CPU process: {tock-tick}s')
plt.imshow(cpu_img)
plt.show()

变换结果：

耗时很长：

#GPU-NUMBA
empty_img = np.zeros(winter.shape,dtype=np.uint8)#定义一个空图，此时在CPU上
gpu_img = cuda.to_device(empty_img)              #把空图拷贝到显存GPU中

gpu_src_img = cuda.to_device(winter)
h,w,c = winter.shape
#指定并行处理的线程数，需要指定block和grid，
threads_per_block = (16,16)                                  #每一个block的线程数，此处规定二维的线程数
blocks_per_grid_x = int(math.ceil(h/threads_per_block[0]))   #每一个block的x维度，图高除第一个维度
blocks_per_grid_y = int(math.ceil(w/threads_per_block[1]))   #每一个block的y维度，图高除第二个维度
blocks_per_grid = (blocks_per_grid_x,blocks_per_grid_y)      #x维度和y维度共同构成blocks，x135 block，y240
@cuda.jit #装饰器
def gpu_process(src,dst):
    threads_idx_x = cuda.blockIdx.x*cuda.blockDim.x+cuda.threadIdx.x
    threads_idx_y = cuda.blockIdx.y*cuda.blockDim.y+cuda.threadIdx.y
    for c_id in range(3):
        color = (src[threads_idx_x,threads_idx_y,c_id]/255)**2.5*255
        dst[threads_idx_x,threads_idx_y,c_id] = np.uint8(color)

tick = time.time()
gpu_process[blocks_per_grid,threads_per_block](gpu_src_img,gpu_img)
tock = time.time()
print(f'GPU process: {tock-tick}s')

用时明显缩短：

但是numba需要指定block和grid，

通过

gpu_process[blocks_per_grid,threads_per_block](gpu_src_img,gpu_img)

进行调用。

编程工具：cupy

cuda版本的numpy。当规模较小时效果可能不如numpy,适用于规模较大的情况。但两者不是可以互相完全替代的，有的三方库不认识cupy

#GPU-cupy

#安装包cupy失败，需要安装对应电脑的cuda版本的，这里是10.1 cupy-cuda101
import cupy as cp
s = time.time()
x_gpu = cp.ones((10,100,1000)) #用cupy就cp.,用numpy就np.
e = time.time()
print(e - s)

这里踩了一个小坑，不能直接pip install cupy,需要安装对应自己cuda版本的cupy，我这里是cuda10.1，所以我执行的指令是 pip install cupy-cuda101
cp:

np:

没体现出速度优势，换了一个运算写法：

import numpy as np
import cupy as cp
s = time.time()
x_gpu = cp.random.randn(100000000)**2
#x_cpu = np.random.randn(100000000)**2
e = time.time()
print(e - s)

CP:

np:这里我从100000一直试到了100000000，np的计算速度终于比cp慢了……可见确实cupy在大规模的时候才能体现出速度优势。

编程工具：pycuda

from pycuda import autoinit
from pycuda import gpuarray
from pycuda.elementwise import ElementwiseKernel
import matplotlib.pyplot as plt
import numpy as np
import time

mandelbrot_kernel = ElementwiseKernel(
"pycuda::complex<float>*lattice, float *mandelbrot_graph, int max_iters, float upper_bound",
"""
mandelbrot_graph[i] = 1;
pycuda::complex<float> c = lattice[i];
pycuda::complex<float> z(0,0);

for (int j=0;j<max_iters;j++){
z = z*z + c;
if (abs(z)>upper_bound){
mandelbrot_graph[i]=0;
break;
}
}
""",
"mandel_kernel")

def gpu_mandelbrot(width=600,height=600,real_low=-2,real_high=2,img_low=-2,img_high=2,\
                    max_iters=50,upper_bound=2):
    real_vals = np.matrix(np.linspace(real_low,real_high,width),dtype=np.complex64)
    imag_vals = np.matrix(np.linspace(img_low,img_high,height),dtype=np.complex64) * 1j
    mandelbrot_lattice = np.array(real_vals + imag_vals.transpose(), dtype=np.complex64)
    
    mandelbrot_lattice_gpu = gpuarray.to_gpu(mandelbrot_lattice)
    mandelbrot_graph_gpu = gpuarray.empty(shape=mandelbrot_lattice.shape,dtype=np.float32)
    mandelbrot_kernel(mandelbrot_lattice_gpu,mandelbrot_graph_gpu,np.int32(max_iters),np.float32(upper_bound))
    mandel_graph_cpu = mandelbrot_graph_gpu.get()
    return mandel_graph_cpu

img = gpu_mandelbrot()
plt.imshow(img)
plt.show()

tick = time.time()
img = gpu_mandelbrot()
tock = time.time()
print(f'PYCUDA mandelbrot algrithom needs {tock-tick}s.')

编程工具：pytorch

指定cuda，无需再声明todevice()

cuda=torch.device('cuda')
cuda0=torch.device('cuda:0')  #第n块显卡
cuda2=torch.device('cuda:2')

import torch
def torch_complex_plane(xrange=(-2,2),yrange=(-2,2),res=1024):
    x = torch.linspace(xrange[0],xrange[1],res).cuda()
    y = torch.linspace(yrange[0],yrange[1],res).cuda()
    real_plane, imag_plane = torch.meshgrid(x,y)
    cplane = tuple([real_plane.transpose(0,1),imag_plane.transpose(0,1)])
    return cplane


def torch_complex_magnitude(real, imag):
    return torch.sqrt(real ** 2 + imag ** 2)


def torch_mandelbrot_method(c, z, max_iters=50, upper_bound=2):
    origin_z = torch.zeros_like(z[0]).cuda()
    for i in range(max_iters):
        mask = torch.lt(torch_complex_magnitude(*z), upper_bound)
        origin_z += mask.to(torch.float32)
        z = (z[0] ** 2 - z[1] ** 2 + c[0], 2 * z[0] * z[1] + c[1])

    return origin_z.cpu()
c = torch_complex_plane()
z_init = torch_complex_plane()
img = torch_mandelbrot_method(c,z_init)
plt.figure(dpi=600)
plt.imshow(img)
plt.show()

踩了一个小坑，报错【pytorch】OMP: Error #15: Initializing libiomp5md.dll, but found libiomp5md.dll already initialized.
解决方法：在前面加了 import os
os.environ[“KMP_DUPLICATE_LIB_OK”]=“TRUE”
参考


加时间戳time.time（）记录了一下耗时：