在上篇文章中，我们基本了解了gpu编程的原理，并简单的绘制了一个分形几何。但是前文的编程都是基于多block单thread的编程，接下我们要去了解多线程编程。也就是每个block中我们有多个thread。
GPU在程序中的排列
诚如上篇博客中提到的，GPU的每个小单元是block，每个block又可以有多个thread来执行操作，这样，整个gpu的布置就划分为block和thread两个层面。同样的，我们需要对“每一套被拷贝的程序”进行定位，以确定他们应该确实执行的动作。这里的定位我们可以借鉴下面的一张图
block_thread.png
整个大图就是一个任务，我们叫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;
}
程序不做过多的解释，唯二要说明的是，
- 在global函数中，我们应该对计算出的位置（上例中的tid）应该加以判断，保证它不越界。事实上，当我们使用多block多thread分配的时候，因为需要处理的数据的长度可能不能被合理的整除，这时候我们分配的总资源就会多于数据长度，如果不验证，就会越界。
- 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显示。结果：
ripple.png
经过简单的改造，我们可以将上面的程序改成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;
}
- 为什么每个子thread计算的结果没有直接存入cache变量中而是通过tmp中间变量然后加进cache变量中呢？正如前文所说，有时候的数据长度过长，我们不能期待一次完成，我们要分多批次来计算。所以相应的cache变量中已经有了上次和式的计算结果，我们这次需要加上去。
- 因为一个block内的所有thread都可以访问共享内存，为了防止意外，我们需要对所有的线程进行同步，同步使用了__syncthreads()函数，程序应当在此处等待，直到所有的thread完成动作后，才继续下一步.
- 然后我用二叉树的形式来对共享内存求和。需要说明的是，该步骤的求和是在一个计算周期内进行的，也就是说每个周期结束的时候，这有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;
}
sin_sqare_img.png
这是截图的一部分，为什么呈现出方块状，程序中的注释已经说明得很清楚了。如果我们采用注释掉的那行代码，我们就可以看到一副二维正弦图。这里采用交换的意义也在程序中说明了。读者可以尝试 注释掉同步的代码，看看会有什么样的图像。
。所以我们假设要整个任务是一幅图片，我们每个位置的定位就成了：
网友评论