GPU编程再瞥

在上篇文章中,我们基本了解了gpu编程的原理,并简单的绘制了一个分形几何。但是前文的编程都是基于多block单thread的编程,接下我们要去了解多线程编程。也就是每个block中我们有多个thread。

GPU在程序中的排列

诚如上篇博客中提到的,GPU的每个小单元是block,每个block又可以有多个thread来执行操作,这样,整个gpu的布置就划分为block和thread两个层面。同样的,我们需要对“每一套被拷贝的程序”进行定位,以确定他们应该确实执行的动作。这里的定位我们可以借鉴下面的一张图

整个大图就是一个任务,我们叫grid。其中的一个大格就是一个block,而每个block中的小格就是一个thread。每个block在grid中的定位用blockIdx.x 和blockIdx.y来确定。在每个block内部,hread在block中的定位用threadIdx.x和threadIdx.y来确定。另外,gridDim.x和gridDim.y指明 整个任务grid上对应的方向上的包含的block数目。blockIdm.x和blockDim.y意义类似,指明包含的thread的数目。所以我们假设要整个任务是一幅图片,我们每个位置的定位就成了:

    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

我们的图片再处理上当作二维来处理,事实上他在内存或显存中是线性排布的,这里的offset就是把二维坐标转化成线性,来操作相应的内存或显存的位置。

当然,如果我们的程序是一维的的操作就能完成,譬如一维向量的运算,我们就只看上图中的第一行就可以了。

多块多线程实例

以下同样是完成一个向量加法的程序,但是,我们使用的是一维的多线程。因为当向量很长的时候,gpu不能够分配到足够的block供我们使用,这个时候就可以使用多thread。以下一共用了1个block,1000个thread。当然,你也可以尝试改为多block多thread的形式。

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

using namespace std;
#define N 1000

__global__ void kernel(int *a, int *b, int *c){
    int tid = threadIdx.x;
    if(tid<N){
        c[tid] = a[tid] + b[tid];
        tid += blockDim.x * gridDim.x;
    }
}

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] = 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);

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

    cudaMemcpy(c, dev_c, N*sizeof(int), cudaMemcpyDeviceToHost);
    for(auto i:c){
        cout<<i<<endl;
    }
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);
    return 0;
}

程序不做过多的解释,唯二要说明的是,

  1. 在__global__函数中,我们应该对计算出的位置(上例中的tid)应该加以判断,保证它不越界。事实上,当我们使用多block多thread分配的时候,因为需要处理的数据的长度可能不能被合理的整除,这时候我们分配的总资源就会多于数据长度,如果不验证,就会越界。
  2. tid += blockDim.x * gridDim.x; 可以理解为,一个周期完成了,执行下一个周期,也就是执行第二行。因为有的数据长度实在太长了,我们一行执行不玩,就分成多行,也就是过个周期。同样的,对与二维数据,我们也可以用
    id += blockDim.x * gridDim.x * blockDim.y *gridDim.y.
    结合上图来理解这两句话,就容易得多。

二维数据的实例

最常见的二维数据当然就是图像了。没错,我们接下来做一个涟漪的程序,为了显示效果,我们将计算的数据从显存中读入内存,然后再调用opencv来显示。

#include<opencv2/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<memory>
#include "cuda.h"
#include <stdio.h>
#include<iostream>
#define DIM 1024
#define PI 3.1415926535897932f

__global__ void kernel( unsigned char *ptr, int ticks ) {
    // map from threadIdx/BlockIdx to pixel position
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    // now calculate the value at that position
    float fx = x - DIM/2;
    float fy = y - DIM/2;
    float d = sqrtf( fx * fx + fy * fy );
    //unsigned char grey = (unsigned char)(x);
    unsigned char grey = (unsigned char)(128.0f + 127.0f *
                                         cos(d/10.0f - ticks/7.0f) /
                                         (d/10.0f + 1.0f));    
    ptr[offset*4 + 0] = grey;
    ptr[offset*4 + 1] = grey;
    ptr[offset*4 + 2] = grey;
    ptr[offset*4 + 3] = 255;
}

int main( void ) {
    cv::Mat_<cv::Vec3b> img(DIM, DIM);
    unsigned char ptrs[4*DIM*DIM];
    unsigned char* dev_bitmap;
    cudaMalloc( (void**)&dev_bitmap, 4*DIM*DIM* sizeof(unsigned char) ) ;
    dim3 blocks(DIM/16,DIM/16);
    dim3 threads(16, 16);

    clock_t begin_ = clock();
    for(int time=0; time<100; time++)
    {
        kernel<<<blocks, threads>>>(dev_bitmap, time);
        cudaMemcpy(ptrs, dev_bitmap, 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<3; ch++)
                    img.at<cv::Vec3b>(i,j)[ch]=ptrs[ 4*(j*DIM+i) + ch];
            }
        }
        cv::imshow("test", img);
        cv::waitKey(1);
    }
    cudaFree(dev_bitmap);
    clock_t end_ = clock();
    std::cout<<"elapsed: "<<end_ - begin_<<std::endl;
    return 0;
}

cmake文件:

cmake_minimum_required(VERSION 2.8)
project(cuda_test)

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

find_package(GLUT)
include_directories(${GLUT_INCLUDE_DIR})

find_package(GLU)
include_directories(${GLU_INCLUDE_PATH})

find_package(CUDA)

CUDA_ADD_EXECUTABLE(g ripple.1.cu)
target_link_libraries(g ${OpenCV_LIBS} ${GLUT_LIBRARY} ${GLU_LIBRARY})

add_executable(c ripple.2.cpp)
target_link_libraries(c ${OpenCV_LIBS})

 

这里我们先计算得到图片坐标x y和显存偏移offset,然后转化成中心坐标fx和fy,然后产生一个余弦波纹,使该波纹成为时间的函数。然后再将图像从显存拷贝至内存中,然后用opencv显示。结果:

经过简单的改造,我们可以将上面的程序改成cpu版本,两者比较,gpu版本块了大约3倍左右。为什么这么少?因为数据在内存和显存之间互相拷贝是很耗时的,如果我们不将数据拷贝至内存,gpu能比cpu版本快200多倍。

书上的显示方式是采用了GL和GLUT等库,所以我当时调试完成后,再cmake中留下了这两个库的寻找和链接,这里并未作移除,如果读者的电脑上找不到这两个库,将其移除便是。

GPU的共享内存及同步

上一小结指出数据再内存和显存之间拷贝是耗时的,因为两者间的带宽有限。所以,接下来要介绍的就是gpu的共享内存shared memory。

共享内存是每个block内的片上cache,速度远快于片外的DRAM。每个block的共享内存对与该block内的thread都是可见的,对于其他block的thread是不可见的。因为本block的所有thread都是可见的,所以,在共享内存中,可以实现线程通信,但是线程通信带来的问题就是同步问题。

通常我们用__shared__来指示一段空间应该分配为共享内存。举个简单的例子,对于向量的內积我们应该很熟悉公式了吧。再gpu的计算中,我们可以把每个乘积项的结果放在显存中,然后再读入,然后在求和。但是更好的办法是,把乘积项的结果放在block的共享内存中,然后每个block对自己处里的部分进行求和,得到一个小和式的结果,最终我们把各个小和式按找二叉树的原则相加,得到整个和式的结果,也就是内积。多说无益,从程序中最能体会其中的味道:

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

using namespace std;
#define N 2560
#define threadsPerBlock 256 
#define blockPerGrid 10

__global__ void kernel(int *a, int *b, int *c){
    __shared__ int cache[threadsPerBlock];  //声明共享内存
    int cacheIndex = threadIdx.x;

    int tid =  threadIdx.x + blockIdx.x * gridDim.x;
    int tmp =0;
    while(tid<N){
        tmp += a[tid] * b[tid];
        tid += blockDim.x * gridDim.x;
    }
    cache[cacheIndex] = tmp;
    __syncthreads();

    int i = threadsPerBlock/2;
    while(i !=0 ){
        if(cacheIndex < i){
            cache[cacheIndex] += cache[cacheIndex+i];
            cache[cacheIndex+i]=0;
        }
        __syncthreads();
        i /=2;
    }
    if(cacheIndex ==0 ){
        c[blockIdx.x] = cache[0];
    }
}

int main(){
    int a[N], b[N], c[blockPerGrid];
    int *dev_a, *dev_b, *dev_c;
    cudaMalloc((void**) &dev_a, N*sizeof(int));
    cudaMalloc((void**) &dev_b, N*sizeof(int));
    cudaMalloc((void**) &dev_c, blockPerGrid*sizeof(int));

    for(int i=0; i<N; i++){
        a[i] = i%10;
        b[i] = i%10;
    }
    cudaMemcpy(dev_a, a, N*sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N*sizeof(int), cudaMemcpyHostToDevice);
    kernel<<<threadsPerBlock, blockPerGrid>>>(dev_a, dev_b, dev_c);

    cudaMemcpy(c, dev_c, blockPerGrid*sizeof(int), cudaMemcpyDeviceToHost);
    static int sum = 0;
    for(auto it:c){
        sum +=it;
    }
    cout<<"sum is: "<<sum<<endl;
    return 0;
}
  1. 为什么每个子thread计算的结果没有直接存入cache变量中而是通过tmp中间变量然后加进cache变量中呢?正如前文所说,有时候的数据长度过长,我们不能期待一次完成,我们要分多批次来计算。所以相应的cache变量中已经有了上次和式的计算结果,我们这次需要加上去。
  2. 因为一个block内的所有thread都可以访问共享内存,为了防止意外,我们需要对所有的线程进行同步,同步使用了__syncthreads()函数,程序应当在此处等待,直到所有的thread完成动作后,才继续下一步。
  3. 然后我用二叉树的形式来对共享内存求和。需要说明的是,该步骤的求和是在一个计算周期内进行的,也就是说每个周期结束的时候,这有cache[0]的数据是有效的,而其他的和式已经被累加到cache[0]中了,应该将他们清0,防止在下一个周期被反复累加。

二维正弦图像

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

using namespace std;
#define DIM 1024
#define PI 3.14159f

__global__ void kernel(unsigned char *p){
    __shared__ float shared[16][16];

    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y* blockDim.x * gridDim.x;

    const float period = 128.0f;
    shared[threadIdx.x][threadIdx.y] = 
        255 * (sinf(x*2.0f*PI/period)+1.0f) *
              (sinf(y*2.0f*PI/period)+1.0f)/4.0f;
    __syncthreads(); 

    p[offset*4 + 0] = 0;
    //这里为什么要用15去减去x和y呢?因为在交换过程中,后面的值未必成功的计算了,这样就可以
    //凸显出__syncthread()函数的意义和效果。如果采用下面注释的那段代码的话,即使去掉前文
    //的syncthread(),也不会有噪声,因为一切都是按顺序计算的。至于为什么从二位空间的正弦图
    //变成了块状图,那是因为在 1/(64X64) 的每个小块内发生了顺序颠倒,但总体上,应该仍然保
    //持正弦趋势。
    p[offset*4 + 1] = (unsigned char)( shared[15-threadIdx.x][15-threadIdx.y]);
    //p[offset*4 + 1] = (unsigned char)( shared[threadIdx.x][threadIdx.y]);
    p[offset*4 + 2] = 0;
    p[offset*4 + 3] = 255;
}

int main(){
    cv::Mat_<cv::Vec3b> img(DIM,DIM);
    unsigned char tmp[4*DIM*DIM];
    unsigned char *dev_bitmap;

    cudaMalloc((void**)&dev_bitmap, 4*DIM*DIM*sizeof(unsigned char));

    dim3 grids(DIM/16, DIM/16);
    dim3 threads(16, 16);
    kernel<<<grids, threads>>>(dev_bitmap);

    cudaMemcpy(tmp, dev_bitmap, 4*DIM*DIM*sizeof(unsigned char), cudaMemcpyDeviceToHost);

    for(int x=0; x<img.rows; x++){
        for(int y=0; y<img.cols; y++){
            for(int ch=0; ch<3; ch++){
                img.at<cv::Vec3b>(x,y)[ch] = tmp[ 4*(y + x*DIM) + ch];
            }
        }
    } 
    cv::imshow("show", img);
    cv::waitKey(0);
    return 0;
}

这是截图的一部分,为什么呈现出方块状,程序中的注释已经说明得很清楚了。如果我们采用注释掉的那行代码,我们就可以看到一副二维正弦图。这里采用交换的意义也在程序中说明了。读者可以尝试 注释掉同步的代码,看看会有什么样的图像。

发表评论

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