基于Cuda开发GPUGPU程序时，最重要的仍然是内核的设计，这是Cuda性能优化的难点，提供了不少岗位，养活了一大批工程师。这里以一个相对简单的的求平方和算法为例，从编程和优化，调试几个维度，介绍利用cuda开发并行计算程序时的关注点。
cuda API
NVIDIA CUDA计算架构为开发者提供了三个层面的API，分别是Cuda Lib, Cuda RT, 和cuda driver。cuda driver是比较底层的API，用法复杂但是性能高，可以深度二次优化，对于研发能力强的用户可以在这个层次上做出高性能的计算方案出来，其次是最常用的cuda runtime，也就是我们常用的cuda API.最上层是cudalib,cudalib提供给研发能力一般，希望快速上手的开发者，包含各类已经预先开发好的数学库和数学函数。
cuda编程模型
在CUDA编程模型中引入主机端和设备端的概念，CPU是主机端，GPU属于设备端，主机端仅有一个，而设备端可以同时有很多（比如NVLINK 8卡互联），CPU负责复杂逻辑处理和运算量少的计算，而GPU负责运行简单但是计算量大的并行计算。
一个kernel函数对应一个grid,每个grid根据需要配置不同的block数量和thread数量。从编程模型可以看出，cuda包含三个逻辑层，grid, block和thread. 一个GRID下的线程共享全局的内存空间。GRID是线程结构的第一层次，GRID下面分BLOCK，一个BLOCK包含很多个线程，线程是第三层次。
CUDA下，一个grid内只能跑一个KERNEL，KERNEL和GRID是一对一的关系，GRID是kernel launch的参数。GPU可以并行Launch多个grid.
warp
GPU编程的一个基本特点是大规模并行，让GPU内数千计的微处理器同时转向要处理的数据，每个线程处理一个数据元素，由SIMT模拟的SIMD（相比较SIMD，SIMT偏向于灵活性而损失了一部分性能）。
一方面是大量需要执行的任务，另一方面是很多等待任务的微处理器，如何让这么多的微处理器有条不紊的把所有任务都执行完毕呢？这里涉及到一个调度粒度的概念。
与军队里把士兵分成一个个小的战斗单位类似，在CUDA中，也把微处理器分成一个个小组。每个组的大小是一样的。NVIDIA 分组的粒度是32个。CUDA个这个组取了个特别的名字：Warp.
在GPU编程中，WARP（线程束）是指一组共同协作的线程，通常为32个线程。这些线程执行相同的指令，但对不同的数据进行计算。这些线程在同一个WARP中会被划分到同一个流处理器中，共享同一个指令单元和寄存器，以实现高效的并行计算。
WARP是GPU中的基本计算单元，因为GPU硬件中的每个流处理器都包含多个WARP处理器。WARP是由GPU硬件自动管理的，程序员通常只需要考虑如何将线程划分为WARP，并让它们协同工作以实现高效的并行计算。
Warp是GPU调度的基本单位，这意味着，当GPU调度硬件资源时，一次分派的执行单元至少是32个，如果每个线程的块大小不足32个，那么也会分配32个，多余的硬件单元处于闲置状态。
Warp一词来源于纺织机，纺织机的核心是织机（Loom),历经书千年的发展，世界各地的人们发明了很多织机，虽然种类很多，但是大多数织机的一个基本原理都是让经线和纬线交织到一起，通常的做法是首先部署好一组经线，然后把系着纬线的梭子穿过经线，如此往复。Warp就是经线。
在纺织中，经线的数量决定了织物的幅度，也可以认为经线的数量决定了并行操作的并行度，在CUDA中，使用Warp来代表同时操作的一批线程，也代表并行度。
在获取数据之后，在SM中以32个线程为一组的线程束(Warp)来调度，来开始处理顶点数据。Warp是典型的单指令多线程（SIMT，SIMD单指令多数据的升级）的实现，也就是32个线程同时执行的指令是一模一样的，只是线程数据不一样，这样的好处就是一个warp只需要一个套逻辑对指令进行解码和执行就可以了，芯片可以做的更小更快，之所以可以这么做是由于GPU需要处理的任务是天然并行的。
Warp调度器会按照顺序分发指令给整个warp，单个warp中的线程会锁步(lock-step)执行各自的指令，如果线程碰到不激活执行的情况也会被遮掩(be masked out)。被遮掩的原因有很多，例如当前的指令是if(true)的分支，但是当前线程的数据的条件是false，或者循环的次数不一样（比如for循环次数n不是常量，或被break提前终止了但是别的还在走），因此在shader中的分支会显著增加时间消耗，在一个warp中的分支除非32个线程都走到if或者else里面，否则相当于所有的分支都走了一遍，线程不能独立执行指令而是以warp为单位，而这些warp之间才是独立的。
Warp中的指令可以被一次完成，也可能经过多次调度，例如通常SM中的LD/ST(加载存取)单元数量明显少于基础数学操作单元。
由于某些指令比其他指令需要更长的时间才能完成，特别是内存加载，warp调度器可能会简单地切换到另一个没有内存等待的warp，这是GPU如何克服内存读取延迟的关键，只是简单地切换活动线程组。为了使这种切换非常快，调度器管理的所有warp在寄存器文件中都有自己的寄存器。这里就会有个矛盾产生，shader需要越多的寄存器，就会给warp留下越少的空间，就会产生越少的warp，这时候在碰到内存延迟的时候就会只是等待，而没有可以运行的warp可以切换。
In a Graphics Processing Unit (GPU), a warp is a group of threads that are executed together as a single unit. The concept of a warp is essential to understanding how GPUs process data in parallel.
Here's a breakdown of what a warp is:
Thread group: A warp is a group of 32 threads (in most modern GPUs) that are executed together. These threads are called a thread group or a warp group.
Thread execution: Each thread within a warp executes the same instruction, but with different data. This is known as single-instruction, multiple-data (SIMD) execution.
Warp scheduling: The GPU schedules warps to execute in a specific order. The scheduler decides which warp to execute next, based on factors like dependencies between warps, memory access patterns, and available resources.
Warp execution: When a warp is executed, all 32 threads within the warp execute the same instruction simultaneously. This is known as a warp execution cycle.
Warp divergence: If threads within a warp take different execution paths (e.g., due to conditional statements), the warp is said to diverge. In this case, the GPU must execute each thread independently, which can lead to reduced performance.
The benefits of warps include:
Parallelism: Warps enable massive parallelism, allowing GPUs to process large amounts of data in parallel.
Resource sharing: Warps share resources like registers, memory, and execution units, which improves resource utilization and reduces overhead.
Cache efficiency: Warps can efficiently use the GPU's cache hierarchy, reducing memory access latency.
In summary, a warp is a group of threads that execute the same instruction in parallel, allowing GPUs to achieve massive parallel processing and efficient resource utilization.算法描述
平方和算法是一种缩减算法，缩减算法指的是从多个数据中提炼出较少的数据的一类算法，在统计中求和，找最值，均值，和方差等应用中，以及在图像处理中求一副图像的总亮度等，都是缩减算法(reduction)。公式为：
nvidia reduction 示意图：
CUDA并行编程方法
得益于数量巨大的核心数量，GPU具有强大的并行计算能力，但是它的局限性也很明显，GPU从单核的结构和ISA性能上讲，计算能力远不如CPU，优势是胜在核多。CPU有复杂的存储器缓冲系统，先进的指令缓存系统和强大的分支预测能力。而GPU中的标量处理器结构相对简单，甚至都不要求是图灵完备的（早期的GPU甚至都不支持条件分支）。CPU支持顺序执行，高效循环和跳转，而GPU相对简单的结构使它较适合处理顺序的，单一的，少循环，少跳转的语句。所以在由GPU和CPU构成的异构系统中，GPU不能独立运行。CUDA编程作为一种实现也不例外。CUDA开发的典型模式为：首先由主机分配主机端和设备端的内存，之后再将计算数据传输给GPU侧，调用设备端核函数得到运行结果，在进行设备到主机端的数据传输将结果回传给主机侧，如下图所示：
GPU Programing
GPU Programing 层次结构 as below:
sp(streaming processor) : 最基本的处理单元，最后具体的指令和任务都是在sp上处理的。GPU进行并行计算，也就是很多个sp同时做处理。sp 就是 cuda core.
sm(streaming multiprocessor):多个sp加上其他的一些资源组成一个sm。
warp:GPU执行程序时的调度单位，目前cuda一个warp有32个thread，同在一个warp的线程，以不同数据资源执行相同的指令。
grid、block、thread：在利用cuda进行编程时，一个grid分为多个block，而一个block分为多个thread.其中任务划分到是否影响最后的执行效果。划分的依据是任务特性和GPU本身的硬件特性。
一个BLOCK只能在一个SM上被调度，SM一般可以调度多个线程块，这要看SM本身的能力。进行划分时，最好保证每个block里的warp比较合理，那样可以一个sm可以交替执行里面的warp，从而提高效率，此外，在分配block时，要根据GPU的sm个数，分配出合理的block数，让GPU的sm都利用起来，提利用率。分配时，也要考虑到同一个线程block的资源问题，不要出现对应的资源不够。
一个KERNEL的线程块可能被分配到一个SM上排队调度，也可能分配到多个SM上同时执行，GRID是逻辑层，而SM才是物理层。BLOCK竞争SM的的资源。SM采用SIMT方式执行。
一个KERNEL的所有线程在物理层面并不一定是同时执行的，而是按照BLOCK排队处理。SM的执行单元是WARP，一个WARP 32个线程，所以BLOCK大小一般要设置为32的倍数。
switch between different thread and block to see the warp change
初始化cuda设备：
初始化cuda设备的代码如下，通过这一步，可以获取到cuda设备的warp size，核心频率，内存大侠以及grid,block维数等信息，对于后续的调试调优有重要意义。
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
static cudaDeviceProp prop;
void init_cuda(void)
{
int count, dev;
int i;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "there is no cuda device.\n");
return;
} else {
cudaGetDevice(&dev);
fprintf(stdout,"there are %d cudda device found, id %d.\n", count, dev);
}
for(i = 0; i < count; i ++) {
printf("===============================Device %d==================================\n", i);
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
printf("%s\n", prop.name);
printf("Total global memory: %ld Bytes\n", prop.totalGlobalMem);
printf("Max shareable memory per block %ld Bytes.\n", prop.sharedMemPerBlock);
printf("Maximum registers per block: %d\n", prop.regsPerBlock);
printf("Wrap Size %d.\n", prop.warpSize);
printf("Maximum threads per block %d.\n", prop.maxThreadsPerBlock);
printf("Maximum block dimensions [%d, %d, %d].\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("Maximum grid dimensions [%d, %d, %d].\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("Total constant memory: %ld.\n", prop.totalConstMem);
printf("Support compute Capability: %d.%d.\n", prop.major, prop.minor);
printf("Kernel Frequency %d kHz.\n", prop.clockRate);
printf("Number of MultProcessors %d.\n", prop.multiProcessorCount);
printf("Is MultiGPU: %s.\n", prop.isMultiGpuBoard ? "Yes" : "No");
printf("L2 Cache Size: %d Bytes.\n", prop.l2CacheSize);
printf("Memory Bus Width: %d.\n", prop.memoryBusWidth);
printf("ECC status: %s.\n", prop.ECCEnabled? "Enable" : "Disable");
}
printf("=========================================================================\n");
}
cudaSetDevice(1);
}
int main(void)
{
init_cuda();
return 0;
}运行结果如下：
此信息和通过gpu-z工具解析到的数据是一致的，可以参考博客：
cuda-z/gpu-z/cpu-z工具分析GPU显卡和CPU算力信息_papaofdoudou的博客-优快云博客_cuda-z
kernel运行时间统计
Kernel运行时间有GPU计时和事件记时两种，顾名思义，GPU记时是由设备端执行计时函数记录时间，相应的函数是clock. CUDA架构的GPU每个多处理器中有一个计数器，用于对核心时钟进行采样计数，clock返回的就是核心频率的计数值。通过计算kernel运行结束和运行开始的时间差值，在除以上一步得到的核心频率，即为kerne时间的运行时间。
由于设备端会同时调度多个Warper同时跑，而每个Warp由32个线程（第一步已经获取到打印出来），所以用clock函数得到的是线程从开始执行内核到执行结束所消耗的时间，并不是线程实际的指令执行时间。
Event方式要通过调用几组cuda API实现：
cudaEventCreate
cudaEventRecord
cudaEventSynchronize
cudaEventElapsedTime
cudaEventDestroy代码实现：
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#define DATA_SIZE (100*1024*1024)
#define BLOCK_NUM 32
#define THREAD_NUM 256
void generate_numbers(int *pnum, int size)
{
int i;
for(i = 0; i < size; i ++)
{
//pnum[i] = rand() % 10;
pnum[i] = 1;
}
}
__global__ static void sum_of_squares(int *pnum, int *pres, clock_t *pclock)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
if(tid == 0) pclock[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM)
{
shared[tid] += pnum[i] * pnum[i];
}
__syncthreads();
if(tid == 0) {
for(i = 1; i < THREAD_NUM; i ++){
shared[0] += shared[i];
}
pres[bid] = shared[0];
}
if(tid == 0) pclock[bid + BLOCK_NUM] = clock();
}
__global__ static void sum_of_squares_eff(int *pnum, int *pres, clock_t *clock_time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
int offset = 1;
int mask = 1;
if(tid == 0) clock_time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM)
{
shared[tid] += pnum[i] * pnum[i];
}
__syncthreads();
while(offset < THREAD_NUM) {
if((tid & mask) == 0){
shared[tid] += shared[tid + offset];
}
offset += offset;
mask += offset;
__syncthreads();
}
if(tid == 0) {
pres[bid] = shared[0];
clock_time[bid + BLOCK_NUM] = clock();
}
}
__global__ static void sum_of_squares_eff_2(int *pnum, int *pres, clock_t *clock_time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
int offset = THREAD_NUM/2;
if(tid == 0) clock_time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM)
{
shared[tid] += pnum[i] * pnum[i];
}
__syncthreads();
while(offset > 0) {
if(tid < offset){
shared[tid] += shared[tid + offset];
}
offset >>= 1;
__syncthreads();
}
if(tid == 0) {
pres[bid] = shared[0];
clock_time[bid + BLOCK_NUM] = clock();
//printf("%s line %ld\n", __func__, clock());
}
}
__global__ static void sum_of_squares_eff_intri(int *pnum, int *pres, clock_t *clock_time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i = 0;
if(tid == 0) clock_time[bid] = clock();
shared[tid] = 0;
#pragma unroll
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += __mul24(BLOCK_NUM, THREAD_NUM))
{
shared[tid] += __mul24(pnum[i], pnum[i]);
}
__syncthreads();
if(tid < 128) {
shared[tid] += shared[tid + 128];
}
__syncthreads();
if(tid < 64) {
shared[tid] += shared[tid + 64];
}
__syncthreads();
// because in the same wrap, the below thread are not need sync like above.
if(tid < 32) shared[tid] += shared[tid + 32];
__syncthreads();
if(tid < 16) shared[tid] += shared[tid + 16];
__syncthreads();
if(tid < 8) shared[tid] += shared[tid + 8];
__syncthreads();
if(tid < 4) shared[tid] += shared[tid + 4];
__syncthreads();
if(tid < 2) shared[tid] += shared[tid + 2];
__syncthreads();
if(tid < 1) shared[tid] += shared[tid + 1];
__syncthreads();
if(tid == 0) {
pres[bid] = shared[0];
clock_time[bid + BLOCK_NUM] = clock();
}
}
static cudaDeviceProp prop;
void init_cuda(void)
{
int count, dev;
int i;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "there is no cuda device.\n");
return;
} else {
cudaGetDevice(&dev);
fprintf(stdout,"there are %d cudda device found, id %d.\n", count, dev);
}
for(i = 0; i < count; i ++) {
printf("===============================Device %d==================================\n", i);
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
printf("%s\n", prop.name);
printf("Total global memory: %ld Bytes\n", prop.totalGlobalMem);
printf("Max shareable memory per block %ld Bytes.\n", prop.sharedMemPerBlock);
printf("Maximum registers per block: %d\n", prop.regsPerBlock);
printf("Wrap Size %d.\n", prop.warpSize);
printf("Maximum threads per block %d.\n", prop.maxThreadsPerBlock);
printf("Maximum block dimensions [%d, %d, %d].\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("Maximum grid dimensions [%d, %d, %d].\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("Total constant memory: %ld.\n", prop.totalConstMem);
printf("Support compute Capability: %d.%d.\n", prop.major, prop.minor);
printf("Kernel Frequency %d kHz.\n", prop.clockRate);
printf("Number of MultProcessors %d.\n", prop.multiProcessorCount);
printf("Is MultiGPU: %s.\n", prop.isMultiGpuBoard ? "Yes" : "No");
printf("L2 Cache Size: %d Bytes.\n", prop.l2CacheSize);
printf("Memory Bus Width: %d.\n", prop.memoryBusWidth);
printf("ECC status: %s.\n", prop.ECCEnabled? "Enable" : "Disable");
}
printf("=========================================================================\n");
}
cudaSetDevice(1);
}
#define USE_HOST_ALLOC
int main(void)
{
int *pdata, *psum;
clock_t *pclock;
int *pgpudata, *pres;
clock_t *pclock_t;
init_cuda();
cudaMallocHost((void**)&pdata, DATA_SIZE * sizeof(int));
if(pdata == NULL) {
fprintf(stderr, "malloc host buffer failure.\n");
return -1;
}
generate_numbers(pdata, DATA_SIZE);
cudaMallocHost((void**)&psum, BLOCK_NUM * sizeof(int));
if(psum == NULL) {
fprintf(stderr, "malloc host buffer failure.\n");
return -1;
}
memset(psum, 0x00, BLOCK_NUM * sizeof(int));
cudaMallocHost((void**)&pclock, sizeof(clock_t) * BLOCK_NUM * 2);
if(pclock == NULL) {
fprintf(stderr, "malloc host buffer failure.\n");
return -1;
}
memset(pclock, 0x00, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMalloc((void**)&pgpudata, sizeof(int) * DATA_SIZE);
if(pgpudata == NULL) {
fprintf(stderr, "malloc device buffer failure.\n");
return -1;
}
#ifndef USE_HOST_ALLOC
cudaMalloc((void**)&pres, sizeof(int) * BLOCK_NUM);
#else
cudaHostGetDevicePointer((void**)&pres, (void*)psum, 0);
#endif
if(pres == NULL) {
fprintf(stderr, "malloc device buffer failure.\n");
return -1;
}
cudaMalloc((void**)&pclock_t,sizeof(clock_t) * BLOCK_NUM * 2);
if(pclock_t == NULL) {
fprintf(stderr, "malloc device buffer failure.\n");
return -1;
}
cudaMemcpy(pgpudata, pdata, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
float tm;
cudaEvent_t start,stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
//sum_of_squares<<< BLOCK_NUM, THREAD_NUM, THREAD_NUM * sizeof(int) >>>(pgpudata, pres, pclock_t);
//sum_of_squares_eff<<< BLOCK_NUM, THREAD_NUM, THREAD_NUM * sizeof(int) >>>(pgpudata, pres, pclock_t);
//sum_of_squares_eff_2<<< BLOCK_NUM, THREAD_NUM, THREAD_NUM * sizeof(int) >>>(pgpudata, pres, pclock_t);
sum_of_squares_eff_intri<<< BLOCK_NUM, THREAD_NUM, THREAD_NUM * sizeof(int) >>>(pgpudata, pres, pclock_t);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&tm,start,stop);
cudaEventDestroy(start);
cudaEventDestroy(stop);
printf("GPU Elapsed time:%.6f ms.\n",tm);
cudaDeviceSynchronize();
//cudaThreadSynchronize();
#ifndef USE_HOST_ALLOC
cudaMemcpy(psum, pres, sizeof(int) * BLOCK_NUM, cudaMemcpyDeviceToHost);
#else
;
#endif
cudaMemcpy(pclock, pclock_t, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost);
int final_result = 0;
int i;
for(i = 0; i < BLOCK_NUM; i ++)
{
//printf("%s line %d psum[%d] = %d.\n", __func__, __LINE__, i, psum[i]);
final_result += psum[i];
}
clock_t clock_min,clock_max;
clock_min = pclock[0];
clock_max = pclock[BLOCK_NUM];
for(i = 0; i < BLOCK_NUM; i ++)
{
//printf("%s line %d %ld.\n", __func__, __LINE__, pclock[i]);
if(clock_min > pclock[i])
clock_min = pclock[i];
if(clock_max < pclock[i + BLOCK_NUM])
clock_max = pclock[i + BLOCK_NUM];
}
printf("%s line %d sum %d, time %ld ms.\n", __func__, __LINE__, final_result, ((clock_max - clock_min) * 1000)/(prop.clockRate * 1000));
cudaFree(pclock_t);
#ifndef USE_HOST_ALLOC
cudaFree(pres);
#endif
cudaFree(pgpudata);
cudaFreeHost(pdata);
cudaFreeHost(psum);
cudaFreeHost(pclock);
return 0;
}运行结果，可以看到，计算100*1024*1024组数据的平方和需要大约9ms,用两种方式统计得到的kernel时间相差无几。
计算包含100M个int类型元素的向量的平方和(sum_of_squares_eff_2函数的实现)在GPU上的执行过程如下图所示：
- group_size = BLOCK_NUM*THREAD_NUM;
- group_start(groupid) = data + groupid *group_size
- group_id = idx/(group_size);
- offst_in_group = idx%group_size;
- blockid = offset_in_group/THREAD_NUM;
- threadid = offset_in_group%THREAD_NUM;
- item1 ≡ item2 mod group_size,则item1和item2在同一个BLOCK的同一个线程上计算。
- task_warp_num = 256(thread)/32 x 32(BLOCK_NUM) = 256 warp,1个SM的最大吞吐是64个WARP，所以需要1个WARP运行4轮，或者4个WARP运行一轮。
- 每个BLOCK有THREAD_NUM/32=256/32=8个warp,任务被分成BLOCK_NUM(看SM容量），依次运行，WARPN内不需要同步，同一个BLOCK内的WARP之间需要BAR同步。
- BLOCK_NUM*THREAD_NUM为一组，它表示的是整个GRID的线程数，向量的元素很多的情况下，按照这个数分组，相当于使用了多个GRID处理，这是一种grid-stride loop方法。
关于实现中的__syncthreads同步，它保证了BLOCK中所有的线程都执行到同一位置，由于WARP中的线程总是同时被激活且一同分配执行任务，因此WARP内的线程不需要同步，同一个BLOCK内可能有多个WARP，__syncthreads进行的是属于同一个BLOCK内的WARP之间的同步。
cudaDeviceSynchronize
在CUDA当中，核函数kernel的执行总是异步的，而cudaMemcpy数据传输总是同步的。特别需要注意的是主机在提交核函数之后，不会阻塞等待核函数执行完毕。运行KERNEL时，一定要记得添加cudaDeviceSynchronize() 同步，或者添加一个数据传输(cudaMemcpy-隐含着同步操作) ，以保证核函数执行结束。证明这个问题不难，统计核函数和cudaDeviceSynchronize函数的执行时间来看，前者执行时间和KERNEL处理的数据量之间没有关系，执行时间为固定长度，而后者随着数据量的增加而增大：
正是因为KERNEL调用是异步调用，所以才有Stream的概念，Stream内是顺序执行的，而Stream之间的顺序是乱序的。
#include "cuda_runtime.h"
#include <iostream>
#include <stdio.h>
#include <math.h>
#define N (1024*1024)
#define FULL_DATA_SIZE N*20
__global__ void kernel(int *a, int *b, int *c)
{
int threadID = blockIdx.x * blockDim.x + threadIdx.x;
if (threadID < N) {
c[threadID] = (a[threadID] + b[threadID]) / 2;
}
}
int main(void)
{
//获取设备属性
cudaDeviceProp prop;
int deviceID;
cudaGetDevice(&deviceID);
cudaGetDeviceProperties(&prop, deviceID);
//检查设备是否支持重叠功能
if (!prop.deviceOverlap) {
printf("No device will handle overlaps. so no speed up from stream.\n");
return 0;
}
//创建一个CUDA流
cudaStream_t stream;
cudaStreamCreate(&stream);
int *host_a, *host_b, *host_c;
int *dev_a, *dev_b, *dev_c;
//在GPU上分配内存
cudaMalloc((void **)&dev_a, N * sizeof(int));
cudaMalloc((void **)&dev_b, N * sizeof(int));
cudaMalloc((void **)&dev_c, N * sizeof(int));
//在CPU上分配页锁定内存
cudaHostAlloc((void **)&host_a, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
cudaHostAlloc((void **)&host_b, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
cudaHostAlloc((void **)&host_c, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
//主机上的内存赋值
for (int i = 0; i < FULL_DATA_SIZE; i++) {
host_a[i] = i;
host_b[i] = FULL_DATA_SIZE - i;
}
//启动计时器
cudaEvent_t start, stop;
float elapsedTime;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
for (int i = 0; i < FULL_DATA_SIZE; i += N) {
cudaMemcpyAsync(dev_a, host_a + i, N * sizeof(int), cudaMemcpyHostToDevice, stream);
cudaMemcpyAsync(dev_b, host_b + i, N * sizeof(int), cudaMemcpyHostToDevice, stream);
kernel <<< N / 1024, 1024, 0, stream >>> (dev_a, dev_b, dev_c);
cudaMemcpyAsync(host_c + i, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost, stream);
}
// wait until gpu execution finish
cudaStreamSynchronize(stream);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&elapsedTime, start, stop);
std::cout << "消耗时间： " << elapsedTime << std::endl;
//输出前10个结果
for (int i = 0; i < 10; i++) {
std::cout << host_c[i] << std::endl;
}
getchar();
// free stream and mem
cudaFreeHost(host_a);
cudaFreeHost(host_b);
cudaFreeHost(host_c);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);
cudaStreamDestroy(stream);
return 0;
}
使用两个流，速度加快
#include "cuda_runtime.h"
#include <iostream>
#include <stdio.h>
#include <math.h>
#define N (1024*1024)
#define FULL_DATA_SIZE N*20
__global__ void kernel(int *a, int *b, int *c)
{
int threadID = blockIdx.x * blockDim.x + threadIdx.x;
if (threadID < N) {
c[threadID] = (a[threadID] + b[threadID]) / 2;
}
}
int main()
{
//获取设备属性
cudaDeviceProp prop;
int deviceID;
cudaGetDevice(&deviceID);
cudaGetDeviceProperties(&prop, deviceID);
//检查设备是否支持重叠功能
if (!prop.deviceOverlap) {
printf("No device will handle overlaps. so no speed up from stream.\n");
return 0;
}
//创建两个CUDA流
cudaStream_t stream1, stream2;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);
int *host_a, *host_b, *host_c;
int *dev_a, *dev_b, *dev_c;
int *dev_a1, *dev_b1, *dev_c1;
//在GPU上分配内存
cudaMalloc((void **)&dev_a, N * sizeof(int));
cudaMalloc((void **)&dev_b, N * sizeof(int));
cudaMalloc((void **)&dev_c, N * sizeof(int));
cudaMalloc((void **)&dev_a1, N * sizeof(int));
cudaMalloc((void **)&dev_b1, N * sizeof(int));
cudaMalloc((void **)&dev_c1, N * sizeof(int));
//在CPU上分配页锁定内存
cudaHostAlloc((void **)&host_a, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
cudaHostAlloc((void **)&host_b, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
cudaHostAlloc((void **)&host_c, FULL_DATA_SIZE * sizeof(int), cudaHostAllocDefault);
//主机上的内存赋值
for (int i = 0; i < FULL_DATA_SIZE; i++) {
host_a[i] = i;
host_b[i] = FULL_DATA_SIZE - i;
}
//启动计时器
cudaEvent_t start, stop;
float elapsedTime;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
for (int i = 0; i < FULL_DATA_SIZE; i += 2 * N) {
cudaMemcpyAsync(dev_a, host_a + i, N * sizeof(int), cudaMemcpyHostToDevice, stream1);
cudaMemcpyAsync(dev_b, host_b + i, N * sizeof(int), cudaMemcpyHostToDevice, stream1);
cudaMemcpyAsync(dev_a1, host_a + i + N, N * sizeof(int), cudaMemcpyHostToDevice, stream2);
cudaMemcpyAsync(dev_b1, host_b + i + N, N * sizeof(int), cudaMemcpyHostToDevice, stream2);
kernel << < N / 1024, 1024, 0, stream1 >> > (dev_a, dev_b, dev_c);
kernel << < N / 1024, 1024, 0, stream2 >> > (dev_a1, dev_b1, dev_c1);
cudaMemcpyAsync(host_c + i, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost, stream1);
cudaMemcpyAsync(host_c + i + N, dev_c1, N * sizeof(int), cudaMemcpyDeviceToHost, stream2);
}
// 等待Stream流执行完成
cudaStreamSynchronize(stream1);
cudaStreamSynchronize(stream2);
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&elapsedTime, start, stop);
std::cout << "消耗时间： " << elapsedTime << std::endl;
//输出前10个结果
for (int i = 0; i < 10; i++) {
std::cout << host_c[i] << std::endl;
}
getchar();
// free stream and mem
cudaFreeHost(host_a);
cudaFreeHost(host_b);
cudaFreeHost(host_c);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);
cudaFree(dev_a1);
cudaFree(dev_b1);
cudaFree(dev_c1);
cudaStreamDestroy(stream1);
cudaStreamDestroy(stream2);
return 0;
}
首先声明一个Stream，可以把不同的操作放到Stream内，按照放入的先后顺序执行。
cudaMemcpyAsync操作只是一个请求，表示在流中执行一次内存复制操作，并不能确保cudaMemcpyAsync函数返回时已经启动了复制动作，更不能确定复制操作是否已经执行完成，可以确定的是放入流中的这个复制动作一定是在其后放入流中的其他动作之前完成的。使用一个流（同时要使用页锁定内存）和使用两个流的结果一致，运算时间分别是87ms和61ms。不使用流的同步方式，其执行时间最长，为109ms.
#include "cuda_runtime.h"
#include <iostream>
#include <stdio.h>
#include <math.h>
#define N (1024*1024)
#define FULL_DATA_SIZE N*20
__global__ void kernel(int *a, int *b, int *c)
{
int threadID = blockIdx.x * blockDim.x + threadIdx.x;
if (threadID < N) {
c[threadID] = (a[threadID] + b[threadID]) / 2;
}
}
int main()
{
int *host_a, *host_b, *host_c;
int *dev_a, *dev_b, *dev_c;
//在GPU上分配内存
cudaMalloc((void **)&dev_a, FULL_DATA_SIZE * sizeof(int));
cudaMalloc((void **)&dev_b, FULL_DATA_SIZE * sizeof(int));
cudaMalloc((void **)&dev_c, FULL_DATA_SIZE * sizeof(int));
//在CPU上分配可分页内存
host_a = (int *)malloc(FULL_DATA_SIZE * sizeof(int));
host_b = (int *)malloc(FULL_DATA_SIZE * sizeof(int));
host_c = (int *)malloc(FULL_DATA_SIZE * sizeof(int));
//主机上的内存赋值
for (int i = 0; i < FULL_DATA_SIZE; i++) {
host_a[i] = i;
host_b[i] = FULL_DATA_SIZE - i;
}
//启动计时器
cudaEvent_t start, stop;
float elapsedTime;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
//从主机到设备复制数据
cudaMemcpy(dev_a, host_a, FULL_DATA_SIZE * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, host_b, FULL_DATA_SIZE * sizeof(int), cudaMemcpyHostToDevice);
kernel <<< FULL_DATA_SIZE / 1024, 1024 >>> (dev_a, dev_b, dev_c);
//数据拷贝回主机
cudaMemcpy(host_c, dev_c, FULL_DATA_SIZE * sizeof(int), cudaMemcpyDeviceToHost);
//计时结束
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
cudaEventElapsedTime(&elapsedTime, start, stop);
std::cout << "消耗时间： " << elapsedTime << std::endl;
//输出前10个结果
for (int i = 0; i < 10; i++) {
std::cout << host_c[i] << std::endl;
}
getchar();
cudaFreeHost(host_a);
cudaFreeHost(host_b);
cudaFreeHost(host_c);
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);
return 0;
}
cuda-gdb调试kernel代码
make debug得到带调试信息的可执行文件，接着用cuda-gdb进行调试：
all:
nvcc -maxrregcount=16 study1.cu
debug:
nvcc -g -G study1.cu
dis:
cuobjdump -sass a.out > disas.s
cuobjdump -ptx a.out > ptxdis.s
cubin:
nvcc --cubin study1.cu
nvdisasm -cfg study1.cubin |dot -ocfg.png -Tpng
fatbin:
nvcc --fatbin study1.cucuda-gdb的用法和gdb类似，gdb上常用的调试方式和调试命令，在cuda-gdb上也能找到对应的影子。
和CPU上不同，GPU以组为单位调度线程，一组内可能有多个线程，这些线程被称为warp.每个warp的运行上下文包括寄存器和PC指针.GPU在一个Block内可以调度多个warp运行，直到把所有的寄存器用完。
我们在kernel入口处打断点，通过info cuda warps命令，可以看到当前cuda warp的调度信息。简单分析一下，发现最右边的threadIdx的第一个维度（因为代码中只定义了1个维度）的开始序号0，32，...160，...等都是32的倍数，这是因为这款GPU的warp内的线程数量是32.通过cuda-gdb 从侧面印证了这一点。
2023/11/30: MX250有3个SM，每个SM同一时刻只能调度一个BLOCK，程序中定义了32个BLOCK，所以可以看到当前SM0调度了BLOCK27，每个BLOCK有256个线程，所以调度到了BLOCK27上的线程160-191对应的WARP. 每个SM上可以同时调度多个BLOCK并发执行。
2024/05/03:block内的线程要整体调度，不能跨SM调度。
注意到active lanes mask和divergent lanes mask表示warp内出现分支的情况，由于我们在kernel函数中针对tid==0的情况做了判断，用来进行记时等操作，所以只有tid==0对应的lane mask在两个mask 中有所区别。
查看设备信息
设备信息如下：
设备号为0，PCIE BDF:02:00.0. 设备有3个SM，每个SM最大支持调度64个WARP，每个WARP用到了32LANE，每个LAN有256个寄存器。 active sms mask为0x7，可能的意思是意味着有3个SMS设备。
设置焦点
当焦点设置到主机线程时，命令将仅适用于该主机线程. 在设备端，焦点总是设置在最低粒度级别——设备线程。
(cuda-gdb) cuda device 0 sm 1
切换到指定的WARP
(cuda-gdb) cuda device 0 sm 2 warp 2 lane 2
由于切换到了SM2上的WARP2的第二个LANE上运行，查看对应WARP上调度的是BLOCK2上的WARP2，由于WARP2上的线程从THREAD64开始，所以LANE2对应线程66，如下图所示：
如果其他的忽略没写的话， 则认为其他变量是在现有的位置之下的变量。
定位当前位置
定位到具体的KERNEL，BLOCK，THREAD，SM，DEVICE等维度。
cuda device sm warp lane block thread
cuda kernel block thread
将KERNEL设置为64Block，32Threads/Block:
此时64个BLOCK分布在三个SM上，根据MASK显示，每个SM上分别有22，21，21个WARP，22+21+21=64个WARP。在BLOCK为32线程情况下，每个BLOCK对应一个WARP。
BLOCK按照顺序平均洒在三个SM上
SM0:
SM1:
SM2:
24 blocks, 256threads/block
24个BLOCK，每个BLOCK 256线程，也就是8个WARP，总共192个WARP，每个SM最大调度能力是64个WARP，正好能打满3个SM。
查看每个SM的调度情况，可以看到64个WARP全部打满：
SM0:
SM1:
SM2:
如果BLOCK数量超过了同时并发执行的BLOCK容量，比如这里24个BLOCK下刚好打满3个SM，如果BLOCK数量设置为32，THREAD_NUM仍然为8会如何呢？
经过测试，这种情况下会首先跑0-23个BLOCK，仍然会将三个SM打满，后面剩下的8个BLOCK会排队等待，当前面24个BLOCK执行完成后，再调度后面8个BLOCK执行。直到完成所有的32个BLOCK。
所以，如果BLOCK中的WARP超过了SM的容量，BLOCK在SM上的执行是排队的。并且以BLOCK为单位整体调度进入SM执行。
48 blocks, 128 threads/block
仍然可以打满三个SM：
96 blocks, 64 threads/block
也可以打满三个SM
192 blocks, 32 threads/block
此时无法打满三个SM了，所以，根据以上的测试可以推论出，如果要发挥GPU的最大性能，每个BLOCK中最好包含多个WARP，这个是软件可以控制的，软件无法控制的是SM的最大调度能力。
如果BLOCK中设置的线程数量过大，会如何？
前面devicequery获取到的信息，每个BLOCK的最大容量为1024线程，每个SM最大的调度能力为64WARP，也就是2048线程。也就是说，如果BLOCK 维度为1024的话，每个SM最多可以调度两个BLOCK，如下图所示，SM0上调度了BLOCK0 和BLOCK 3。
并且，BLOCK0和BLOCK3既然显示被调度，则一定是BLOCK内的所有线程都被调度了，不存在只调度BLOCK中的部分线程的情况，所以定义BLOCK的形状时要注意，不能超过硬件的限制。
单线程调试
#include <cuda_runtime.h>
#include <stdio.h>
__global__ void add(int a, int b, int *c)
{
*c = a + b;
}
int main(void)
{
int c;
int *gpu_c;
cudaMalloc((void **)&gpu_c, sizeof(int));
add<<<1,1>>>(3,6,gpu_c);
cudaMemcpy(&c, gpu_c, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(gpu_c);
printf("3 + 9 eques %d.\n", c);
return 0;
}
因为只发起了一个线程，所以SM只启动了一个LANE,其余31个LANE都没有启用。
block维数设置为32:
此时32LANE全部活跃,唯一的一个BLOCK内的32个线程分布在SM的32个LANE上。
LANE寄存器调试
运行kernel的时候执行 info register命令，查看到当前LANE上有256个R寄存器，一个PC寄存器，64个UR寄存器和8个UP寄存器：
R寄存器每个LANE独立的还是共享同一份?
在LANE0中设置R0为0xdeadbeef,切换到LANE1上去查看，如果R0不同值，说明每个LANE有独立的R寄存器：
按照上图的调试序列，发现R寄存器每个LANE有独一份。
CUDA线程块内同步机制：
在CUDA编程模型中，__syncthreads();可以用于同一线程块内线程的同步操作，它对应的PTX指令为bar指令，该指令会在其所在的程序计数器（PC）位置产生一个同步栅栏(barrier),并要求线程块内所有的线程都达到这一栅栏位置才能继续执行，这可以通过监控线程的PC来实现，如下图所示，在线程同步的要求下，即便有一些线程束运行比较快而先到达BAR指令处也要暂停，直到整个线程块的所有线程都达到BAR指令才能整体继续执行。__syncthreads函数保证了先于该语句的执行结果对线程块内所有的线程可见，由于线程块内部按照线程束组织，而且同一个线程束内的线程是自然同步的，所以__syncthreads实际上保证了线程块内所有的线程束都执行到统一位置再继续。如下图所示（图示有一个问题容易引起误解，图中横线代表的BAR指令穿越了两个线程块，BAR用于线程块内部的同步，不能跨线程块同步，所以直线在两个BLOCK中间的部分应该抹掉，另外，线程束内部的线程是lock-step锁步执行的，应该是齐头并进，并没有先后关系）。
测试，在KERNEL中添加对__syncthreads的调用，之后编译并反编译为PTX文件
PTX文件中出现了BAR指令：
BR支持同样的指令，并且将其重新宏定义为一个新的符号，宏名更加名副其，__sync_block_threads.
编译后，对应的汇编指令为bar指令
如果遇到计算单元STUCK住的问题，比如SPC或者SM一直卡着不执行，很有可能卡在了BAR指令这里，因为GPU线程执行的逻辑都非常简单，基本上都是顺序指令的依赖，如果没有其它访存之类的异常（这类异常通常以其它方式上报），基本上很少会有机会卡住执行流，BAR指令是最有可能的一个情况。
GPU Devices info
MX250 include 3 SMs, 64warps each SM, 32 thread each warp. you can get this information from cuda gdb by this command.
detail infomation about SMs
the mask corrspond bit set to "1" means the warp is active.
warp按照SM独立划分：
我的笔记本显卡型号是MX250，包含三个SM，按照上面cuda-gdb的输出信息，三个SM总共有32x64x3=6144个线程。
2023/11/30:这里的 64 Warp/SM可能并不是说每个SM可以同时并行运行64个WARP，而是指能并发执行调度这么多。就和SUPA中1个EU可以调度8个WARP一样（1个CU就是32个WARP），也不是并行执行，而是并发调度这么多。所以上面的图GDB调试时，只有一个WARP有*号表示正在运行的WARP，也就是硬件资源只够每个SM同一时刻执行一个WARP，GPU有三个SM，也就是GPU同一时刻最多有三个WARP在执行，是硬件的分时复用造成了一个SM可以运行64个WARP的假象。
而按照上面的例子对计算网格的划分，包含32个256个thread的BLOCK，也就是有32x256=8192个线程，显然即便动员三个SM上的所有线程，也不满足。
怎么办呢？解决办法仍然是分时复用，3个SM同时刻只能调度进24个block(24x256 = 6144).
2023/12/17: 根据cuda deviceQuery的信息，每个SM最多支持2048个线程，对应64个warp,每个SM有128个CUDA CORE，三个SM有384个CUDA CORE。而每个SM貌似有32 LANE，是不是说每LANE有4个CUDA CORE？
2024/05/03:每个CUDA Core是scaler单元，1个LANE四个Cuda core表示1个LANE有四路SIMD。1SM有32个LANE就是有128个 Cuda core.每个CUDA Core上刷16个线程，就是2048个线程，32个线程为一个WARP，就是64WARP/SM。
2024/06/25:MX250 有384个cuda core(SP)也是百度数据，并不一定可靠，主要原因是CUDA CORE是后面演进出的概念，可能和LANE是对等的一个概念，并不能用在一起。
可能是调度原因，前面的调试没有看到完整的64个WARP，下面则抓到一次：
组合使用如下cuda-gdb命令可以查看不同的warp&block在三个sm上的分布情况
$ info cuda warps
$ cuda block(0,0,0) thread(0,0,0)
$ cuda block(1,0,0) thread(0,0,0)
$ cuda block(2,0,0) thread(0,0,0)
$ cuda block(21,0,0) thread(0,0,0)
$ cuda block(22,0,0) thread(0,0,0)
$ cuda block(23,0,0) thread(0,0,0)一开始运行的时候，全部24个block投入运行将会占满整个GPU的算力，后续8个block还没有投入运行，此时是无法通过cuda-gdb调试这些线程的。比如，当使用cuda block(28,0,0) thread(0,0,0)切换线程时，会发现线程不存在。
当有warp完成后退出，没有执行的block才允许被调度进来，我们用条件断点来观察：
break 127 if blockIdx.x==28
lane的复用
下面两张图分别为block (2,0,0)下的WARP2和WARP6 两个WARP的执行现场，可以看到他们复用了SM2的32个LANE，也就是说LANE是分时复用的。
warp和SIMD
现代处理器一般包括cpu core和SIMD.在统计核数的时候，SIMD的通道数是不考虑在内的。与GPU的线程数对比，CPU的这种统计方式有些吃亏，因为如果按照GPU的统计方式，GPU的核数至少需要乘以SIMD的通道数。而如果按照CPU的统计方式，GPU的核数需要除以warp数。
打一个不恰当的比喻，CPU统计的核数是相当于统计有多少条公路，而GPU统计的核数相当于统计公路的车道数。
计算单元的层次结构
cuda-gdb对程序的影响
当使用gdb运行程序时，所有的功能都会变成同步的，所以kernel运行时间会明显边长。以此为例，增加-g -G选项编译后，执行时间从先前的8ms骤然增加到70多ms.
zl@czl-RedmiBook-14:~/workspace/cudastudy/study1$ make
nvcc -maxrregcount=16 study1.cu
study1.cu(23): warning #177-D: function "sum_of_squares" was declared but never referenced
study1.cu(50): warning #177-D: function "sum_of_squares_eff" was declared but never referenced
study1.cu(87): warning #177-D: function "sum_of_squares_eff_2" was declared but never referenced
czl@czl-RedmiBook-14:~/workspace/cudastudy/study1$ ./a.out
there are 1 cudda device found, id 0.
===============================Device 0==================================
NVIDIA GeForce MX250
Total global memory: 2099970048 Bytes
Max shareable memory per block 49152 Bytes.
Maximum registers per block: 65536
Wrap Size 32.
Maximum threads per block 1024.
Maximum block dimensions [1024, 1024, 64].
Maximum grid dimensions [2147483647, 65535, 65535].
Total constant memory: 65536.
Support compute Capability: 6.1.
Kernel Frequency 1582000 kHz.
Number of MultProcessors 3.
Is MultiGPU: No.
L2 Cache Size: 524288 Bytes.
Memory Bus Width: 64.
ECC status: Disable.
=========================================================================
GPU Elapsed time:8.153568 ms.
main line 319 sum 104857600, time 8 ms.
czl@czl-RedmiBook-14:~/workspace/cudastudy/study1$ make debug
nvcc -g -G study1.cu
study1.cu(23): warning #177-D: function "sum_of_squares" was declared but never referenced
study1.cu(50): warning #177-D: function "sum_of_squares_eff" was declared but never referenced
study1.cu(87): warning #177-D: function "sum_of_squares_eff_2" was declared but never referenced
czl@czl-RedmiBook-14:~/workspace/cudastudy/study1$ ./a.out
there are 1 cudda device found, id 0.
===============================Device 0==================================
NVIDIA GeForce MX250
Total global memory: 2099970048 Bytes
Max shareable memory per block 49152 Bytes.
Maximum registers per block: 65536
Wrap Size 32.
Maximum threads per block 1024.
Maximum block dimensions [1024, 1024, 64].
Maximum grid dimensions [2147483647, 65535, 65535].
Total constant memory: 65536.
Support compute Capability: 6.1.
Kernel Frequency 1582000 kHz.
Number of MultProcessors 3.
Is MultiGPU: No.
L2 Cache Size: 524288 Bytes.
Memory Bus Width: 64.
ECC status: Disable.
=========================================================================
GPU Elapsed time:66.803780 ms.
main line 319 sum 104857600, time 72 ms.
czl@czl-RedmiBook-14:~/workspace/cudastudy/study1$
反编译：
执行make dis进行反编译，可以反编译fabin段，得到NVIDIA GPU的二进制指令码，可以看到，类似于CPU，NVIDIA GPU也有自己的指令集，每条指令都是8个字节长度，这一点也非常类似于RISC处理器指令长度统一的规整模式。从这一点看，GPU处理器没有什么特别的。
$ cuobjdump -sass a.out > disas.s
$ cuobjdump -ptx a.out > ptxdis.s
生成流程图
nvdisasm 工具可以生成流程描述文件，再通过dog工具将描述文件转换流程图，执行make cubin命令：
$ nvcc --cubin study1.cu
$ nvdisasm -cfg study1.cubin |dot -ocfg.png -Tpng
texture
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#define DATA_SIZE 8
texture<int, 1, cudaReadModeElementType>reftex;
cudaArray *cuarray;
void init_cuda(void)
{
int count, dev;
int i;
cudaDeviceProp prop;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "there is no cuda device.\n");
return;
} else {
cudaGetDevice(&dev);
fprintf(stdout,"there are %d cudda device found, id %d.\n", count, dev);
}
for(i = 0; i < count; i ++) {
printf("===============================Device %d==================================\n", i);
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
printf("%s\n", prop.name);
printf("Total global memory: %ld Bytes\n", prop.totalGlobalMem);
printf("Max shareable memory per block %ld Bytes.\n", prop.sharedMemPerBlock);
printf("Maximum registers per block: %d\n", prop.regsPerBlock);
printf("Wrap Size %d.\n", prop.warpSize);
printf("Maximum threads per block %d.\n", prop.maxThreadsPerBlock);
printf("Maximum block dimensions [%d, %d, %d].\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("Maximum grid dimensions [%d, %d, %d].\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("Total constant memory: %ld.\n", prop.totalConstMem);
printf("Support compute Capability: %d.%d.\n", prop.major, prop.minor);
printf("Kernel Frequency %d kHz.\n", prop.clockRate);
printf("Number of MultProcessors %d.\n", prop.multiProcessorCount);
printf("Is MultiGPU: %s.\n", prop.isMultiGpuBoard ? "Yes" : "No");
printf("L2 Cache Size: %d Bytes.\n", prop.l2CacheSize);
printf("Memory Bus Width: %d.\n", prop.memoryBusWidth);
printf("ECC status: %s.\n", prop.ECCEnabled? "Enable" : "Disable");
}
printf("=========================================================================\n");
}
cudaSetDevice(1);
}
__global__ void convolution(unsigned unsizedata, int *pres)
{
const int idx = threadIdx.x;
pres[idx] = unsizedata + 1 - tex1D(reftex, idx);
}
int main(void)
{
unsigned unsizedata = DATA_SIZE;
unsigned data = 0;
int *sampler;
int i;
init_cuda();
cudaMallocHost((void**)&sampler, unsizedata * sizeof(int));
if(sampler == NULL){
fprintf(stderr, "malloc host buffer failure.\n");
return -1;
}
for(i = 0; i < unsizedata; i ++) sampler[i] = ++data;
for(i = 0; i < unsizedata; i ++) {
printf("%d\t", sampler[i]);
}
printf("\n");
// prepare texture
cudaChannelFormatDesc cuDesc = cudaCreateChannelDesc<int>();
cudaMallocArray(&cuarray, &cuDesc, unsizedata);
cudaMemcpyToArray(cuarray, 0, 0, sampler, unsizedata * sizeof(int), cudaMemcpyHostToDevice);
cudaBindTextureToArray(reftex, cuarray);
int *pres;
cudaMalloc((void**)&pres, unsizedata * sizeof(int));
convolution <<<1, unsizedata >>>(unsizedata, pres);
cudaMemcpy(sampler, pres, unsizedata * sizeof(int), cudaMemcpyDeviceToHost);
for(i = 0; i < unsizedata; i ++) {
printf("%d\t", sampler[i]);
}
printf("\n");
cudaUnbindTexture(reftex);
cudaFreeHost(sampler);
cudaFreeArray(cuarray);
cudaFree(pres);
return 0;
}flow:
验证warp数量，设置为1 block 8线程处理8笔数据
1 block 32 线程, and 32笔数据.
matrix Add
下图实现两个8X16的二维矩阵加法，按照warp大小对计算进行切分，示意图如下：
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#define MATRIX_M 8
#define MATRIX_N 16
void init_cuda(void)
{
int count, dev;
int i;
cudaDeviceProp prop;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "there is no cuda device.\n");
return;
} else {
cudaGetDevice(&dev);
fprintf(stdout,"there are %d cudda device found, id %d.\n", count, dev);
}
for(i = 0; i < count; i ++) {
printf("===============================Device %d==================================\n", i);
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
printf("%s\n", prop.name);
printf("Total global memory: %ld Bytes\n", prop.totalGlobalMem);
printf("Max shareable memory per block %ld Bytes.\n", prop.sharedMemPerBlock);
printf("Maximum registers per block: %d\n", prop.regsPerBlock);
printf("Wrap Size %d.\n", prop.warpSize);
printf("Maximum threads per block %d.\n", prop.maxThreadsPerBlock);
printf("Maximum block dimensions [%d, %d, %d].\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("Maximum grid dimensions [%d, %d, %d].\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("Total constant memory: %ld.\n", prop.totalConstMem);
printf("Support compute Capability: %d.%d.\n", prop.major, prop.minor);
printf("Kernel Frequency %d kHz.\n", prop.clockRate);
printf("Number of MultProcessors %d.\n", prop.multiProcessorCount);
printf("Is MultiGPU: %s.\n", prop.isMultiGpuBoard ? "Yes" : "No");
printf("L2 Cache Size: %d Bytes.\n", prop.l2CacheSize);
printf("Memory Bus Width: %d.\n", prop.memoryBusWidth);
printf("ECC status: %s.\n", prop.ECCEnabled? "Enable" : "Disable");
}
printf("=========================================================================\n");
}
cudaSetDevice(1);
}
void init_matrix(int *data, int m, int n, int val)
{
int i = 0, j = 0;
for(i = 0; i < m; i ++){
for(j = 0; j < n; j ++){
*(data + i * MATRIX_N +j) = val;
}
}
}
void print_matrix(int *data, int m, int n)
{
int i, j;
for(i = 0; i < m; i ++){
printf("\n");
for(j = 0; j < n; j ++){
printf("%d\t", *(data + i * MATRIX_N +j));
}
}
printf("\n");
}
__global__ static void sum_of_matrix(int *matrixA, int *matrixB, int *matrixC)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
printf("%s line %d,i = %d, j = %d.\n", __func__, __LINE__, i, j);
if(i < MATRIX_N && j < MATRIX_M){
*(matrixC + j * MATRIX_N + i) = *(matrixA + j * MATRIX_N + i) + *(matrixB + j * MATRIX_N + i);
}
}
int main(void)
{
int *pdataA;
int *pdataB;
int *pdataC;
int *pdata_gpuA;
int *pdata_gpuB;
int *pdata_gpuC;
init_cuda();
cudaMallocHost((void**)&pdataA, MATRIX_M * MATRIX_N * sizeof(int));
cudaMallocHost((void**)&pdataB, MATRIX_M * MATRIX_N * sizeof(int));
cudaMallocHost((void**)&pdataC, MATRIX_M * MATRIX_N * sizeof(int));
if(pdataA == NULL || pdataB == NULL || pdataC == NULL) {
fprintf(stderr, "fatal error, malloc host buffer failure.\n");
return -1;
}
cudaHostGetDevicePointer((void**)&pdata_gpuA, (void*)pdataA, 0);
cudaHostGetDevicePointer((void**)&pdata_gpuB, (void*)pdataB, 0);
cudaHostGetDevicePointer((void**)&pdata_gpuC, (void*)pdataC, 0);
if(pdata_gpuA == NULL || pdata_gpuB == NULL || pdata_gpuC == NULL) {
fprintf(stderr, "fatal error, get device ptr failure.\n");
return -1;
}
init_matrix(pdataA, MATRIX_M, MATRIX_N, 1);
init_matrix(pdataB, MATRIX_M, MATRIX_N, 2);
init_matrix(pdataC, MATRIX_M, MATRIX_N, 0);
print_matrix(pdataA, MATRIX_M, MATRIX_N);
print_matrix(pdataB, MATRIX_M, MATRIX_N);
print_matrix(pdataC, MATRIX_M, MATRIX_N);
dim3 bloksize(8, 4, 1);
dim3 gridsize(2, 2, 1);
sum_of_matrix<<< gridsize, bloksize >>>(pdata_gpuA, pdata_gpuB, pdata_gpuC);
cudaDeviceSynchronize();
print_matrix(pdataA, MATRIX_M, MATRIX_N);
print_matrix(pdataB, MATRIX_M, MATRIX_N);
print_matrix(pdataC, MATRIX_M, MATRIX_N);
cudaFree(pdataA);
cudaFree(pdataB);
cudaFree(pdataC);
return 0;
}
此时四个BLOCK在三个SM上的分布情况如下图所示,BLOCK是横向分配，最大限度利用GPU上的所有SM，并且由于BLOCK不能跨多个SM，所以你会看到如下的BLOCK分布。
同时每个BLOCK有32个线程，对应一个WARP，恰好也有32个LANE，每个LANE对应一个线程：
也就是说，这个计算任务为计算两个8x16矩阵的加法，矩阵中有128个元素，正好对应了划分中的4个BLOCK，128个线程，每个线程计算一个元素的加法就可以了,如下图：
实际的线程块和线程ID分布（blockid << 16 | threadIdx.x):
在矩阵中的位置分布：
指令流程如下图所示：
matrix multiply
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#define MATRIX_M 16
#define MATRIX_N 16
void init_cuda(void)
{
int count, dev;
int i;
cudaDeviceProp prop;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "there is no cuda device.\n");
return;
} else {
cudaGetDevice(&dev);
fprintf(stdout,"there are %d cudda device found, id %d.\n", count, dev);
}
for(i = 0; i < count; i ++) {
printf("===============================Device %d==================================\n", i);
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
printf("%s\n", prop.name);
printf("Total global memory: %ld Bytes\n", prop.totalGlobalMem);
printf("Max shareable memory per block %ld Bytes.\n", prop.sharedMemPerBlock);
printf("Maximum registers per block: %d\n", prop.regsPerBlock);
printf("Wrap Size %d.\n", prop.warpSize);
printf("Maximum threads per block %d.\n", prop.maxThreadsPerBlock);
printf("Maximum block dimensions [%d, %d, %d].\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("Maximum grid dimensions [%d, %d, %d].\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("Total constant memory: %ld.\n", prop.totalConstMem);
printf("Support compute Capability: %d.%d.\n", prop.major, prop.minor);
printf("Kernel Frequency %d kHz.\n", prop.clockRate);
printf("Number of MultProcessors %d.\n", prop.multiProcessorCount);
printf("Is MultiGPU: %s.\n", prop.isMultiGpuBoard ? "Yes" : "No");
printf("L2 Cache Size: %d Bytes.\n", prop.l2CacheSize);
printf("Memory Bus Width: %d.\n", prop.memoryBusWidth);
printf("ECC status: %s.\n", prop.ECCEnabled? "Enable" : "Disable");
}
printf("=========================================================================\n");
}
cudaSetDevice(1);
}
void init_matrix(int *data, int m, int n, int val)
{
int i = 0, j = 0;
for(i = 0; i < m; i ++){
for(j = 0; j < n; j ++){
*(data + i * MATRIX_N +j) = val;
}
}
}
void print_matrix(int *data, int m, int n)
{
int i, j;
for(i = 0; i < m; i ++){
printf("\n");
for(j = 0; j < n; j ++){
printf("%d\t", *(data + i * MATRIX_N +j));
}
}
printf("\n");
}
__global__ static void sum_of_matrix(int *matrixA, int *matrixB, int *matrixC)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
printf("%s line %d,i = %d, j = %d.\n", __func__, __LINE__, i, j);
if(i < MATRIX_N && j < MATRIX_M){
*(matrixC + j * MATRIX_N + i) = *(matrixA + j * MATRIX_N + i) + *(matrixB + j * MATRIX_N + i);
}
}
__global__ static void multiply_of_matrix(int *matrixA, int *matrixB, int *matrixC)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
int sum = 0;
int cnt = 0;
if(i < MATRIX_N && j < MATRIX_M){
for(cnt = 0; cnt < 16; cnt ++) {
sum += *(matrixA + i * MATRIX_N + cnt) * (*(matrixB + cnt * MATRIX_N + j));
}
*(matrixC + j * MATRIX_N + i) = sum;
}
}
int main(void)
{
int *pdataA;
int *pdataB;
int *pdataC;
int *pdata_gpuA;
int *pdata_gpuB;
int *pdata_gpuC;
init_cuda();
cudaMallocHost((void**)&pdataA, MATRIX_M * MATRIX_N * sizeof(int));
cudaMallocHost((void**)&pdataB, MATRIX_M * MATRIX_N * sizeof(int));
cudaMallocHost((void**)&pdataC, MATRIX_M * MATRIX_N * sizeof(int));
if(pdataA == NULL || pdataB == NULL || pdataC == NULL) {
fprintf(stderr, "fatal error, malloc host buffer failure.\n");
return -1;
}
cudaHostGetDevicePointer((void**)&pdata_gpuA, (void*)pdataA, 0);
cudaHostGetDevicePointer((void**)&pdata_gpuB, (void*)pdataB, 0);
cudaHostGetDevicePointer((void**)&pdata_gpuC, (void*)pdataC, 0);
if(pdata_gpuA == NULL || pdata_gpuB == NULL || pdata_gpuC == NULL) {
fprintf(stderr, "fatal error, get device ptr failure.\n");
return -1;
}
init_matrix(pdataA, MATRIX_M, MATRIX_N, 1);
init_matrix(pdataB, MATRIX_M, MATRIX_N, 2);
init_matrix(pdataC, MATRIX_M, MATRIX_N, 0);
print_matrix(pdataA, MATRIX_M, MATRIX_N);
print_matrix(pdataB, MATRIX_M, MATRIX_N);
print_matrix(pdataC, MATRIX_M, MATRIX_N);
dim3 bloksize(8, 4, 1);
dim3 gridsize(2, 4, 1);
//sum_of_matrix<<< gridsize, bloksize >>>(pdata_gpuA, pdata_gpuB, pdata_gpuC);
multiply_of_matrix<<< gridsize, bloksize >>>(pdata_gpuA, pdata_gpuB, pdata_gpuC);
cudaDeviceSynchronize();
print_matrix(pdataA, MATRIX_M, MATRIX_N);
print_matrix(pdataB, MATRIX_M, MATRIX_N);
print_matrix(pdataC, MATRIX_M, MATRIX_N);
cudaFree(pdataA);
cudaFree(pdataB);
cudaFree(pdataC);
return 0;
}
计算结构，二维的线程矩阵与二维的结果矩阵元素一一对应，每个线程计算matrixC矩阵中的一个元素。
矩阵乘法的例子block在SM上的分配图
分块矩阵乘法实现：
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#define MATRIX_M 16
#define MATRIX_N 16
void init_cuda(void)
{
int count, dev;
int i;
cudaDeviceProp prop;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "there is no cuda device.\n");
return;
} else {
cudaGetDevice(&dev);
fprintf(stdout,"there are %d cudda device found, id %d.\n", count, dev);
}
for(i = 0; i < count; i ++) {
printf("===============================Device %d==================================\n", i);
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
printf("%s\n", prop.name);
printf("Total global memory: %ld Bytes\n", prop.totalGlobalMem);
printf("Max shareable memory per block %ld Bytes.\n", prop.sharedMemPerBlock);
printf("Maximum registers per block: %d\n", prop.regsPerBlock);
printf("Wrap Size %d.\n", prop.warpSize);
printf("Maximum threads per block %d.\n", prop.maxThreadsPerBlock);
printf("Maximum block dimensions [%d, %d, %d].\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("Maximum grid dimensions [%d, %d, %d].\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("Total constant memory: %ld.\n", prop.totalConstMem);
printf("Support compute Capability: %d.%d.\n", prop.major, prop.minor);
printf("Kernel Frequency %d kHz.\n", prop.clockRate);
printf("Number of MultProcessors %d.\n", prop.multiProcessorCount);
printf("Is MultiGPU: %s.\n", prop.isMultiGpuBoard ? "Yes" : "No");
printf("L2 Cache Size: %d Bytes.\n", prop.l2CacheSize);
printf("Memory Bus Width: %d.\n", prop.memoryBusWidth);
printf("ECC status: %s.\n", prop.ECCEnabled? "Enable" : "Disable");
}
printf("=========================================================================\n");
}
cudaSetDevice(1);
}
void init_matrix(int *data, int m, int n, int val)
{
int i = 0, j = 0;
for(i = 0; i < m; i ++){
for(j = 0; j < n; j ++){
*(data + i * MATRIX_N +j) = val;
}
}
}
void print_matrix(int *data, int m, int n)
{
int i, j;
for(i = 0; i < m; i ++){
printf("\n");
for(j = 0; j < n; j ++){
printf("%d\t", *(data + i * MATRIX_N +j));
}
}
printf("\n");
}
__global__ static void sum_of_matrix(int *matrixA, int *matrixB, int *matrixC)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
printf("%s line %d,i = %d, j = %d.\n", __func__, __LINE__, i, j);
if(i < MATRIX_N && j < MATRIX_M){
*(matrixC + j * MATRIX_N + i) = *(matrixA + j * MATRIX_N + i) + *(matrixB + j * MATRIX_N + i);
}
}
__global__ static void multiply_of_matrix(int *matrixA, int *matrixB, int *matrixC)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
int sum = 0;
int cnt = 0;
if(i < MATRIX_N && j < MATRIX_M){
for(cnt = 0; cnt < 16; cnt ++) {
//sum += *(matrixA + i * MATRIX_N + cnt) * (*(matrixB + cnt * MATRIX_N + j));
sum += matrixA[i * MATRIX_N + cnt] * matrixB[cnt * MATRIX_N + j];
}
//*(matrixC + j * MATRIX_N + i) = sum;
matrixC[j * MATRIX_N + i] = sum;
}
}
__global__ static void block_mul_matrix(int *matrixA, int *matrixB, int *matrixC)
{
#define BLOCK_SIZE 4
__shared__ int Mds[BLOCK_SIZE][BLOCK_SIZE];
__shared__ int Nds[BLOCK_SIZE][BLOCK_SIZE];
int col = threadIdx.x + blockIdx.x * blockDim.x;
int row = threadIdx.y + blockIdx.y * blockDim.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
int block_sum = 0;
for(int i = 0; i < 16/BLOCK_SIZE; i ++)
{
Mds[ty][tx] = matrixA[row * 16 + tx + i * BLOCK_SIZE];
Nds[ty][tx] = matrixB[col + (ty + i * BLOCK_SIZE) * 16];
__syncthreads();
for(int j = 0; j < BLOCK_SIZE; j ++)
{
block_sum += Mds[ty][j] * Nds[j][tx];
__syncthreads();
}
}
matrixC[row*MATRIX_N + col] = block_sum;
}
int main(void)
{
int *pdataA;
int *pdataB;
int *pdataC;
int *pdata_gpuA;
int *pdata_gpuB;
int *pdata_gpuC;
init_cuda();
cudaMallocHost((void**)&pdataA, MATRIX_M * MATRIX_N * sizeof(int));
cudaMallocHost((void**)&pdataB, MATRIX_M * MATRIX_N * sizeof(int));
cudaMallocHost((void**)&pdataC, MATRIX_M * MATRIX_N * sizeof(int));
if(pdataA == NULL || pdataB == NULL || pdataC == NULL) {
fprintf(stderr, "fatal error, malloc host buffer failure.\n");
return -1;
}
#if 1
cudaHostGetDevicePointer((void**)&pdata_gpuA, (void*)pdataA, 0);
cudaHostGetDevicePointer((void**)&pdata_gpuB, (void*)pdataB, 0);
cudaHostGetDevicePointer((void**)&pdata_gpuC, (void*)pdataC, 0);
#else
pdata_gpuA = pdataA;
pdata_gpuB = pdataB;
pdata_gpuC = pdataC;
#endif
if(pdata_gpuA == NULL || pdata_gpuB == NULL || pdata_gpuC == NULL) {
fprintf(stderr, "fatal error, get device ptr failure.\n");
return -1;
}
printf("%s line %d, pdataA = %p, pdata_gpuA = %p.\n", __func__, __LINE__, pdataA, pdata_gpuA);
if(pdataA != pdata_gpuA) {
printf("%s line %d, not uma address space.\n", __func__, __LINE__);
}
if(pdataB != pdata_gpuB) {
printf("%s line %d, not uma address space.\n", __func__, __LINE__);
}
if(pdataC != pdata_gpuC) {
printf("%s line %d, not uma address space.\n", __func__, __LINE__);
}
init_matrix(pdataA, MATRIX_M, MATRIX_N, 1);
init_matrix(pdataB, MATRIX_M, MATRIX_N, 2);
init_matrix(pdataC, MATRIX_M, MATRIX_N, 0);
//print_matrix(pdataA, MATRIX_M, MATRIX_N);
//print_matrix(pdataB, MATRIX_M, MATRIX_N);
//print_matrix(pdataC, MATRIX_M, MATRIX_N);
dim3 bloksize(4, 4, 1);
dim3 gridsize(4, 4, 1);
//sum_of_matrix<<< gridsize, bloksize >>>(pdata_gpuA, pdata_gpuB, pdata_gpuC);
//multiply_of_matrix<<< gridsize, bloksize >>>(pdata_gpuA, pdata_gpuB, pdata_gpuC);
block_mul_matrix<<< gridsize, bloksize >>>(pdata_gpuA, pdata_gpuB, pdata_gpuC);
cudaDeviceSynchronize();
//print_matrix(pdataA, MATRIX_M, MATRIX_N);
//print_matrix(pdataB, MATRIX_M, MATRIX_N);
print_matrix(pdataC, MATRIX_M, MATRIX_N);
cudaFree(pdataA);
cudaFree(pdataB);
cudaFree(pdataC);
return 0;
}其中，循环内的__syncthreads作用说明如下：
第一次WARP计算：
第二次WARP计算：
第三次WARP计算：
第四次WARP计算：
四此计算中，必须加入__syncthread同步。
run result:
From the code, we see pdata_gpuA Is equal pdataA, pdata_gpuB is equal pdataB, pdata_gpuC is equal C, so it necessary calling cudaHostGetDevicePointer? we comment it and see the result.
you can find it still fine,so it seems not necessary to call cudaHostGetDevicePointer, is this universal for all scenario?
下图是nvidia关于函数cudaHostGetDevicePointer的官方文档,没看到说这个函数可以不用的说明，可能这个函数的作用是pin memory,让memroy不swapout而是常驻内存，所以不影响功能，纯属猜测。
CUDA ISA架构
中间指令PTX：
在上面的makefile写法中，反编译阶段的-ptx指导反编译器产生ISA中间指令代码。它非常类似于一种ISA,但实际上又不是，更像是一层抽象的IR语言，PTX面向的是virtual GPU architecture。NVIDIA 历代GPGPU的底层架构不尽相同，运行的机器指令也存在一定差异，因此CUDA定义了一套稳定的，和具体架构无关的指令集以及变成系统，称为PTX（parallel thread execution). PTX与底层GPGPU硬件架构基本无关，能够跨越多种硬件，在运行时再转化为更底层的机器指令集SASS执行。从而更具通用性。
机器指令集SASS
-sass(Shader Assembly?)选项则指引反编译器产生二进制指令。如果PTX面向的是虚拟GPU架构，那么SASS面向的就是具体的某款GPU实现了，他是由PTX代码编译而成。对应不同的GPGPU底层架构。如下图中，最右边列的二进制编码就是SASS指令，当前每条指令都是64BIT，不同的架构下指令编码长度不同，比如Touring架构下每个指令是128BIT编码。第一列代表的是地址。
特殊寄存器
kernel函数中有几个关键的builtin变量，比如blockIdx,threadIdx,他们是绑定到架构中的特殊寄存器中的，如下图所示：
S2R和CS2R，把特殊register的值载入到GPR。常用的特殊寄存器有: SR_TID.X, SR_TID.Z, SR_TID.Z, SR_CTAID.X, SR_CTAID.Y, SR_CTAID.Z, SR_LANEID, SR_*MASK（LT，LE，GT，GE等）, SR_CLOCKHI/SR_CLOCKLO等等。具体可以参考PTX关于特殊寄存器的文档。CS2R与S2R的区别是CS2R的latency是固定的，S2R则不固定。不过CS2R好像就只看到支持SRZ（就是RZ？）和SR_CLOCKHI/SR_CLOCKLO。这里附带提一句，blockIdx对应SR_CTAID，threadIdx对应SR_TID，两者的xyz都有special register，那blockDim和gridDim呢？其实很简单，只有变动或是不确定的数才有放到special register内的价值，blockDim和gridDim在kernel运行时都是全局定值，现在都会被放在constant memory里。编译器会自动把相应访问转为对constant memory相应地址的访问。
构建方式
反编译
和GCC相比，NVCC不太智能的是，对于.cu文件中定义但是没有调用的kernel函数，NVCC在编译时不会自动优化掉，而是编译进最后的可执行文件中。如下图所示：
另外需要注意的是， 每个Kernel的地址都是从0开始的，这个不太符合HOST端连接器的对代码和数据的布局和分配机制，在HOST端，每个代码或者数据都是有唯一的运行地址的，不能冲突，程序加载时，也必须要加载到正确的运行地址，否则会跑飞。
但是从对于kernel却没问题，这是因为编译器产生的运行地址对DEVICE侧的运行里说毫无意义，远端也不看这个信息。Kernel实际的运行地址是临时分配的device memory,也就是DEVICE端是地址无关代码，拷贝到哪里都可以被GPU执行。
存储模型
从GPU的结构来看，是分成Grid->Block->Thread，所以，GPU也是针对这个层次来设计存储结构的，从上面的表格中可以看到。register和local memory是针对Thread来设计的，shared memory是针对Block来设计，而其他的三个(constant, texture, global)是针对Grid来进行设计的。
共享存储是BLOCK私有的，无法跨BLOCK共享，尽管如此，调试发现，在同一个KERNEL中，不同的BLOCK在运行时访问的共享存储地址是相同的,比如如下的KERNEL定义，其printed和shared是共享存储，在不同的BLOCK中其地址相同,但是对应的是不同的物理存储单元：
当私有本地寄存器（TLR，Thead Local Register)不足时，可以从device 存储区申请一块TLM（Thread Local Memory)，用TLM当做TLR用，存储kernel里的临时变量,TLR和TLM分别对应下表中的私有本地存储器和本地存储器。
| CUDA | Scope | Life Time | Managed By | 实现设备 |
| TLM(Local Memory) | Thread | Thread | Application | 显存 |
| L1 | all threads in the block | thread | hardware | SRAM |
| Share Memory | all threads in the block | Thread block | Application | 显存 |
| L2 | All threads | Application | Application | SRAM |
| Global Memory | All threads | Application | Application | 显存 |
| Constant Memory | All threads | Application | Application | 显存 |
TLM,也就是本地存储器，之所以称之为本地存储器，是因为它的作用域同私有本地存储器，也就是寄存器一样，只能被单个线程访问，这个名称不意味着它能提供高速的数据传输速度，实际上，本地存储器位于设备存储器上，因此访问本地存储器消耗的时间同访问全局存储器是一样的。本地存储器唯一的用途即为存放部分自动变量，自动变量是那种在全局函数和设备函数中不适用限定符说明的变量，自动变量由系统决定它们放在哪里。一般情况下，自动变量都放在寄存器TLR中，当线程使用的寄存器大小超过了实际的寄存器容量时，会发生溢出，此时，溢出的数据会被暂时写入到TLM。可能会被存放在TLM里的自动变量常为较大的结构体或者数组，以及被编译器认为具有动态下标的数组,如下图，对于kernel声明的巨大数组，TLR占不下，nvcc编译器将其置于TLM中。
cuobjdump -ptx a.out 反编译为PTX汇编，查看PTX中kernel中的local描述,确实被放到了TLM中。
通过nvcc matrix.cu --ptxas-options=-v 命令，还可以得到每个kernel的寄存器以及TLM分配细节信息，如下图：
zlcao@zlcao-RedmiBook-14:~/cuda/cuda$ nvcc matrix.cu --ptxas-options=-v
matrix.cu(78): warning #177-D: function "sum_of_matrix" was declared but never referenced
Remark: The warnings can be suppressed with "-diag-suppress <warning-number>"
ptxas info : 42 bytes gmem
ptxas info : Compiling entry function '_Z16block_mul_matrixPiS_S_' for 'sm_52'
ptxas info : Function properties for _Z16block_mul_matrixPiS_S_
0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info : Used 32 registers, 128 bytes smem, 344 bytes cmem[0], 4 bytes cmem[2]
ptxas info : Compiling entry function '_Z18multiply_of_matrixPiS_S_' for 'sm_52'
ptxas info : Function properties for _Z18multiply_of_matrixPiS_S_
0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info : Used 17 registers, 344 bytes cmem[0], 4 bytes cmem[2]
ptxas info : Compiling entry function '_Z13sum_of_matrixPiS_S_' for 'sm_52'
ptxas info : Function properties for _Z13sum_of_matrixPiS_S_
40032 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas info : Used 20 registers, 344 bytes cmem[0]
- 0 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads，是说本地TLM使用为0。
- Used 17 registers,是说使用了17个TLR。
- 40032 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads，是说使用了本地TLM 40032字节。
另外，关于存储器访问，如果一个warp中的32个thread访问连续的32字节的内存数据时，会触发内存访问合并，此时仅需要一次内存事务（即一条指令）就能完成访存操作。所以，我们希望全局内存访问能够合并，这样一个warp内对global memory的读写操作合并到尽可能少的内存事务中，使得内存带宽被尽可能充分利用，最终Kernel的执行速度更快。当warp中内存访问地址对齐32字节的整数倍且thread间访问内存连续时，内存访问效率最高。
TLM分配大小
TLM大小和总的线程数量有关，而总线程数量是由各级计算单元决定的。
totalTlmSize=properties.numHbmSections(16)*properties.numCUsPerSPC(4)*properties.numEUsPerCU(4)*properties.maxWarpPerEU(8)*properties.numThreadsPerWarp(32)*tlm_size_per_thread(4096) = 16x4x4x8x32x4096=256M.
nvprof
最后使用nvprof对程序进行profile,可以看到，大部分时间都花费在了HOST到DEVICE之间的拷贝上，实际的KERNEL运行时间不长。
Nsight profile
Nsight是NVIDIA面相开发者提供的开发工具套件，能提供深入的跟踪、调试、评测和分析，以优化跨 NVIDIA GPU和CPU的复杂计算应用程序。Nsight主要包含Nsight System、Nsight Compute、Nsight Graphics三部分。安装CUDA后，默认都已经安装。
设置好profile参数，CUDA程序路径，点击Start开始Profile:
可以看到主机端和GPU端 CUDA API执行的timeline:
host侧：
SM调度BLOCK的能力上限
从调试来看，每个SM对BLOCK调度能力是有上限的，将数据形状设置为如下参数进行调试：
可以看到，在每个BLOCK 32线程的参数配置下，虽然我们设置了256个BLOCK，并且每个SM的调度容量有64WARP，但是由于每个SM最大只能调度32个BLOCK，导致SM并没有打满，硬件利用率只有一半。
调整参数，每个BLOCK两个WARP，这样如果每个SM的调度上限是32个BLOCK，正好对应64个WARP，打满整个SM。
可以看到，此时正好打满整个SM。根据分析，可以绘画这样的ROOFLINE：
Kernel运行时间统计
kernel在DEVICE侧的运行时间通过cudaEventCreate/cudaEventRecord/cudaEventSynchronize/cudaEventElapsedTime/cudaEventDestroy接口统计。
cuda event事件是直接在GPU上实现的，因此时间不包括同时包含设备代码和主机代码的混合代码计时。
计时原理是在GPU工作负载前后插入两个计时command packet，start和stop,当stop执行完毕后（通过cudaEventSynchronize确保），说明start, stop以及中间执行的COMMAND PKT都已经执行完毕，这样，stop-start得到的时间就是实际的工作负荷时间。
总结：
优化时遵守的一些原则：
1.Grid一定要给足block.
2.Block内一定要给足thread,目的是提高并发WARP的数目，隐藏延迟.
3.Block内线程的数目一定是warpsize的整数倍.
每个warp的执行上下文（execution context，如程序计数器 和 寄存器等）在warp的整个生命周期内都被保存在片上内存（on-chip memory）。因此从一个执行上下文切换到另一个执行上下文是无开销的。在每个指令执行时间里，warp scheduler都会选择一个warp，这个warp中的所有thread都准备好执行它的下一个指令，warp scheduler会向这些thread发送指令并执行。那么等待这个warp准备好的这段时间（即时钟周期数）以及warp本身执行需要的数据ready之前，就是执行延迟（latency），为了掩盖这个延迟，warp scheduler可以在延迟期间发指令给其他的warp，因此一个SM内的warp数越高，延迟掩盖会越充分，利用率和性能也会越高。而一个SM上同时能并发存在多少个warp和block，是跟该SM上的寄存器数量和shared memory大小相关的（所以为了性能更高，需要把warp设置的多一些）。
参考资料
Cuda Stream:https://developer.download.nvidia.cn/CUDA/training/StreamsAndConcurrencyWebinar.pdf
https://www3.nd.edu/~zxu2/acms60212-40212-S16/Lec-11-GPU.pdf
https://cuda-tutorial.github.io/part2.pdf
基于ARM平台(Jetson)的CUDA编程入门_哔哩哔哩_bilibili
扫一扫
1095
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
