GPU编程一瞥

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

CUDA的安装

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

初识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上完成一个变量的加减,但是我们跳过以上步骤,接下来完成一个向量的加法。

#include<iostream>
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集。在复数域,如果对于一个初始点\(Z_0\),我们对他使用公式$$Z_{n+1}=Z^2_n+C$$对它进行迭代,如果迭代值不发散,那么我就说\(Z_0\)处再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})

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

未完待续…

发表评论

电子邮件地址不会被公开。 必填项已用*标注