2020-02-20

CUDA

CUDA-再看规约-一段规约

一阶段规约

使用原子操作和shared memory的组合可以避免第二个kernel的调用。使用一个（global还是shared）变量记录哪个block已经做完了自己的工作。一旦所有的block都完成后，使用一个block执行最有的规约。

【CUDA专家手册 269页】

更简洁的一段规约

在之前的二段归约中，每一阶段的第一步是每个block的每个thread得到部分和，由于block间是不能同步的，所以才有了第二阶段。

现在是使用一个寄存器partialSum来存储每个thread对应的自己的部分和，最后使用原子操作对这些部分和再求和：

__global__ void reductionKernel(int* out, int* in, size_t N){
    int tid = threadIdx.x + blockDim.x * blockIdx.x;
    int stride = blockDim.x * gridDim.x;

    // 每个thread得到自己的部分和，
    int partialSum = 0;
    for (size_t i = tid; i<N; i+=stride){
        partialSum += in[i];
    }

    // 使用原子操作，求部分和的和
    atomicAdd(out, partialSum);
}

void func(int* res, int* in, size_t N){
    // 初始化
    cudaMemset(res, 0, sizeof(int));
    reductionKernel<<<numBlock, threadsPerBlock>>>(answer, in ,N);
}

注意一点，要写入out数组，那么out数组的初始值必须设为0；

2020-02-20

CUDA

CUDA-再看规约-循环展开

循环展开

循环展开减少了循环中的变量的判断，自加等其他操作，是提升性能一个手段。

1
2
3

for (int i=0; i<1000; i++){
    a[i] = b[i] + c[i];
}

将循环次数减少一半：

for (int i=0; i<1000; i+=2){
    a[i] = b[i] + c[i];
    a[i+1] = b[i+1] + c[i+1];
}

下面的两段规约法，展示使用完全的循环展开。

再看规约-两段规约

类似与点积中的两段法。两段规约同样使用shared memory。

过程见下图：

实现上图过程：

/*
* in: 待读取的数组。out: 待写入的数组。N: in数组的长度
*/
__global__ void reductionKernel(int* out, int* in, size_t N){
    // 此处是一个未分配大小的shared memory
    __shared__ int buffer[];

    // 第一阶段，
    int sum=0;
    int tid = threadIdx.x + blockDim.x * blockIdx.x;
    int stride = blockDim.x * gridDim.x;
    // 使用stride loop
    for (size_t i = tid; i<N; i+=stride){
        sum += in[i];
    }
    buffer[threadIdx.x] = sum;
    // 等待得到所有threads的sum值
    // block与block之间怎样同步呢？这里不需要同步
    __syncthreads();   

    // 第二阶段，规约求每个block的最终值
    for (int i = blockDim.x >> 1; i!=0; i>>=1){
        if (threadIdx.x < i){
            buffer[threadIdx.x] += buffer[threadIdx.x + i];
        }
        __syncthreads();
    }
    // 最后一步，将每个block对应的buffer中的第一个元素放入out数组
    // 此时out数组中存放的是每个block的最终值。
    if (threadIdx.x == 0)
        out[blockIdx.x] = buffer[0];
}

分析说明：

强调：第一阶段的blockDim.x应该是32的倍数，每个block对应的Shared memory大小也要是32的倍数。这样才能保证在第二阶段所处理的数据是整齐的。
第二阶段结束后block与block间结果怎样同步？kernel函数执行完毕，表示两个block都得到了结果，即同步了。所以一种方式的到最终结果是两次调用这个kernel函数。见下面code。
对out数组的最终处理，可以考虑将out传回host后在CPU中计算，就像向量点积。也可以使用一个block，第二次调用上述kernel函数，得到最终值。（我感觉不如传回CPU）。
shared memory 的访问索引是threadIdx.x。因为一个block共用一段Shared memory，所以Shared memory的不同地址可用这个block的threadIdx.x唯一标识。
使用blockDim.x >> 1和i>>=1，这属于指令优化。
kernel中的Shared memory可以不分配大小，调用kernel函数时再指定大小。

在Host中调用kernel函数：

void reduction(int* res, int* partial, int* in, 
                size_t N, int numBlocks, int numThreads){
    int sharedSize = numThreads*sizeof(int);
    // 第一个kernel得到所有block对应的buffer，
    // 这些block之间没有同步，所以只有所有的block的buffer都得到了
    // 这个kernel才时执行结束，也就是说，为了同步所有的block的buffer，
    // 只能再调用一次这个和函数了。
    reductionKernel<<<numBlocks, numThreads, sharedSize>>>
                    (partial, in, N);
    reductionKernel<<<1, numThreads, sharedSize>>>
                    (res, partial, numBlocks);
}

看到了吧：

block间的同步，其实是当该kernel函数执行结束后，自动同步。
第二次调用相同的kernel，只使一个block。
注意：此处在<< , , >>>中出现了第三个参数，由于在kernel中，没有指定shared memory的大小，所以在调用kernel函数时要制定大小了。
这个kernel中的第二个for循环中，循环的之后阶段，当threads数小于等于32时，可以特殊处理，此时的所有活动threads在一个warp中，所以不需要同步。细节间32特殊处理规约。
补充：每个block中的warp时按照lockstep的方式执行每条指令。

再看规约-第二阶段循环展开

优化意识：当threads数小于等于32时，也就是threadIdx.x从0到31时，可以特殊处理。这种处理就是循环展开。

/*
* in: 待读取的数组。out: 待写入的数组。N: in数组的长度
*/
__global__ void reductionKernel(int* out, int* in, size_t N){
    // 此处是一个未分配大小的shared memory
    __shared__ int buffer[];

    // 第一阶段，
    int sum=0;
    int tid = threadIdx.x + blockDim.x * blockIdx.x;
    int stride = blockDim.x * gridDim.x;
    // 使用stride loop
    for (size_t i = tid; i<N; i+=stride){
        sum += in[i];
    }
    buffer[threadIdx.x] = sum;
    // 等待得到所有threads的sum值
    __syncthreads();

    // 第二阶段，规约求每个block的最终值, 直到i=32.
    // 之后就是特殊处理了
    for (int i = blockDim.x >> 1; i>32; i>>=1){
        if (threadIdx.x < i){
            buffer[threadIdx.x] += buffer[threadIdx.x + i];
        }
        __syncthreads();
    }

    // 第三步，当元素个数小于等于32时：
    // 这个展开避免了循环控制的开销，和thread同步的逻辑
    if(threadIdx.x < 32){
        volatile int* wsSum = buffer;
        if(blockDim.x > 32)
            wsSum[threadIdx.x] += wsSum[threadIdx.x + 32];
        // 一下5条指令，一定都会被执行，不过也可以加上if条件。
        wsSum[threadIdx.x] += wsSum[threadIdx.x + 16];
        wsSum[threadIdx.x] += wsSum[threadIdx.x + 8];
        wsSum[threadIdx.x] += wsSum[threadIdx.x + 4];
        wsSum[threadIdx.x] += wsSum[threadIdx.x + 2];
        wsSum[threadIdx.x] += wsSum[threadIdx.x + 1];
        if (threadIdx.x == 0){
            volatile int* wsSum = buffer;
            out[blockIdx.x] = wsSum[0];
        }
    }
}

注意：当编写warp同步的code时，必须对shared memory 的指针使用volatile关键字。这个关键字告诉编译器必须将wsSum[threadIdx.x]的值存回global中。如果没有volatile，为什么不能得到正确结果？

第三步为什么正确，看下图之前，假设：

下图将大小为32的warp假设为4，即warp size=4：
所展示的情况是从上一个图片中的第一阶段的结果开始，即blockDim.x=8.

对于code中的第三步，见下图：

从上述code第二阶段的循环开始：

blockDim.x=8;
for(i=8>>1; i>4; i>>=1){
    i=4, i!>4, 退出循环;
}

对于 threadIdx.x<4，blockDim.x>4，即 0 1 2 3. 执行下面操作：
wsSum[0 1 2 3] += wsSum[(0 1 2 3) + 4];
wsSum[0 1 2 3] += wsSum[(0 1 2 3) + 2];
wsSum[0 1 2 3] += wsSum[(0 1 2 3) + 1];

说明：

其中执行单元只有这个block 中的第一个warp。
最终需要的是每个block对应buffer是第一个值，所以除第一个位置以外的位置，其值是无意义的。但是由于warp中的线程执行相同的指令，所以warp中其他threads也工作了，只不过这些结果我并不需要。
再强调：第一阶段的blockDim.x应该是32的倍数，每个block对应的Shared memory大小也要是32的倍数。这样才能保证在第二阶段所处理的数据是整齐的。

再次证明，用假设小的数值，用小的例子，小的数据有助于理解，

再看规约-第二阶段进一步循环展开

使用template把上述实现的第二阶段的for循环，做循环展开，这个展开是完全的展开。思路与上述实现一样，warp同步的优化方案就更近一步了。code如下：

template<int threadsPerBlock>
__global__ void reductionKernel(int* in, int* out, size_t N){
    __shared__ int buffer[];

    // 第一步，
    int tid = threadIdx.x;
    int sum = 0;
    int id = tid + threadsPerBlock * blockIdx.x;
    int stride = threadsPerBlock * gridDim.x; 
    for (size_t i = id; i< N; i += stride){
        sum += in[i];
    }
    buffer[tid] = sum;
    __syncthreads();
    
    // 第二步，循环展开, 从1024开始，
    // 每个block最多threads个数是1024（就目前的GPU架构）
    if(threadsPerBlock >= 1024){
        if(tid < 512)   // tid 小于512 的threads开始工作
            buffer[tid] += buffer[tid + 512];
        __syncthreads();
    }
    if(threadsPerBlock >= 512){
        if(tid < 256)  // tid 小于256 的threads开始工作
            buffer[tid] += buffer[tid + 256];
        __syncthreads();
    }
    if(threadsPerBlock >= 256){
        if(tid < 128)  // tid 小于128 的threads开始工作
            buffer[tid] += buffer[tid + 128];
        __syncthreads();
    }
    if(threadsPerBlock >= 128){
        if(tid < 64)  // tid 小于64 的threads开始工作
            buffer[tid] += buffer[tid + 64];
        __syncthreads();
    }
    // 第三步。warp同步操作，不需要__syncthreads()
    if( tid<32 ){     // tid 小于32 的threads开始工作
        volatile int* wsSum = buffer;
        // 下面的6条if语句都会被执行，
        // 看第一个if的64与第二步最后的64连接上了。
        if (threadsPerBlock >= 64) wsSum[tid] += wsSum[tid + 32];
        if (threadsPerBlock >= 32) wsSum[tid] += wsSum[tid + 16];
        if (threadsPerBlock >= 16) wsSum[tid] += wsSum[tid + 8];
        if (threadsPerBlock >= 8) wsSum[tid] += wsSum[tid + 4];
        if (threadsPerBlock >= 4) wsSum[tid] += wsSum[tid + 2];
        if (threadsPerBlock >= 2) wsSum[tid] += wsSum[tid + 1];
        if (tid == 0) out[blockIdx.x] = wsSum[0];
    }
}

完全展开，避免了循环控制的开销，第三步的展开还避免了thread同步的逻辑。

这里的执行有个特点，看if判断中的数字：1024，512，256，128，64，32，16，8，4，2，1，从逻辑看，一定是从大到小排列的，而且执行的时候，如果有一个if满足，那么这个if及其之后的if语句都会被执行。过程看下图。

第三次强调：第一阶段的blockDim.x应该是32的倍数，每个block对应的Shared memory大小也要是32的倍数。这样才能保证在第二阶段所处理的数据是整齐的。

所以设threadPerBLock为512时的示意图如下：

第三步中，不需要__syncthreads()，一个warp中的threads天然同步。

这个kernel是个模板kernel，在被host函数调用时，要指定threadsPerBlock的值。

调用kernel的host函数也必须是一个模板函数：

template<int threadsPerBlock>
void callReduction(int* res, int* partial, 
                    int* in, size_t N, 
                    int numBlocks){
    // 第一阶段，得到所有block对应的shared memory中的值
    reductionKernel<numBlocks><<<numBlocks, 
                                threadsPerBlocks, 
                                threadsPerBlocks*sizeof(int)
                                >>>(partial, in, N);
    // 两个阶段串行
    // 第二阶段，对上一步中partial中的所有值，执行规约
    reductionKernel<numBlocks><<<1, 
                                threadsPerBlocks, 
                                threadsPerBlocks*sizeof(int)
                                >>>(res, partial, numBlocks);             
}

// 使用switch语句，按照不同 block大小出发对应的模板函数。
void func(int* out, int* partial, 
            int* in, size_t N, 
            int numBlocks, int threadsPerBlock){
    switch( theadsPerBlock ){
        case 1: return callReduction< 1 >(...);
        case 2: return callReduction< 2 >(...);
        case 4: return callReduction< 4 >(...);
        case 8: return callReduction< 8 >(...);
        case 16: return callReduction< 16 >(...);
        case 32: return callReduction< 32 >(...);
        case 64: return callReduction< 64 >(...);
        case 128: return callReduction< 128 >(...);
        case 256: return callReduction< 256 >(...);
        case 512: return callReduction< 512 >(...);
        case 1024: return callReduction< 1024 >(...);
    }
}

之所以不能在一个kernel中获得最后结果，是因为，block与block之间是没有同步的。而第二阶段的规约需要所有block的buffer中第一个值，集合在一起，作为前提，如第一幅图所示。可不可以只一阶段就完成任务呢。

2020-02-20

CUDA

CUDA-ComputeCapacity-SM_version6.1的参数值

Compute Capability

对于CPU，不同架构的CPU有者不同的功能和指令集（MMX，SSE，SSE2）。而对于GPU，不同的功能由不同的Compute Capability表示。

NVIDIA GPU支持的Compute Capability有1.0，1.0，... 2.0，... 6.0，6.1，7.5。高版本的Compute Capability是低版本Compute Capability的超集，就是俄罗斯套娃式的嵌套结构。比如6.1版本支持1.0版本的所有功能。

比如，原子操作时硬件在内存上执行的。而只有1.1和1.1之后的版本才支持global memory上的原子操作，而只有1.2和1.2之后的版本才支持shared memory上的原子操作

在应用中，通常会指定最低Compute Capability版本，如2.3，告诉编译器，如果硬件支持的Compute Capability版本低于2.3，那么将无法执行这个和函数。做法是使用nvcc时增加一个选项nvcc -arch=sm_23，所以Compute Capability有称作SM版本。

SM版本的形式是X.Y，X表示核心的架构，7表示Volta，6表示Pascal。Y表示硬件的特性版本。

不同的SM版本所支持的特性，和技术参数，间这个页面的table14和table15：
https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#compute-capabilities

SM version=6.1 的关键参数

指标	值
Maximum dimensionality of grid of thread blocks	3
Maximum x- or y-dimension of a block	1024
Maximum z-dimension of a block	64

指标	值
Maximum number of threads per block	1024
Warp size	32
Maximum number of resident grids per device	32
Maximum number of resident blocks per multiprocessor	32
Maximum number of resident warps per multiprocessor	64
Maximum number of resident threads per multiprocessor	2048

指标	值
Number of 32-bit registers per multiprocessor	64 K
Maximum number of 32-bit registers per thread block	64 KB
Maximum number of 32-bit registers per thread	255

指标	值
Maximum amount of shared memory per multiprocessor	96 KB
Maximum amount of shared memory per thread block	48 KB
Number of shared memory banks	32
Amount of local memory per thread	512 KB
Constant memory size	64 KB
Maximum number of instructions per kernel	512 million

这些指标要心里有数

2020-02-20

CUDA

CUDA-原子操作-例-直方图

原子操作

atomicAdd(), atomicSub(), atomicXor()...

原子操作要排队，所以，能不用就不要使用。

原子操作-直方图

前面说过了，原子操作能不用就不使用。但是有些情况只能使用原子操作，当成千上万个threads同时修改同一个内存地址时，大规模并行系统会带来负担，在硬件中支持的原子操作可以减轻这种负担。比如计算“直方图”：计算一个数组中元素的出现频率。许多算法需要计算直方图，比如在图像处理使用过。

这里给出计算直方图的3种方法，CPU，未优化的GPU，优化后的GPU。

假如有一个很大的字符数组，实际中可能是像素的颜色值，或者是音频的采样数据。计算这个字符数组的直方图。由于每个随机的字符（一个char型占1 byte）都有2x2x2x2x2x2x2x2=256个不同的可能取值。所以直方图中要包含256个元素。

CPU版本

#define SIZE (100*1024*1024)

int main{
    // 1) 生成100MB的随机数据
    char* buffer = (char*)big_random_block( SIZE );
    // 2) 将buffer数组初始化为0
    int histo[256];
    for (int i=0;i<256;i++){
        histo[i] = 0;
    }
    // 3) 计算每个字符出现的频率
    for (int i=0;i<SIZE;i++){
        histo[buffer[i]]++;
    }
    // 4) 回收资源
    free(buffer);
    return 0;
}

不必解释。

未优化的GPU

计算直方图的特点是，读取buffer中的每一个元素，在histo数组中找到这个元素，后这个元素的频数加一，直到得到最后的直方图。如果使用并行计算，看到了吧，问题在哪？每一个thread读取buffer中不同的地址，取元素，没毛病。但是当写入histo数组时，会同时有多个threads写入histo的同一个位置，只必然造成Race condition（竞争）。所以多个threads写入同一个地址时，将写操作串行。先给出Host的主体：

int main( void ) {
    char *buffer = (unsigned char*)big_random_block( SIZE );

    char *dev_buffer;
    int *dev_histo;
    HANDLE_ERROR( cudaMalloc( (void**)&dev_buffer, SIZE ) );
    HANDLE_ERROR( cudaMemcpy( dev_buffer, buffer, SIZE,cudaMemcpyHostToDevice ) );

    HANDLE_ERROR( cudaMalloc( (void**)&dev_histo, 256 * sizeof( long ) ) );
    // dev_histo的初始化使用cudaMemset(), 将所有元素值设为0
    HANDLE_ERROR( cudaMemset( dev_histo, 0, 256 * sizeof( int ) ) );

    kernel<<<>>>();  //！！！

    int histo[256];
    HANDLE_ERROR( cudaMemcpy( histo, dev_histo, 256 * sizeof(int), cudaMemcpyDeviceToHost ) );

    // 验证GPU上计算的正确性，执行逆向运算，
    // 判断最后histo数组是否全为零（看是否可以会到初始值）
    for (int i=0; i<SIZE; i++)
        histo[buffer[i]]--;  
    for (int i=0; i<256; i++) {
        if (histo[i] != 0)
            printf( "Failure at %d!\n", i );
    }
    cudaFree( dev_histo );
    cudaFree( dev_buffer );
    free( buffer );
    return 0;
}

现在只剩下kernel函数。结果直方图数组有256个元素，所以考虑每个block含有256个threads。其实可以有其他的配置，比如100MB数据共有104,857,600个字节，所以可以启动一个block，让每个thread处理409,600个数据，或者可以启动409600个blocks，每个thread处理1个元素。

最有的配置方案介于这两种极端情况之间。实验发现，当block的数量是芯片SM数量的2倍时，将达到最有性能。

所以要先查硬件信息得到所使用GPU的SM数量：

cudaDeviceProp prop;
HANDLE_ERROR( cudaGetDeviceProperties( &prop, 0 ) ); // 使用0号GPU
int numSM = prop.multiProcessorCount;
kernel<<<numSM*2, 256>>>(dev_buffer, SIZE, dev_histo); // 未实现

下面实现kernel。

kernel的参数有三个：一个指向输入数组的指针，输入数组的长度，一个指向输出直方图数组的指针。

__global__ void histo_kernel( char *buffer,
                                long size,
                                int *histo ) {
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    
    // stride loop
    while (i < size) {
        atomicAdd( &(histo[buffer[i]]), 1 );
        i += stride;
    }
}

其中atomicAdd( &(histo[buffer[i]]), 1 )展示了如何使用原子操作：atomicAdd( add, y ), 这个操作读取地址add中值，将y加到add的值上。这个例子中add就是直方图中相应元素的位置。

性能结果怎样呢？实验得GPU版本比CPU版本性能差很多倍。分析原因：

第一，kernel函数实际上只包含了非常少的计算工作，而对global的访问远远多于计算，并且对global的访问相当慢。计算少，访存多，撞上性能瓶颈。
第二，当上万个threads访问少量（256）的内存位置时，发生严重的race condition。而且原子操作为了保证结果的正确性，对于相同的内存位置都将串行化。这导致排了非常长的队伍，因此抵消了并行带来的性能提升。

所以尝试优化GPU

使用shared memory优化原子操作

主要问题并不是使用了过多的原子操作，而是有太多的threads写入太少的地址造成的竞争。考虑使用两段操作。

为方便画图表示，假如元素只有x，y，z三种字符。那么自然地，每个block使用3个threads干活，对应的Shared memory大小也为3。假设block个数是2。两阶段如下图：

过程如下：

第一段，让每个并行的block计算它所处理的数据，将结果写入各自block对应的shared memory中。这样做的好处有：

第一，所有blocks之间的计算是并行的。
第二，写入shared memory效率远高于写入global memory。

第三，这种方式仍然需要原子操作，但是此时是256个threads在256个地址上的竞争，竞争程度显著减少。

所以，第一阶段计算得到每个block的临时直方图，比如，有三个blocks, 那么遍历完所有元素后会得到三个临时直方图。当得到所有三个临时直方图后（要同步），开始第二个计算阶段。

第二段，将得到的全部临时直方图合并为最终直方图，也应该是原子操作。这一段是，会有block数量个threads在256个地址上的竞争，竞争程度又显著减少。

__global__ void histo_kernel( char* buffer,
                                long size,
                                int* histo ) {
    // 初始化
    __shared__ unsigned int temp[256];
    temp[threadIdx.x] = 0;
    __syncthreads();

    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int offset = blockDim.x * gridDim.x;
    // 第一段计算
    while (i < size) {
        atomicAdd( &temp[buffer[i]], 1 );
        i += offset;
    }
    __syncthreads();

    //第二段计算
    atomicAdd( &(histo[threadIdx.x]), temp[threadIdx.x] );
}

// kernel 的配置如下
kernel<<<numSM*2, 256>>>(dev_buffer, SIZE, dev_histo);

在实际SP和threads执行中，上述代码（指令）被复制很多份，有多少个threads就复制多少份。这些拷贝的不同在于，threadIdx.x，blockIdx.x这些值会被自动赋值，0, 1, 2, 3, ...

性能如何？比CPU版本的快了一个数量级。

总结：使用了一种两阶段的算法，从而降低了在global memory上访存的竞争程度。这种降低竞争的策略总能获得不错的效果，所以当使用原子操作时，记住这种操作。

2020-02-20

CUDA

CUDA-常量内存提升性能-例-RayTracing

constants 常量必须在函数外声明
cudaMemcpyToSymbol()

这个函数是特殊版本的cudaMemcpy(), 唯一的不同是，cudaMemcpyToSymbol()将数据复制到常量内存constant memory，而cudaMemcpy()将数据复制到global memory。

常量内存为什么能提升性能-Ray tracing

__constant__ 把变量的访问限制为read-only，有了这个限制，必然获得某种回报。这个回报是：与从global memory中读取数据相比，从constant memory中读取相同的数据可以节约内存带宽，有两个原因：

对constant memory的单次读操作可以广播到其他临近的threads。

constant memory的数据将会缓存起来，因此相同地址的连续读操作不会产生额外的内存通信。

half-warp的广播是一把双刃剑。当16 个threads都读取constant memory相同的地址的时，性能极大地提升。但是当16个threads分别读取从constant memory不同的地址时，读操作会被串行化，性能会下降，这种情况就不如从global中读取。

从Ray-Tracing实例中体会：

要使用到的数据结构：

#define INF 2e10f
struct Sphere {
    float r,b,g;
    float radius;
    float x,y,z;
    // 假设有一个观察平面，那么（ox，oy）是这个观察平面中一个像素的坐标，
    // 这个方法将计算从这个像素中发射出的光线是否与这个球面相交。
    // 如果相交，则计算从这个像素点到这个球面的距离。当光线命中多个球时，
    // 只有最近的球面才会被看到。
    __device__ float hit( float ox, float oy, float *n ) {
        float dx = ox - x;
        float dy = oy - y;
        if (dx*dx + dy*dy < radius*radius) {
            float dz = sqrtf( radius*radius - dx*dx - dy*dy );
            *n = dz / sqrtf( radius * radius );
            return dz + z;
        }
        return -INF;
    }
};

kernel函数输入一个bitmap：

__global__ void kernel( unsigned char *ptr ) {
    // 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;
    float ox = (x - DIM/2);
    float oy = (y - DIM/2);
    float r=0, g=0, b=0;
    float maxz = -INF;
    // 每个线程考察所有的球，得到每个距离每个像素最近的球面。
    for(int i=0; i<SPHERES; i++) {
        float n;
        float t = s[i].hit( ox, oy, &n );
        if (t > maxz) {
            float fscale = n;
            r = s[i].r * fscale;
            g = s[i].g * fscale;
            b = s[i].b * fscale;
        }
    }
    ptr[offset*4 + 0] = (int)(r * 255);
    ptr[offset*4 + 1] = (int)(g * 255);
    ptr[offset*4 + 2] = (int)(b * 255);
    ptr[offset*4 + 3] = 255;
}

for循环中，每一个thread循环球集合Sphere *s中的所有球s[i]，得到每个像素的一个颜色值。这里的模式是：所有threads都会只读相同的存储地址。考虑上述的关于constant的特点，使用constant memory 进行优化。

优化前的main函数如下：其中s是声明在global memory中的。用constant memory优化只需要将下面code中的<1><2>分别修改为

<1> __constant__ Sphere s[SPHERES];

<2> HANDLE_ERROR( cudaMemcpyToSymbol( s, temp_s, sizeof(Sphere) * SPHERES) );

#include "cpu_bitmap.h"
#define rnd( x ) (x * rand() / RAND_MAX)
#define SPHERES 20
Sphere *s; // <1>

int main( void ) {
    CPUBitmap bitmap( DIM, DIM );
    unsigned char *dev_bitmap;
    // allocate memory on the GPU for the output bitmap
    HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() ) );
    // allocate memory for the Sphere dataset
    HANDLE_ERROR( cudaMalloc( (void**)&s, sizeof(Sphere) * SPHERES ) ); // <2>

    Sphere *temp_s = (Sphere*)malloc( sizeof(Sphere) * SPHERES );
    for (int i=0; i<SPHERES; i++) {
        temp_s[i].r = rnd( 1.0f );
        temp_s[i].g = rnd( 1.0f );
        temp_s[i].b = rnd( 1.0f );
        temp_s[i].x = rnd( 1000.0f ) - 500;
        temp_s[i].y = rnd( 1000.0f ) - 500;
        temp_s[i].z = rnd( 1000.0f ) - 500;
        temp_s[i].radius = rnd( 100.0f ) + 20;
    }

    HANDLE_ERROR( cudaMemcpy( s, temp_s, sizeof(Sphere) * SPHERES,
                                cudaMemcpyHostToDevice ) );
    free( temp_s );

    // generate a bitmap from our sphere data
    dim3 grids(DIM/16,DIM/16);
    dim3 threads(16,16);
    kernel<<<grids,threads>>>( dev_bitmap );

    // copy our bitmap back from the GPU for display
    HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap, bitmap.image_size(),
                                cudaMemcpyDeviceToHost ) );
    bitmap.display_and_exit();
    // free our memory
    cudaFree( dev_bitmap );
    cudaFree( s );
}

上述过程使用constant memory保存只读对象。感受这个模式：每个threads都访问相同的只读数据时，将获得额外的性能提升. 前面说了的两个原因：第一，这种模式将读取操作在半个warp中广播，第二，芯片上包含了常量内存缓存。

在许多算法中，内存带宽都是瓶颈，因此要时刻想着改善这种情况，

2020-02-20

CUDA

CUDA-死锁-例-点积

好好体会-Dot product

（CUDA by example 55页）

这个例子值得好好感悟

点积：两个长度相同的向量A和B，对应元素相乘后相加。

CUDA实现思路：每个线程分别读取A和B中对应位置的元素，紧接着执行相乘操作，最后将相乘结果存入shared memory的对应位置。

注意，当向量元素个数远远超过一个block中的threads数量时的处理。

const int N = 33*1025;
const int threadsPerBlock = 256;
const int blocksPerGrid = (N+threadsPerBlock-1)/threadsPerBlock;

__global__ void dotProduction(float* a, float* b){
/// 第一步：
    // thread计算得到a和b对应元素的乘积，
    // 存入这个thread对应的shared memory中的位置
    // 每个block对应的shared memory的大小为这个block中的threads数量
    __shared__ float cache[threadsPerBlock];

    int tid = threadIdx.x + blockDim.x*blocIdx.x;

    // 每个block的shared memory的索引就是threadIdx.x，与blockIdx.x 无关,
    // 这个block和那个block中的thread的ID是一摸一样的。
    int cacheIndex = threadsIdx.x;

    // 执行乘法操作，更新threads ID，
    // 同stride-loop
    // 看图一的过程
    float tmp = 0.0f;
    while (tid<N){
        tmp += a[tid]*b[tid];
        tid += blockDim.x*gridDim.x;
    }
    cache[cacheIndex] = tmp;

    // 同步这个block中的所有threads
    __syncthreads();
    // 确保所有threads完成工作之后，执行后续指令

/// 第二步：
    // 对于每个block对应的shared memory中的元素，进行规约求和。
    // 其中threadPerBlock必须是2的指数。
    // 同for-loop
    int i = blockDim.x/2;
    while(i != 0){
        if (cacheIndex <i )
            cache[cacheIndex] += cache[cacheIndex + i];
        // 确保上一轮所有和得到，所以要同步
        __syncthreads();
        i = i >> 2;
    }
    // 把每个shared memory中的第一个元素，也就是这个shared memory中
    // 所有元素之和，写入c中对应的位置。
    if(cacheIndex == 0)
        c[blockIdx.x] = cache[0];
}
...

// 在主函数中，把c从device拷贝到host，以及之后：
cudaMemcpy(h_c, c, cudaMemcpyDeviceToHost);

// 在CPU上将h_c中的结果求和，于此同时GPU上可以后其他任务执行。
// CPU和GPU并行执行。隐藏延时
float sum = 0;
for (int i=0;i < blocksPerGrid; i++){
    sum += h_c[i];
}
// sum即是最终点积结果。

其中这一部分：

float tmp = 0.0f;
while (tid<N){
    tmp += a[tid]*b[tid];
    tid += blockDim.x*gridDim.x; // 更新tid，自加不是1，而是所有threads数量。
}
cache[cacheIndex] = tmp;

当所有threads的数目小于a或b中的元素个数时（不论有多少个blocks），上述保证正确，与stride更新效果相同。当threads个数等于元素个数时，也正确。所以这样写，分析见下图：

上图中只使用了一个block，所以在第二步归约计算时，就可以在第0个位置上得到最终结果。当使用多个blocks时，第二步得到每个block的第0个元素，而这些若干个第0个元素保存于c（Global）中，最终还要将c中元素求和。

而下面这种实现：

1
2
3

if (tid<N){
    cache[threadIdx.x] = a[tid]*b[tid];
}

只适用于thread个数等于元素个数时（如下图）。但通常元素个数会远大于threads数。所以不适用此法。

技能：多个blocks中的各个shared memory 缓存同时被操作。
为什么要将最后的结果传回host计算？

因为，事实证明，想GPU这种大规模并行机器在执行最后的规约步骤时，通常会浪费计算资源，因为此时的数据集非常小。比如，当使用480 个threads将32 个数相加时，将难以充分使用每一个threads。

总结一下，使用shared memory优化dot-production为什么有效，因为减少了写入global memory的次数，并且复制回host的数据量减少。性能增加。

敲黑板注意threadIdx.x 与threads ID的区别，前者相对ID后者绝对ID。访存Shared memory时，一定使用threadIdx.x。

__syncthreads() 放错位置会导致死锁

规约程序中：

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

如果将__syncthreads()放入if语句，会产生死锁：

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

解释一下，当出现线程发散时，发散的分支会使得某些threads处于空闲状态，而其他threads将执行分支中的代码。而对于__syncthreads()而言，CUDA架构确保，一个block中的所有threads都执行到__syncthreads()，才能执行__syncthreads()之后的语句。这样一来，上述代码块，只要有一个threads没有执行if语句，也就不能够执行__syncthreads()，其他执行了if语句的threads会等待哪一个thread，一直等下去，造成死锁。

所以，对于__syncthreads()要谨慎使用。

2020-02-20

CUDA

CUDA-存储优化-例-矩阵转置

存储优化

重叠内存传输和计算，组团传输减少小块数据的频繁传输。

现代GPU对访存做了优化，使得不是严格的Coalescing Access也是可以接受的。但是

避免凌乱无规律的访存
避免一个线程访问连续的一段空间。

将global memory中的数据放入shared memory中，使得数据的位置相邻后写入global memory，此时相邻的threads就可以访问相邻的地址，满足Coalescing Access。shared memory的设计目的之一就是通过对它的编程，来规则化访存模式。

回忆：shared memory的架构是连续32 bits(即4 Bytes)的地址被分配到连续的bank中。见下图：

原图来自这里。里面有不少其他资源可参考。

对于bank冲突，注意理解地址的编号和bank的编号的不同。

注意

Coalescing access是就global memory而言的。

存储优化，优化实列：矩阵转置

未优化的矩阵转置

对global memory的读和写总有一个，对global存储的地址访问不是连续的，这就不能最大化coalescing access。看下图：

每个block负责矩阵的一个子块（tile），所有子块并行执行。TILE_DIM=blockDIm.x。

__global__ void transpose(float* a, float* b){
    // 每个线程的id
    int idx = threadIdx.x + TILE_DIM * blockIdx.x;
    int idy = threadIdx.y + TILE_DIM * blockIdx.y;

    // 每个数据转置前后，在矩阵中的索引
    int index_a = idx + WIDTH * idy;
    int index_b = idy + HEIGHT * idx;

    // 转置操作
    b[index_b] = a[index_a];
}

所以：

使用shared memory优化的矩阵转置

在写入global memory之前，先从global中将所读取的元素存入shared memory中，当shared memory中有了所有元素，将此时的shared memory转置一下，最后将结果再写（合并地）入global memory。

如此一来对global的读和写，地址都是被连续访问的。看下图：

实现：

__global_ void transpose(float* a, float* b){

    __shared__ float tile[TILE_DIM][TILE_DIM];

    // 每个线程的id
    int idx = threadIdx.x + TILE_DIM * blockIdx.x;
    int idy = threadIdx.y + TILE_DIM * blockIdx.y;
    // idx 与 idy 线性组合 得到原矩阵的Index
    int index_in = idx + idy * WIDTH;  

    // 每个数据转置前后，在矩阵中的索引
    idx = threadIdx.x + TILE_DIM * blockIdx.y;
    idy = threadIdx.y + TILE_DIM * blockIdx.x;
    // // idx 与 idy 线性组合 得到转置后矩阵的Index
    int index_out = idx + idy * HEIGHT;

    // 从a中写入到bank中，会产生冲突。(马上解决)
    tile[threadIdx.y][threadIdx.x] = a[index_in];

    // 等待这个block对应的tile中的 元素都有了之后，再执行下后续操作。
    __syncthreads();

    // 转置操作
    b[index_out] = tile[threadIdx.x][threadIdx.y];
}

在Shared memory中存在bank冲突，如何解决。

解决上述过程中的bank冲突

对于上述实现过程中tile[threadIdx.y][threadIdx.x] = a[index_in];，如果tile的大小是16x16，那么这一句会产生16路的bank冲突。！！！如何理解bank的编号！！！

如何解决bank冲突，看下图：

实现中只需更tile的定义为__shared__ float tile[TILE_DIM][TILE_DIM+1];，列中的数据存于相同的bank。

如此一来，不论是对tile行或列的访存，都不会产生bank冲突。相同的颜色在不同的行和列。

2020-02-20

CUDA

CUDA-例-规约

优化实列：并行规约（parallel reduction）

规约的更为一般的形式：对一个输入数组进行某种操作，会产生一个更小的结果数组。比如点积算子，累加，min，max，平方和，逻辑与，逻辑或等等。规约的成立前提是，这些算子中的二元操作符合结合率。

未优化的规约（为便于图中展示，假设warp大小为2）

过程看图：

其中

id既是存储地址的id，也是threads的id。（因为threads的id此处没有更新）。
n个元素需要lg(n)次平行。
上图从上到下，是在时间上展开，并没有新的空间被使用。

warp大小为2，从一开始的第一次并行，就存在divergence。

过程实现（假如在shared memory中实现）：

__global__ void func(){

    __shared__ float partialSum[];
    // ...将数据放入shared memory中

    // 因为只使用一个block，所以只需要threads 的id
    int x = threadIdx.x;

    // 循环计算每一层，stride 为 1，2，4
    for(int stride=1; stride < blockDim.x; stride *= 2){
        // 对指定的thread 进行加操作，与id和stride有关，拿笔画画就找到规律。
        if (x % (2*stride) == 0)
            partialSum[x] += partialSum[x + stride];

        // 这一层都求和结束后，才能进行下一步
        __syncthreads();
    }
}

上述过程中，（看图）每次循环会规律的有线程没有做实际的工作，这些threads也在工作（因为是在同一个warp中），但是没有实际操作数。每一轮实际所需的线程数在减半。

优化后规约（为便于图中展示，假设warp大小为2）

由上一个实现中，看到了，每一轮的实际工作线程数在减半，但是实际上所有的threads都在工作，很多是没有意义的工作。

这样为什么不好？因为违反了Coalescing Access：相邻的线程处理相邻位置的数据.

所以改进如图：

其中：

id既是存储地址的id，也是threads的id。（因为threads的id此处没有更新）。
n个元素需要lg(n)次平行。
当数据元素更多时，橙色虚线框对应的threads与其他threads（很多时候）不属于同个warp，所以threads所用资源较前一个实现提前释放。

warp大小为2，图中所有次并行计算都不存在divergence。

实现：

__global__ void func(){

    __shared__ float partialSum[];
    // ...将数据放入shared memory中

    // 因为只是用一个block，所以只需要threads 的id
    int x = threadIdx.x;

    // 循环计算每一层，stride 为 4，2，1
    for(int stride=blockDim.x/2; stride > 0; stride /= 2){
        // 对特定的thread 进行加操作，与id和stride有关，拿笔画画就找到规律。
        if (x < stride)
            partialSum[x] += partialSum[x + stride];

        // 这一层都求和结束后，才能进行下一步
        __syncthreads();
    }
}

为什么此法好？

每一轮都有一半的 thread不需要工作，资源释放掉。让warp提前完工，释放资源。
block中的warp，没有了分支发散，或者说是最小化了分支发散。回忆warp：
- 在一个block中，连续的32个threads一组构成一个warp;
- warp 是最基本的调度单元
- warp中的threads在同步执行相同的指令（SIMT）
- warp中threads需要执行不同的路径时（分支发散），warp中每个threads都要执行所有的分支，因为是同步的。比如一个宿舍的6个学生可以是一个warp，今天有的想先上厕所，后吃饭，而有的不需要上厕所，此时所有的同学都会一起先上厕所，后一起吃饭。
- warp之间是没有关系的。
- warp间的切换时没有代价的。多warp工作可以隐藏延时
- warp的分割，连续32个threadIdx.x 为一组（一个warp）0到31，32到63，…
- warp中的分支发散不总是问题，但是如有很多分支语句的话，每个thread就需要执行所有的分支。
- 在设计程序时，不应个这样：if (threadIdx.x > 15){...} 而应该这样if (threadIdx.x > 32-1) {...}。后者，表示说，第一个warp干他的事，第二个warp执行这个分支。

就这个优化，当元素的个数小于32（warp的大小）时，不可避免的会产生分支发散。

2020-02-20

CUDA

CUDA-例-矩阵相乘

矩阵相乘CPU版本

思路：三重循环，

外层：遍历结果矩阵所有行，
中层：对结果矩阵每一行，遍历所有列，
内层：结果矩阵某一行某一列元素的得到，需要A矩阵对应行，B矩阵对应列，的元素相乘后累加，这需要一个循环。内层循环结束后，就得到结果矩阵中的一个元素。

看图：

实现：

void matrixMultiplicationCPU(float* M, float* N, float* P, 
                                int A_height, int Awidth_Bheight, int B_width){

    for (int i=0; i < A_height; i++){
        for (int j=0; j < B_width; j++){
            float sum = 0;
            for (int k=0; k < Awidth_Bheight; k++){
                float a = M[i * Awidth_Bheight + k];
                float b = N[k * B_width + j];
                sum += a * b;
            }
            P[i * B_width + j] = sum;  
        }
    }
}

其中，M的大小为（A_height x Awidth_Bheight），N的大小为（Awidth_Bheight x B_width）. 注意中括号内，索引的计算。

未优化的矩阵相乘kernel

参数：

AColBRow：AxB中A的列数，也是B的行数
BCol：B的列数
a，b矩阵的元素索引按行标注索引

__global__ void simpleMultiply2(float *a, float* b, float *c){

    int row = threadIdx.y;
    int col = threadIdx.x;
    float sum = 0.0f;
    for (int i = 0; i < AColBRow; i++) {
        sum += a[row * AColBRow + i] * b[i * BCol + col];
    }
    c[row * BCol + col] = sum;
}

这种方法中每一个线程负责乘累加结果的一个元素，这个元素的得到是通过循环读取a矩阵的一行，b矩阵的一列，后乘累加得到。这个过程其实就是把CPU版本中的外循环和中层循环去掉，替代这两个循环的是threadIdx.x和threadIdx.y，这两个值自加。对于a和b的访存索引的计算是不变的。因为不涉及到结果间的依赖，所以不需要同步机制。

缺点：

观察row和col的定义，发现这个实现只使用了一个block。对于大的矩阵相乘，一个block不足以覆盖所有元素。
矩阵相乘的过程中，同一个Global memory地址被频繁访问，所以速度慢。

多blocks分块处理矩阵相乘

参数：

AColBRow：AxB中A的列数，也是B的行数
BCol：B的列数
a，b矩阵的元素索引按行标注索引

__global__ void simpleMultiply2(float *a, float* b, float *c){

    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    float sum = 0.0f;
    for (int i = 0; i < AColBRow; i++) {
        sum += a[row * AColBRow + i] * b[i * BCol + col];
    }
    c[row * BCol + col] = sum;
}

此情况的线程结构是2D block，2D thread，注意线程组织的index，事实上是这样的：

注意：横向，index的第一个分量递增；纵向，index的第二个分量递增。

与上一种方法相比，优势是，使用多个blocks，每个block负责结果矩阵中的一个矩形块中的元素。
每个thread的操作细节看下图：

注意，C/C++中矩阵元素是按行排列的。假设线程组织是2D thread 2D block，假设tile大小是2x2的，假设使用4个blocks。那么跟踪上述code可得到结果矩阵中的第14个元素：

row = 3;
col = 2;
i=0: a[row * 3 + i] * b[i * 4 + col] = a[9] *b[2];
i=1: a[row * 3 + i] * b[i * 4 + col] = a[10]*b[6];
i=2: a[row * 3 + i] * b[i * 4 + col] = a[11]*b[10];
GET sum;
WRITE sum INTO c[row * 4 + col] = c[14];

看出逻辑及结果是正确的。

有了kernel，如何调用kernel：

// block的个数
dim3 grid((A的行数+1)/tile的行数，(B的列数+1)/tile的列数)
// 每个block中的线程数
dim3 block(tile的行数，tile的列数)
simpleMultiply2<<<grid, block>>>(a, b, c)

比如图中

定义tileRow=2，tileCol=2.
A的行数为ARow=4，列数ACol=3；
B的行数为BRow=3，列数BCol=4；

那么，有如下对kernel的配置，自然地每个tile的大小为2x2，每个block大的大小也是2x2. block与block之间相互不影响，所以所有blocks并行执行。

// block的个数，共有4个
dim3 grid(2，2)
// 每个block中的线程数，有4个
dim3 block(2，2)
simpleMultiply2<<<grid, block>>>(a, b, c)

其中tile是一个block处理区域的大小，也就是一个block的大小，比如2x2，3x4. 是开发者根据问题情况自己定义的（通常是32 的倍数）。通过结果矩阵的行列数和tile的行列数，可以计算出，所需block数量个blocks阵大小。

总结：这个优化算然可以处理任何大小的矩阵相乘，但是对global memory的读写并没有做优化。

Global memory慢

假如一个GPU的性能峰值时346.5 Gflops。所以需要346.5x4bit=1386 GB/s的带宽才能达到峰值。但是这个GPU的存储实际带宽可能只有86.4GB/s，也就是86.4GB/s / 4bit=21.6 Gflops的性能。而实际上代码运行的速度可能只有15 Gflops。只有峰值的1/20. 所以要尽量减少对global memory的访存。

对于矩阵相乘，可以使用shared memory来减少对global memory的访存。始终要有这个意识。

优化实列：使用shared memory优化矩阵相乘

Shared memory这么使用：

思路：

因为使用了若干个block来分块结果矩阵，每一个block有自己单独的shared memory，这个存储空间由这个block内的所有threads共享，而且就矩阵相乘这个问题，A矩阵的每一行，和B矩阵的每一列在计算这个块时会不止一次被使用。所以：

将结果阵中这一块对应的A中行和B中列存入这个block对应的shared memory中，如此避免了对global memory的频繁访问。如下图：
进一步，由于A和B会很大，存入shared memory中的A行，和B列也是分段的，避免了shared memory不够用。如下图：

其中橘色是这个block的Shared memory空间，循环所有的Asub和Bsub，先绿色后紫色，将每次的元素放入Shared memory，求得结果累加到Cvalue。循环结束后将最终的Cvalue写入最终位置，见图中箭头，

kernel函数包含二重循环：

外层，循环所有的Asub和Bsub，并将其放入Shared memory中。
内层，循环求Asub每部分行和Bsub每部分列的部分乘累加和Cvalue。

Bare in mind，每一个thread处理一个结果阵中的一个元素，也就是说，kernel函数是一个thread得到一个结果矩阵中的一个元素。

官方文档中有官方的解释，可以参考。

注意，整个过程中有两处同步操作，因为涉及到线程间的同步。

程序如下：

__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C)
{
    // block的行ID和列ID
    int blockRow = blockIdx.y;
    int blockCol = blockIdx.x;

    // 每个block计算结果矩阵中的一块
    Matrix Csub = GetSubMatrix(C, blockRow, blockCol);

    // 每个thread计算结果矩阵中的一个元素，
    // 这个结果通过累加存入Cvalue
    float Cvalue = 0;

    // thread的行ID和列ID
    int row = threadIdx.y;
    int col = threadIdx.x;

    // Loop over all the sub-matrices of A and B that are
    // required to compute Csub
    // Multiply each pair of sub-matrices together
    // and accumulate the results to Cvalue
    for (int m = 0; m < (A.width / BLOCK_SIZE); ++m) {

        // 得到这个block在这次循环中的A的子矩阵
        Matrix Asub = GetSubMatrix(A, blockRow, m);

        // 得到这个block在这次循环中的B的子矩阵
        Matrix Bsub = GetSubMatrix(B, m, blockCol);

        // 保存A的子矩阵和B的子矩阵
        __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

        // Load Asub and Bsub from device memory to shared memory
        // 每个thread加载A的子矩阵中一个元素，和B的子矩阵中的一个元素
        As[row][col] = GetElement(Asub, row, col);
        Bs[row][col] = GetElement(Bsub, row, col);

        // 将A中这部分所有的元素，和B中元素都写入到shared memory中后，
        // 才能开始下面的操作。
        __syncthreads();  

        // 计算当前A的子矩阵与B的子矩阵，相乘后累加结果。
        // 当外层循环结束，表示累加每一块累加的结果完成，
        // 得到最后的C中这个位置的结果。
        for (int e = 0; e < BLOCK_SIZE; ++e)
            Cvalue += As[row][e] * Bs[e][col];

        // Synchronize to make sure that the preceding
        // computation is done before loading two new
        // sub-matrices of A and B in the next iteration
        // 等待这一阶段的 加类乘完成后，才能够进行下一次循环，更新shared memory中的内容。
        __syncthreads();
    }

    // Write Csub to device memory
    // Each thread writes one element
    SetElement(Csub, row, col, Cvalue);
}

那么对global memory 的访问减少了多少倍呢?

假如一个block计算2x2的结果C的子矩阵，那么未优化的实现要访问Global memory (2*2*2)8次行/列（定义1次行/列：访问A的完整一行或B的完整一列）。优化后只访问4次行/列。8/4=2，就是tile的宽。
假如一个block计算3x3的结果C的子矩阵，那么未优化的实现要访问Global memory (2*3*3)18次行/列。优化后只访问6次行/列。18/6=3，就是tile的宽。
假如一个block计算100x100的结果C的子矩阵，那么未优化的实现要访问Global memory (2*100*100)20,000次行/列。优化后只访问200次行/列。20,000/200=100，就是tile的宽。

所以对于很大的矩阵，优化前后的性能差别是巨大的。相差tile宽，这么多倍。

相差tile宽，这么多倍。极限情况下，如果shared memory有足够的空间，可以把A，B都放入每个block的shared memory中。此情况下对global memory的访问最少。

上述优化方法中确定块tile的宽度要考虑shared memory的大小：

尽量最大化tile，block的大小

假如一个GPU每个SM有16kB的shared memory，并且每个SM有8个blocks。可以算出每个block可以使用2kB的shared memory。

其中1kB存矩阵A的相关元素，1kB存B的相关元素。假设数据类型是float，占4B。所以可以计算出1kB=16x16x4。block的大小为16x16. 也就是一个block有16x16个threads。

这种情况充分使用了SM的threads和shared memory资源，最大化资源占用率。

2020-02-20

CUDA

CUDA-线程同步-线程调度

线程同步（一个block内的同步）

一个block内的所有threads有时候是需要同步的（如使用shared memory优化的矩阵相乘中），方法是在kernel函数中的适当位置加上__syncthreads()。

当一个thread执行到__syncthreadas()时，这个thread会看它所在的block内的其他所有threads情况，如果发现还有其他threads没有执行到这个位置，则这个thread等待其他threads。直到block中所有active的threads都执行到此，接着向下执行。

注意这个同步只是这个block内的threads同步。而非是全局同步，CUDA中没有全局通过不的原因是，全局同步的系统开销会很大。

并行系统中负载均衡：要求线程的执行时间尽量接近。如果不均衡，在需要线程同步时，所有线程会等待最慢的那个thread，此时整体的速度就时最慢的那个thread的速度，其他速度快的threads 的优势完全没有了。

block与block之间是异步的，不存在相互等待（当）。而同一个warp内的threads天然同步。

注意线程同步，可能导致死锁，比如下情况：
1
2
3
4
5
6
if (func()){
__syncthreads();
}
else{
__syncthreads();
}
产生死锁，执行不同分支的threads相互等待，谁也等不到谁。
block间的同步

block与block之间是异步的，当前kernel执行结束后，block之间自然同步了。言外之意是kernel的执行时间是由最慢的哪个block决定，所以，上文强调负载均衡。

线程调度
1. SP和threads
  
  为什么GPU的threads数量远远多于物理执行单元（SP）。因为每个SM中与CPU一样也含有上下文空间，用于执行上下文切换，从而实现多线程。（？提升吞吐量？）
  
  这里有个SP和threads的关系。以GPU G80为例：
  
  G80 的硬件信息：16个SM，每个SM含有8个SP，（共有16x8=128个SP），每个SM最多驻扎768个threads，总共同时执行12288个threads。（所以可以通过每个SM中最多可以驻扎的threads数，除以每个SM中的SP数，就得到了）
  
  解释：
  - 16，表是芯片实际含有16个SM
  - 8，表示每个SM含有8个SP（Streaming Processors），真正执行指令的工人。
  - 12288，表示这个芯片上可以同时有12288个threads进行调度。调度不意味着一定要实际执行。
  - 128个SP，表示每个时钟周期内实际并行执行的指令流为128个。但总共有12288个指令流间不停的切换。切换的目的是达到延时隐藏效果。
    
    warp的调度是零开销的。因为warp的上下文是存在与物理空间中的，需要了这个warp干活时，程序切换到这个warp上去即可。
    
    warp中的所有threads执行相同的指令。当有分支时，由于warp内threads天然的同步，所以含有分支时的执行会有性能下降。
2. warp和SP
  
  每个warp含有32个threads，假如每个SM只有8个SP，此时一个SM如何调度一个warp?
  - 第一个周期内，有8个threads进入SP
  - 第二，三，四个周期SP各进入8个threads
  - 如此循环，直到所有指令执行完毕
  - 所以此情况一个SM调度一个warp，需要4个周期（4个周期只是调度完成，指令执行完成需要4的倍数个周期）
    
    现代GPU的SP数已经远远大于32了。所以不存在上述问题了。
3. 调度warp实现延时隐藏
  
  一个实例：有一个kernel含有
  - 一个对global memory的读操作，这个操作耗时200个时钟周期。
  - 4个独立的乘或加操作，一个乘或加操作耗时4个时钟周期。
    
    那么需要多少个warp才可以将对global memory访问的延迟隐藏掉?
    
    首先每个warp需要执行4个独立的乘或加操作，共耗时4x4=16个始终周期。要覆盖200个时钟周期，就需要200/16=12.5，即13个warp串行，才能有效隐藏对global memory的访问延时。
    
    回忆：每个SM一次只能执行一个warp。待确认。。。。