CUDA基本原理及概念
本文主要包含如下内容：
Introduction
Scalability and Portability 可扩展性和可移植性
相同的应用程序可以在一个内核上有效运行
相同的应用程序可以在更多的相同内核上有效运行
相同的应用程序可以在不同类型的内核上有效运行
相同的应用程序可以在具有不同组织和接口的系统上高效运行
GPU Programming Languages
CUDA内核由线程网格（Grid）组成；网格中的所有线程都运行相同的内核代码（单程序多个数据）；每个线程都有用于计算内存地址并进行控制决策的索引。
Block:
block可以是一维、二维、三维的，每个block可以有1到512个threads，block中所有的threads将执行相同的线程程序。Maximum number of threads per block: 512
Warps
现在的主流Streaming Multiprocessor都是基于32threads的，每一个Thread Blocks被分为32个thread Warps，即32个threads被分为1个Warps
将线程数组划分成多个块,块内的线程通过共享内存,不同块中的线程不进行交互
// Compute vector sum C = A + B
void vecAdd(float *h_A, float *h_B, float *h_C, int n)
{
int i;
for (i = 0; i<n; i++) h_C[i] = h_A[i] + h_B[i];
}
int main()
{
// Memory allocation for h_A, h_B, and h_C
// I/O to read h_A and h_B, N elements
…
vecAdd(h_A, h_B, h_C, N);
} 上面这个程序是在CPU上运行的求和程序，接下来，调用GPU编程。
cudaMalloc()：在
GPU的全局内存中分配对象，即分配一个和原始数据一样大的内存，方便读取与存储。
- 两个参数；指向分配对象的指针的地址；分配对象的大小（以字节计）
cudaFree()：从全局内存释放对象
- 一个参数：指向释放对象的指针
cudaMemcpy()：内存数据传输
- 四个参数：指向目的地；指向源的指针；复制的字节数；类型/转移方向
cudaThreadSynchronize(); 同步，等待GPU内核执行结束
__global__：GPU kernel 代码，由CPU发起，返回void，不能由其他kernel调用
__device__：由GPU kernel调用的函数，不能由其他kernel调用
__host__：在CPU上执行的函数
__shared__：变量位于共享存储器
// vecAdd CUDA Host Code
#include <cuda.h>
void vecAdd(float *h_A, float *h_B, float *h_C, int n)
{
int size = n * sizeof(float);
float *d_A, *d_B, *d_C;
// Part 1
// Allocate device memory for A, B, and C
// copy A and B to device memory
cudaMalloc((void **) &d_A, size); //在全局内存中分配对象
cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice); //内存数据传输
cudaMalloc((void **) &d_B, size);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
cudaMalloc((void **) &d_C, size);
// Part 2
// Kernel launch code – the device performs the actual vector addition
dim3 DimGrid((n-1)/256 + 1, 1, 1); //总共3维，这里使用了一维
dim3 DimBlock(256, 1, 1); //每个block有256个线程
vecAddKernel<<<DimGrid,DimBlock>>>(d_A, d_B, d_C, n); //表明内核的参数
cudaThreadSynchronize(); //同步，等待执行结束
// Part 3
// copy C from the device memory
// Free device vectors
cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost); //内存数据传输
cudaFree(d_A); cudaFree(d_B); cudaFree (d_C); //从全局内存释放对象
}
__global__ //定义核心函数
void vecAddKernel(float* A, float* B, float* C, int n) //Device code必须返回void
{
int i = threadIdx.x+blockDim.x*blockIdx.x; //查找指针真正的位置,即索引
if(i<n) C[i] = A[i] + B[i]; //条件语句防止数据溢出，最后一个block可能存在空值
}Multidimensional Grids
注意：在Device中横向是x，纵向是y，在计算行和列的时候注意变化关系
...
dim3 DimGrid((n-1)/16+1,(m-1)/16+d1,1);
dim3 DimBlock(16,16,1);
PictureKernel<<<DimGird,DimBlock>>>(d_Pin, d_Pout, m, n)
...
_global__
void PictureKernel(float* d_Pin, float* d_Pout, int height, int width)
{
// Calculate the row # of the d_Pin and d_Pout element
int Row = blockIdx.y*blockDim.y + threadIdx.y; //获取原始图像的位置
// Calculate the column # of the d_Pin and d_Pout element
int Col = blockIdx.x*blockDim.x + threadIdx.x; //获取原始图像的位置
// each thread computes one element of d_Pout if in range
if ((Row < height) && (Col < width)) {
d_Pout[Row*width+Col] = 2.0*d_Pin[Row*width+Col]; //转化为一维坐标
}
}Example:RGB to Grayscale Conversion Code
#define CHANNELS 3 // we have 3 channels corresponding to RGB
// The input image is encoded as unsigned characters [0, 255]
__global__ void colorConvert(unsigned char * grayImage, unsigned char * rgbImage, int width, int height)
{
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
if (x < width && y < height)
{
// get 1D coordinate for the grayscale image
int grayOffset = y*width + x;
// one can think of the RGB image having
// CHANNEL times columns than the gray scale image
int rgbOffset = grayOffset*CHANNELS;
unsigned char r = rgbImage[rgbOffset ]; // red value for pixel
unsigned char g = rgbImage[rgbOffset + 2]; // green value for pixel
unsigned char b = rgbImage[rgbOffset + 3]; // blue value for pixel
// perform the rescaling and store it
// We multiply by floating point constants
grayImage[grayOffset] = 0.21f*r + 0.71f*g + 0.07f*b;
}
}Example:Image Blurring
__global__
void blurKernel(unsigned char * in, unsigned char * out, int w, int h)
{
int Col = blockIdx.x * blockDim.x + threadIdx.x;
int Row = blockIdx.y * blockDim.y + threadIdx.y;
if (Col < w && Row < h)
{
int pixVal = 0;
int pixels = 0;
// Get the average of the surrounding 2xBLUR_SIZE x 2xBLUR_SIZE box
for(int blurRow = -BLUR_SIZE; blurRow < BLUR_SIZE+1; ++blurRow)
{
for(int blurCol = -BLUR_SIZE; blurCol < BLUR_SIZE+1; ++blurCol)
{
int curRow = Row + blurRow;
int curCol = Col + blurCol;
// Verify we have a valid image pixel
if(curRow > -1 && curRow < h && curCol > -1 && curCol < w)
{
pixVal += in[curRow * w + curCol];
pixels++; // Keep track of number of pixels in the accumulated total
}
}
}
// Write our new pixel value out
out[Row * w + Col] = (unsigned char)(pixVal / pixels);
}
}Registers、Shared memory、Global memory
学会合理分配CUDA内存，可以加快数据访问速度，从而加快程序运行速度。
Registers
存储在寄存器存储器中的数据只对写入它的线程可见，并且只持续该线程的生命周期。
Global memory
对应用程序内的所有threads都是可见的，持续时间为主机分配。全局内存较大，速度较慢。
Shared memory
对每一个block内的所有线程都是可见的，持续时间为block的持续时间。允许线程在彼此之间进行通信和共享数据。速度非常快，如果没有冲突，速度和Registers一样快，但是最大大小仅为48KB。适合情况：在数据需要重复读取时，若数据只用读取一次，则Shared memory并无实质性的提高。
Registers中存储的是自动变量，数据一般可以存储在shared memory中，比全局内存具有更高速度（延迟和吞吐量）的访问，访问和共享的范围：线程块
生命周期 ：线程块，相应的线程终止执行后，内容将消失。
//A Basic Matrix Multiplication
__global__ void MatrixMulKernel(float* M, float* N, float* P, int Width)
{
// Calculate the row index of the P element and M
int Row = blockIdx.y*blockDim.y+threadIdx.y;
// Calculate the column index of P and N
int Col = blockIdx.x*blockDim.x+threadIdx.x;
if ((Row < Width) && (Col < Width))
{
float Pvalue = 0;
// each thread computes one element of the block sub-matrix
for (int k = 0; k < Width; ++k) {
Pvalue += M[Row*Width+k]*N[k*Width+Col];
}
P[Row*Width+Col] = Pvalue;
}
} 以上程序访问的是Global memory,速度较慢
已知矩阵乘法是由一行乘以一列完成的，由于shared memory大小有限，需要将每个阶段的线程块的数据访问集分解为在M的一个tile和N的一个tile，其中，每个tile含有BLOCK_SIZE个元素
__syncthreads()：同步所有的操作，即块内所有的内存，用以防止RAW/WAR/WAW危害
__global__ void MatrixMulKernel(float* M, float* N, float* P, Int Width)
{
__shared__ float ds_M[TILE_WIDTH][TILE_WIDTH];
__shared__ float ds_N[TILE_WIDTH][TILE_WIDTH];
int bx = blockIdx.x; int by = blockIdx.y;
int tx = threadIdx.x; int ty = threadIdx.y;
int Row = by * blockDim.y + ty; //查找指针真正的位置,即索引
int Col = bx * blockDim.x + tx;
float Pvalue = 0;
// Loop over the M and N tiles required to compute the P element
for (int p = 0; p < n/TILE_WIDTH; ++p)
{
// Collaborative loading of M and N tiles into shared memory
ds_M[ty][tx] = M[Row*Width + p*TILE_WIDTH + tx]; //for_loop 循环加载shared_memory
ds_N[ty][tx] = N[(t*TILE_WIDTH+ty)*Width + Col];
__syncthreads(); //确保上一步操作完成
for (int i = 0; i < TILE_WIDTH; ++i)Pvalue += ds_M[ty][i] * ds_N[i][tx]; //矩阵乘法
__synchthreads(); //确保上一步操作完成
}
C[Row*Width+Col] = Cvalue;
} 本段程序仅仅考虑Shared memory的应用
Example:1D Stencil
一维数组元素，每个输出元素是半径内的输入元素的总和，半径为３。
__global__ void stencil_1d(int *in, int *out)
{
__shared__ int temp[BLOCK_SIZE + 2 * RADIUS];
int gindex = threadIdx.x + blockIdx.x * blockDim.x + RADIUS; //查找指针真正的位置,即索引
int lindex = threadIdx.x + RADIUS;
// Read input elements into shared memory
temp[lindex] = in[gindex];
if (threadIdx.x < RADIUS)
{
temp[lindex – RADIUS] = in[gindex – RADIUS];
temp[lindex + BLOCK_SIZE] = in[gindex + BLOCK_SIZE];
}
// Synchronize (ensure all the data is available)
__syncthreads(); //确保上一步操作完成
// Apply the stencil
int result = 0;
for (int offset = -RADIUS ; offset <= RADIUS ; offset++)
result += temp[lindex + offset];
// Store the result
out[gindex - RADIUS] = result; //存储结果
} 考虑到该函数和周边像素求和，定义了RAIDUS，注意该种定义对源程序的影响。
Example: edge detection algorithm
下面一段程序为host调用GPU内核函数，故函数声明为__global__,并且涉及二维操作，将block和thread定义为二维的。
dim3 blocksPerGrid ( N/TILE_WIDTH , N/TILE_WIDTH , 1); //共三维，这里使用了两维
dim3 threadsPerBlock(TILE_WIDTH,TILE_WIDTH,1); //同样的，thread定义为两维
inverseEdgeDetect2D<<< blocksPerGrid, threadsPerBlock >>>(d_output, d_input, d_edge); //定义GPU内核函数的参数 GPU内核函数，由于是二维的，注意Shared memory的使用，并且考虑到该函数和周边像素做差分，定义了RAIDUS，注意该种定义对源程序的影响。注意一维和二维的联系与区别
__global__ void inverseEdgeDetect2D(float *d_output, float *d_input, \
float *d_edge)
{
__shared__ float ds_input[TILE_WIDTH+ 2*RAIDUS][TILE_WIDTH+2*RAIDUS];
int bx = blockIdx.x; int by = blockIdx.y;
int tx = threadIdx.x + RAIDUS; int ty = threadIdx.y + RAIDUS;
int Row = by * blockDim.y + ty ;
int Col = bx * blockDim.x + tx ;
int idx;
int numcols = N + 2;
// Read input elements into shared memory
ds_input[ty][tx] = d_input[Row*numcols + Col];
if (threadIdx.x < RAIDUS)
{
ds_input[ty][tx - RAIDUS] = d_input[Row*numcols + Col- RAIDUS];
ds_input[ty][tx + TILE_WIDTH] = d_input[Row*numcols + Col + TILE_WIDTH];
}
if (threadIdx.y < RAIDUS)
{
ds_input[ty- RAIDUS][tx] = d_input[(Row- RAIDUS)*numcols + Col];
ds_input[ty + TILE_WIDTH][tx] = d_input[(Row + TILE_WIDTH)*numcols + Col];
}
// Synchronize (ensure all the data is available)
__syncthreads();
// Synchronize (ensure all the data is available)
idx = Row * numcols + Col;
d_output[idx]=(ds_input[ty+RAIDUS][tx] + ds_input[ty-RAIDUS][tx] \
+ ds_input[ty][tx+RAIDUS] + ds_input[ty][tx-RAIDUS] - \
d_edge[idx]) * 0.25;
} Memory Coalescing
在访问全局内存时，当所有线程连续访问数据的内存时，会有峰值的性能利用率。因此，要尽可能连续访问内存。
扫一扫
4256
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
