GPU编程一瞥

之前因为SLAM中计算描述子的缘故，想到通过GPU加速编程来提高SIFT描述子的计算速度，从而达到实时的效果。于是前段时间就了解了一下GPU编程。本来该篇博客在更早的时候就应该写的，但是由于自己的拖延症，一直挨到了今天，心想，实在不能再拖下去了。本系列基本以<a href="https://developer.download.nvidia.com/books/cuda-by-example/cuda-by-example-sample.pdf">该书</a>为基础，对书中的部分代码作了大幅度改动，修正了其中的一些运行结果，并引入了OpenCV和cmake，使得gpu编程更加工程化。系列中的代码都可以在我的github主页中找到。

CUDA的安装

CUDA的安装应该说不算有难度，其安装方式在<a href="https://developer.nvidia.com/cuda-downloads">官网</a>已经给出，并按照说明安装好就OK了。我使用的是Ubuntu16.04、cuda8.0 （为什么不用Ubuntu18.04？为什么不要cuda9.0？因为跑得太快的话，坑太多，我就是从18.04回滚回来的），唯一需<span style="font-size: 1rem;">要说明的是，安装好后的二进制文件不在系统的默认搜索文件中，我们需要export PATH=$PATH:/usr/local/cuda-8.0/bin 来引入二进制命令文件。这句话也可以加入.zshrc中。打开终端测试nvcc命令， 看其是否起作用，这个命令是我们用来编译.cu文件的命令，地位等同于gcc和g++。</span>

初识CUDA

gpu编程与普通c族语言编程最大最核心的区别就是： cuda会把应该在gpu上运行的程序“拷贝”若干份，每个gpu的小单元都会得到一份“程序”，然后执行。 这样，对于cpu来说是一套耗时的重复简单操作，可以gpu上得到大规模的并行化执行。现在我们就来完成一个极简单的gpu小程序。

#include<iostream>
using namespace std;

__global__ void kernel(void){
}

int main(){
    kernel<<<1,1>>>();
    cout<<"Hello world"<<endl;
    return 0;
}

不要吐槽我的大括号风格 ，我两种风格都喜欢。该程序再gpu上什么也没做，然后打印hello world， 然后结束。需要说明的是，

1. global_ 是用来指示该程序是在gpu上运行的主程序入口。
2. kernel<<<1,1>>>() 是指示该程序再gpu上的部署方式，第一个数字，指明改程序分配多少个block，第二程序指明每个block中的thread个数（block和thread在以后会说明）。 三个尖括号大约是为了区别c族语言中的一些符号。
   然后我们将其命名为 hello.cu 。然后用nvcc -arch=sm_30 hello.cu 来编译它，这里的-arch=sm_30是指明价格，抑制警告。最后我们得到可执行二进制文件，执行之。我们以后的程序大都类似于这个架构，一个gpu函数，一个cpu上的main函数。

稍微复杂一点的程序

按照书，接下来的程序应该是获取gpu的一些属性，然后在gpu上完成一个变量的加减，但是我们跳过以上步骤，接下来完成一个向量的加法。

using namespace std;
#define N 100

__global__ void add(int *a, int *b, int *c){
    int tid = blockIdx.x;
    if(tid<N){
        c[tid] = a[tid] +b[tid];
    }
}
int main(){
    int a[N], b[N], c[N];
    int *dev_a, *dev_b, *dev_c;

    for(int i=0; i<N; i++){
        a[i] = i;
        b[i] = 2;
    }
    cudaMalloc((void**) &dev_a, N*sizeof(int));
    cudaMalloc((void**) &dev_b, N*sizeof(int));
    cudaMalloc((void**) &dev_c, N*sizeof(int));

    cudaMemcpy(dev_a, a, N*sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N*sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_c, c, N*sizeof(int), cudaMemcpyHostToDevice);

    add<<<N,1>>>(dev_a, dev_b, dev_c);

    cudaMemcpy(c, dev_c, N*sizeof(int), cudaMemcpyDeviceToHost);

    for(auto item:c){
        cout<<item<<endl;
    }
    cudaFree(dev_a);    
    cudaFree(dev_b);    
    cudaFree(dev_c);    

    return 0;
}

1. global的意义前文已经说明，关键是函数中的tid变量。正如前文所说，这里的kernel函数并不是只有一份，它可能会被复制很多很多份，然后每份由单独一个一个gpu单元去执行。向量的加法中，我们需要每个gpu小单元去计算一个加法运算。譬如，第一个单元算a[0]+b[0]，第二个单元算a[1]+a[2]……那么我们需要知当前处于哪个单元，然后让它计算相应的加法。blockIdx.x就是每个小单元的ID，基于此，我们就能让每个单元做相应的加法，而不混乱。
2. 当运行再gpu上的程序是没法直接使用内存中的数据的，他只能访问显存，因此，我们需要吧要计算的加数复制到显存中去。那么首先我们需要显存分配空间，这就是cudaMalloc的作用，第一个是分配的地址，第二个分配的大小。
3. 接下来就是吧内存中的数据复制进显存，于是我们使用cudaMemcpy()，第一个参数是目标位置，第二是源位置，第三个是复制的大小，第四个参数是复制方向。y因为计算的结果dev_c也是显存中的地址，所以后面我们又将其复制到内存中去。
4. 如上文所说，三个尖括号是来指示部署方式，本程序是使用了N的block，每个block一个thread。结合第一条。……第N的block刚好计算a[N]+b[N]。
5. 最后，我们释放显存指针。如果不释放的话，下次执行的时候，可能受到上次执行分配的内存的影响。
   程序解释完了，我们引入cmake来对.cu程序进行编译。cmake编译的好处自不多言。

cmake_minimum_required(VERSION 2.8)
project(chapter01)

set(CMAKE_CXX_FLAGS "-std=c++11 -Wall")

find_package(CUDA)
CUDA_ADD_EXECUTABLE(add vector_add.cu)

当然，这样需要你的 /usr/share/cmake-3.5/Modules有FindCUDA.cmake时，它才会work。理论上是，安装的时候，该模块就嵌入了cmake的默认cmake的相应路径中了。以后我们使用cmake的时候，很多的库的模块并没有嵌入，其中最简单的方法就是，吧相应的Findxxx.cmake复制到上述的路径中就可以了。

接下来自然是cmake 和make，得到我们的执行程序，可以试一试看看结果是否正确。

到此，你已经掌握了gpu编程的基本方法。

第一个GPU编程的实例

分形几何，分形几何是一个很有意思的数学分支。今天我们就来简单的画一个分形图案--julia集。在复数域，如果对于一个初始点，我们对他使用公式对它进行迭代，如果迭代值不发散，那么我就说处再julia集中。

思路是，我们给出一副图，我们对它的每一个像素坐标进行考察，因为像素左边只能取整，我们乘以一个尺度，将其缩放，然后对其进行1000轮的验证，检验是否发散，如果不发散，该点我们就认为其在Julia中。程序如下：

#include<opencv2/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<memory>
#include "cuda.h"
#include<iostream>
#include<memory>

using namespace std;
#define DIM 1000

class cuComplex{
public:
    double r,i;

    __device__ cuComplex(double a, double b): r(a),i(b){}
    __device__ cuComplex operator* (const cuComplex& a){
        return cuComplex(r*a.r - i*a.i, r*a.i + i*a.r);
    }
    __device__ cuComplex operator+ (const cuComplex& a){
        return cuComplex(r + a.r, i + a.i);
    }
    __device__ double magnitude2(void){
        return r*r+ i*i;
    }
};

__device__ int julia(int x, int y){
    double cx = 1.5*(double)(x-DIM/2)/(DIM/2);
    double cy = 1.5*(double)(y-DIM/2)/(DIM/2);

    cuComplex a(cx,cy);
    cuComplex c(-0.8, 0.156);

    for(int i=0; i<200; i++){
        a = a*a + c;
        if(a.magnitude2() > 1000)//二范数
            return 0;
    }
    return 1;
}

__global__ void kernel(unsigned char *ptr){
    int x = blockIdx.x;
    int y = blockIdx.y;

    int offset = y*gridDim.x + x;

    int jV = julia(x, y);
    ptr[offset*4 + 0] = 0;
    ptr[offset*4 + 1] = 0;
    ptr[offset*4 + 2] = 255*jV;
    ptr[offset*4 + 3] = 255;
}

int main(){
    cv::Mat_<cv::Vec3b> img(DIM, DIM);
    unsigned char ptr[4*DIM*DIM];
    unsigned char *dev_img;
    cudaMalloc((void**) &dev_img, 4*DIM*DIM*sizeof(unsigned char));

    dim3 grid(DIM,DIM);
    kernel<<<grid, 1>>>(dev_img);

    cudaMemcpy(
        ptr, 
        dev_img, 
        4*DIM*DIM*sizeof(unsigned char), 
        cudaMemcpyDeviceToHost);
    for(int i=0; i<img.rows; i++){
        for(int j=0; j<img.cols; j++){
            for(int ch=0; ch<4; ch++){
                img.at<cv::Vec3b>(i,j)[ch] = ptr[ 4*(j*DIM+i) + ch];
            }
        }
    }
    cv::imshow("julia", img);
    cv::waitKey(0);

    return 0;
}

程序说明：

1. device用以指示该函数运行在gpu上，我们知道gpu上的代码不能直接访问内存，当然也就不能调用普通的c簇函数的了，因此用device来指示。
2. 我们先声明了一个复数类，并在其中定义了复数的基本算法，而且指明这些函数用device修饰，因为我们的迭代验证是再GPU上进行，当然复数需要能够被gpu调用。
3. 我们这里使用opencv来显示图像。
   cmake的内容如下：

cmake_minimum_required(VERSION 2.8)
project(chapter01)

set(CMAKE_CXX_FLAGS "-std=c++11 -Wall")
find_package(OpenCV REQUIRED) 
include_directories(${OpenCV_INCLUDE_DIRS})
find_package(CUDA)

CUDA_ADD_EXECUTABLE(add vector_add.cu)
CUDA_ADD_EXECUTABLE(g julia_gpu.cu)
target_link_libraries(g ${OpenCV_LIBS})

执行出来的结果也是相当的漂亮的：

julia.png

未完待续...