并行系统学习之路（三） ---- CUDA学习
NVIDIA CUDA 编程学习
本篇我们来学习CUDA编程，首先我们明确一下几个名词的相互关系：
-
GPU:
- Graphics processing unit
- GPU和CPU一样也是计算机的计算单位，但是他们的设计目标有很大不同，GPU专注于计算，它有有更多的ALU，而CPU为了处理各种计算机指令需要许多其他的复杂结构
- GPU适合大规模并行计算，而CPU适合逻辑控制串行运算
- CPU与GPU架构对比
-
GPGPU：
- General Purpose computation using GPU in applications(Other than 3D graphics)
- 通常我们说起GPU，总是会想起3D图形渲染相关的事物，其实我们也可以运用GPU的并行计算特性可以完成很多其他的工作，我们把这种方法称为GPGPU
-
CUDA
- Compute Unified Device Architecture
- CUDA是一个GPU编程模型，相当于一个开发者与GPU之间的接口，开发者使用CUDA可以写出在NVIDIA GPU上运行的并行程序
- CPU称为Host，GPU称为Device
CUDA编程模型
- 相当于CPU的协作处理器
- 拥有自己的DRAM
- SIMD架构
- GPU线程开销小，可以轻松执行多个线程
线程分配
- 每一个Kernel函数将会在包含多个线程block的grid中执行
- 一个block中的所有线程共享memory
- 一个block中的线程可以通过以下方法互相交流
__syncthreads()同步操作- 通过shared memory进行数据分享
- 两个不同block中的线程不能相互交流（除非使用全局memory，但是比较慢）
- 线程和block都有属于自己的ID
CUDA中的三种函数
| 执行地点 | 调用地点 | |
|---|---|---|
__device__ float DeviceFunc() |
device | device |
__global__ void KernelFunc() |
device | host |
__host__ float HostFunc() |
host | host |
CUDA 模型下的内存空间分配
- 每个线程可以
- R/W per-thread registers
- R/W per-thread local memory
- R/W per-block shared memory
- R/W per-grid global memory
- Read only per-grid constant memory
- Read only per-grid texture memory
- host只可以读写global，constant以及texture的memory
- Global memory
- host与device之间读写数据的主要方式
- 所有线程可见，只读
- Texture and Constant Memories
- host初始化的常量
- 所有线程可见，只读
调用Kernel函数实例
调用Kernel函数是必须指定config：
__global__ void KernelFunc(...);
dim3 DimGrid(100, 50); // 5000 thread blocks
dim3 DimBlock(4, 8, 8); // 256 threads per block size_t SharedMemBytes = 64; // 64 bytes of shared memory KernelFunc<<< DimGrid, DimBlock, SharedMemBytes >>>(...);
- 所有Kernel函数调用都是异步的（除非显式堵塞）
- Kernel函数允许递归
下面是一个实例，该实例把一个数组里的所有值自增了1
__global__ void incrArrayOnDevice(float *a, int N)
{
int idx = blockIdx.x*blockDim.x + threadIdx.x;
if (idx<N) a[idx] = a[idx]+1.f;
}
int main() {
float *a_h, *a_d; int i, N=10;
size_t size = N*sizeof(float);
a_h = (float *)malloc(size);
for (i=0; i<N; i++) a_h[i] = (float)i;
// allocate array on device
cudaMalloc((void **) &a_d, size);
// copy data from host to device
cudaMemcpy(a_d, a_h, sizeof(float)*N, cudaMemcpyHostToDevice);
// do calculation on device:
// Part 1 of 2. Compute execution configuration
int blockSize = 4;
int nBlocks = N/blockSize + (N%blockSize == 0?0:1);
// Part 2 of 2. Call incrementArrayOnDevice kernel
incrArrayOnDevice <<< nBlocks, blockSize >>> (a_d, N);
// Retrieve result from device and store in b_h
cudaMemcpy(a_h, a_d, sizeof(float)*N, cudaMemcpyDeviceToHost);
// cleanup
free(a_h);
cudaFree(a_d);
}
硬件执行模型
- 在grid里，每一个block会被分成warp，每一个warp将由一个多处理器(SM)执行
- device每次只执行一个grid
- 每一个block将由一个多处理器执行
- 因此shared memory空间将位于on-chip shared memory，每个block可以共享
- 一个多处理器可以同时执行多个block
- Shared memory 和寄存器空间是由根据所有同时运行的block来进行分配的
- 因此减少shared memory和寄存器的使用(一个block中的)会增加可以同时运行的block的数量
Thread, Warps, Block 之间关系
- 在一个warp中最多可以有32个线程(小于等于32)
- 一个block中最多有32个warp
- 每个block在一个SM中执行(因此每个warp也一样)
- GF110有16个SMs
- 至少16个block可以保证”填满”device
- 越多效率提升越高
- 如果有足够资源(寄存器，线程空间，Shared memory),一个SM可以执行多个block
总结
- device = GPU = 所有SM
- 多处理器 = 许多处理器和对应shared memory的组合
- Kernel = GPU 程序
- Grid = 一组执行kernel的block
- Block = 一组SIMD的线程，他们可以执行kernel function以及通过shared memory交流
实例 – 使用Leibniz公式求π值
#include <stdio.h>
#include <math.h>
#include "mytime.h"
#define THREADS 512
#define MAX_BLOCKS 64
// find the nearest 2 power
__device__ int NearestPowerOf2 (int n) {
if (!n) return n; //(0 == 2^0)
int x = 1;
while(x < n)
{
x <<= 1;
}
return x;
}
// GPU kernel, we know: THREADS == blockDim.x
__global__ void integrate(int *n, int *blocks, double *gsum) {
const unsigned int bid = blockDim.x * blockIdx.x + threadIdx.x;
const unsigned int tid = threadIdx.x;
double sum;
int i;
int nTotalThreads;
__shared__ double ssum[THREADS];
sum = 0.0;
for (i = bid ; i < *n; i += blockDim.x * *blocks) {
sum += powf(-1, i) / (2*i + 1);
}
ssum[tid] = sum * 4;
// block reduction
__syncthreads();
nTotalThreads = NearestPowerOf2(blockDim.x/2);
for (i = nTotalThreads; i > 0; i >>= 1) { /* per block */
if (tid < i)
if(tid+i < blockDim.x)
ssum[tid] += ssum[tid + i];
__syncthreads();
}
gsum[bid] = ssum[tid];
}
// number of threads must be a power of 2
__global__ static void global_reduce(int *n, int *blocks, double *gsum)
{
__shared__ double ssum[THREADS];
const unsigned int tid = threadIdx.x;
unsigned int i;
int nTotalThreads;
if (tid < *blocks)
ssum[tid] = gsum[tid * THREADS];
else
ssum[tid] = 0.0;
__syncthreads();
nTotalThreads = NearestPowerOf2(blockDim.x/2);
for (i = nTotalThreads; i > 0; i >>= 1) { /* per block */
if (tid < i)
if(tid+i < blockDim.x)
ssum[tid] += ssum[tid + i];
__syncthreads();
}
gsum[tid] = ssum[tid];
}
int main(int argc, char *argv[]) {
int n, blocks;
int *n_d, *blocks_d; // device copy
double PI25DT = 3.141592653589793238462643;
double pi;
double *mypi_d; // device copy of pi
struct timeval startwtime, endwtime, diffwtime;
// Allocate memory on GPU
cudaMalloc( (void **) &n_d, sizeof(int) * 1 );
cudaMalloc( (void **) &blocks_d, sizeof(int) * 1 );
cudaMalloc( (void **) &mypi_d, sizeof(double) * THREADS * MAX_BLOCKS );
while (1) {
printf("Enter the number of intervals: (0 quits) ");fflush(stdout);
scanf("%d",&n);
printf("Enter the number of blocks: (<=%d) ", MAX_BLOCKS);fflush(stdout);
scanf("%d",&blocks);
gettimeofday(&startwtime, NULL);
if (n == 0 || blocks > MAX_BLOCKS)
break;
// copy from CPU to GPU
cudaMemcpy( n_d, &n, sizeof(int) * 1, cudaMemcpyHostToDevice );
cudaMemcpy( blocks_d, &blocks, sizeof(int) * 1, cudaMemcpyHostToDevice );
integrate<<< blocks, THREADS >>>(n_d, blocks_d, mypi_d);
if (blocks > 1)
global_reduce<<< 1, 512 >>>(n_d, blocks_d, mypi_d);
// copy back from GPU to CPU
cudaMemcpy( &pi, mypi_d, sizeof(double) * 1, cudaMemcpyDeviceToHost );
gettimeofday(&endwtime, NULL);
MINUS_UTIME(diffwtime, endwtime, startwtime);
printf("pi is approximately %.16f, Error is %.16f\n",
pi, fabs(pi - PI25DT));
printf("wall clock time = %d.%06d\n",
diffwtime.tv_sec, diffwtime.tv_usec);
}
// free GPU memory
cudaFree(n_d);
cudaFree(mypi_d);
return 0;
}