NVIDIA CUDA 高度并行处理器编程（四）：性能优化习题

1.性能优化里面的两个归约 kernel 函数在使用线程时很浪费，每个块中有一半线程根本没执行。修改 kernel 来消除这种浪费现象。给出 kernel 启动时与之相关的配置参数的值。在运算操作方面是否需要更多额外成本？这种修改能消除哪种类型的资源限制？

答：修改后的第一个：

unsigned int t = threadIdx.x * 2;
    for(int stride = 1;stride <= blockDim.x;stride <<= 1){
        if(t % (2*stride) == 0){
            partialSum[t] +=partialSum[t + stride];
        }
        __syncthreads();
    }

修改后的第二个：

__shared__ double partialSum[];
unsigned int t = threadIdx.x;
for (unsigned int stride = blockDim.x; stride > 0; stride >>=1 ) 
{
 	__syncthreads();
 	if (t < stride)
 		partialSum[t] += partialSum[t+stride];
}

启动配置：

reduction_sum<<<ceil((double)n/THREAD_LENGTH/2), THREAD_LENGTH>>>(d_A, n);

第一个将线程 ID 的二倍赋给 t，同时循环结束条件由 blockDim.x/2 变成了 blockDim.x ，所以多了两个额外操作。

第二个将循环中初始条件：stride乘 2，所以多了一个循环的额外操作。

2.将上一题修改后的kernel函数进行比较，那种修改方案引入的运算操作较少？

第二种比第一种少了线程 ID 的乘法。

3.在习题1的基础上编写一个完整的kernel函数：（1）添加几条语句，实现把输入数组的部分数据从全局存储器加载到共享存储器中；（2）利用变量 blockIdx.x 让多个块作用于数组的不同部分；（3）根据 blockIdx.x 将这每部分的归约值写入到一个位置中。

kernel：

__global__ void reduction_sum(double *X, size_t input_size){
    __shared__ double partialSum[2 * THREAD_LENGTH];
    int i = 2 * blockIdx.x * blockDim.x + threadIdx.x;
	if(i < input_size) partialSum[threadIdx.x] = X[i];
    else partialSum[threadIdx.x] = 0.0;
    if(i + blockDim.x < input_size) partialSum[threadIdx.x + blockDim.x] = X[i + blockDim.x];
    else partialSum[threadIdx.x + blockDim.x] = 0.0;
    __syncthreads();
    unsigned int t = threadIdx.x;
    for(int stride = blockDim.x; stride >= 1; stride /= 2){
        if(t < stride)
            partialSum[t] += partialSum[t + stride];
            __syncthreads();
    } 
    // unsigned int t = threadIdx.x * 2;
    // for(int stride = 1;stride <= blockDim.x;stride <<= 1){
    //     if(t % (2*stride) == 0){
    //         partialSum[t] +=partialSum[t + stride];
    //     }
    //     __syncthreads();
    // }
    if(t == 0){
        X[blockIdx.x] = partialSum[t];
    }
}

4.以3中的客人了为基础设计一个归约程序。主机代码包括：（1）将大输入数组传入全局存储器中；（2）利用反复循环调用习题3中编写的kernel函数并调整配置参数的值，以便产生整个输入数组的归约结果。

#include<cuda.h>
#include<stdlib.h>
#include<stdio.h>
#define THREAD_LENGTH 1024
//the same as before 
__global__ void reduction_sum(double *X, size_t input_size){
    __shared__ double partialSum[2 * THREAD_LENGTH];
    int i = 2 * blockIdx.x * blockDim.x + threadIdx.x;
	if(i < input_size) partialSum[threadIdx.x] = X[i];
    else partialSum[threadIdx.x] = 0.0;
    if(i + blockDim.x < input_size) partialSum[threadIdx.x + blockDim.x] = X[i + blockDim.x];
    else partialSum[threadIdx.x + blockDim.x] = 0.0;
    __syncthreads();
    //without console stream
    unsigned int t = threadIdx.x;
    for(int stride = blockDim.x; stride >= 1; stride /= 2){
        if(t < stride)
            partialSum[t] += partialSum[t + stride];
            __syncthreads();
    } 
    //with console stream
    // unsigned int t = threadIdx.x * 2;
    // for(int stride = 1;stride <= blockDim.x;stride <<= 1){
    //     if(t % (2*stride) == 0){
    //         partialSum[t] +=partialSum[t + stride];
    //     }
    //     __syncthreads();
    // }
    if(t == 0){
        X[blockIdx.x] = partialSum[t];
    }
}
//host code
double reduceArray(double* array, unsigned int length){
    double *d_A;
    int size = length*sizeof(double);
    cudaMalloc(&d_A, size);
    cudaMemcpy(d_A, array, size, cudaMemcpyHostToDevice);
    int num_blocks = (length - 1)/THREAD_LENGTH/2 + 1;
    while(num_blocks >= 1){
        reduction_sum<<<num_blocks, THREAD_LENGTH>>>(d_A, length);
        if(num_blocks == 1)
            break;
        length = num_blocks; 
        num_blocks = (num_blocks - 1)/THREAD_LENGTH/2 + 1;
    }
    double result(0);
    cudaMemcpy(&result, d_A, sizeof(double), cudaMemcpyDeviceToHost);
    cudaFree(d_A);
    return result;
}
//test
int main(int argc, char **argv){
    int n = atoi(argv[1]);
    double *A = (double *)malloc(n * sizeof(double));
    for(int i = 0; i < n;++i){
        A[i] = 1.0;
    }
    double result = reduceArray(A, n);
    printf("%lf\n", result);
    free(A);
    return 0;
}

运行结果：

5.对于性能优化中的矩阵乘法 kernel 函数，以一个$16\times 16$的小型矩阵为例画出同一个 warp 中所有线程在第 9 行和第 10 行中的访问模式。计算同一个 warp 中每个线程的 tx 和 ty 的值，并在第 9 行和第 10行中计算 d_M 和 d_N 的索引值时使用 tx 和 ty 的值。说明在每次迭代过程中线程确实访问全局内存中连续的d_M和d_N位置。

假设线程块大小为 $16\times 16$，第一个warp就是包含列索引为 0 1 的两行thread。
d_M的索引：
 $(blockIdx.y\times TILE\_WIDTH+threadIdx.y)\times Width+m\times TILE\_WIDTH+threadIdx.x$
=$((blockIdx.y\times Width+m)\times TILE\_WIDTH)+threadIdx.y\times Width+threadIdx.x$

d_N的索引：
 $(m\times TILE\_WIDTH+threadIdx.y)\times Width+blockIdx.x\times TILE\_WIDTH+threadIdx.x$
=$((m\times Width+blockIdx.x)\times TILE\_WIDTH)+threadIdx.y\times +threadIdx.x$

括号内的值对于每次循环时单个block内的线程是相同的，所以只证明括号外的值：$threadIdx.y\times 16+threadIdx.x$产生了相邻的地址即可。

tx	ty	$threadIdx.y\times Width+threadIdx.x$
0	0	0
1	0	1
2	0	2
3	0	3
4	0	4
5	0	5
6	0	6
7	0	7
8	0	8
9	0	9
10	0	10
11	0	11
12	0	12
13	0	13
14	0	14
15	0	15
0	1	16
1	1	17
2	1	18
3	1	19
4	1	20
5	1	21
6	1	22
7	1	23
8	1	24
9	1	25
10	1	26
11	1	27
12	1	28
13	1	29
14	1	30
15	1	31

得证。

但所用线程块较小时，只是同一行访问是连续的，不同行的访问并不能合并。

10.对于下图的设计，编写矩阵乘法 kernel 函数。

#define BLOCK_SIZE 16
__global__ void
matrixMul(float *Pd, float *Md, float *Nd, int widthM, int widthN)
{
    unsigned int bx = blockIdx.x;
    unsigned int by = blockIdx.y;
    unsigned int tx = threadIdx.x;
    unsigned int ty = threadIdx.y;
    __shared__ float Ms[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ float Ns[BLOCK_SIZE][BLOCK_SIZE * 2];
    // M 的第一个子矩阵左上元素索引
    int mBegin = widthM * BLOCK_SIZE * by;
    // M 的最后一个子矩阵右上元素索引，用来控制循环何时结束
    int mEnd = mBegin + widthM - 1;
    // M 子矩阵每次循环需要加上的距离
    int mStep = BLOCK_SIZE;
    // N 的第一个子矩阵左上元素索引
    int nBegin = BLOCK_SIZE * bx;
    // N 子矩阵每次循环需要加上的距离
    int nStep = BLOCK_SIZE * widthN;
	//Psub 用来存储结果矩阵每块元素的部分和
    float Psub1 = 0.0f;
    float Psub2 = 0.0f;
    // 遍历所有子矩阵
    for (int m = mBegin, n = nBegin; m <= mEnd; m += mStep, n += nStep)
    {
        // 将数组加载到共享存储器，每个线程加载 3 个元素
        Ms[ty][tx] = Md[m + widthM * ty + tx];
        Ns[ty][tx] = Nd[n + widthN * ty + tx];
        Ns[ty][tx + blockDim.x] = Nd[n + widthN * ty + tx + blockDim.x];

        __syncthreads();
        // 一次计算两个矩阵
        for (int k = 0; k < BLOCK_SIZE; ++k)
        {
            Psub1 += Ms[ty][k] * Ns[k][tx];
            Psub2 += Ms[ty][k] * Ns[k][tx + blockDim.x];
        }

        __syncthreads();
    }
    // 写回计算结果，每个线程写一个
    int p = widthN * BLOCK_SIZE * by + BLOCK_SIZE * bx;
    p += widthN * ty + tx;
    Pd[p] = Psub1;
    Pd[p + blockDim.x] = Psub2;
}

12.一个年轻的工程师为了提高性能使用如下的kernel进行归约。（A）你认为性能会提升吗？（B）这个工程师应该受奖还是受罚？为啥？

__shared__ float partialSum[];
unsigned int tid = threadIdx.x;
for (unsigned int stride = n >> 1; stride >= 32; stride >>= 1)
{
    __syncthreads();
    if (tid < stride)
        shared[tid] += shared[tid + stride];
}
__syncthreads();
if (tid < 32)
{ // unroll last 5 predicated steps
    shared[tid] += shared[tid + 16];
    shared[tid] += shared[tid + 8];
    shared[tid] += shared[tid + 4];
    shared[tid] += shared[tid + 2];
    shared[tid] += shared[tid + 1];
}

会提升，会受奖赏

当执行到只有 32 个线程在计算时，每次循环都只有一个warp在使用，其他warp没有使用，会造成资源浪费。将最后一次计算时的warp单独展开，减少了要执行指令的数量，从而可以提高性能。