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\times16$ 的小型矩阵为例画出同一个 warp 中所有线程在第 9 行和第 10 行中的访问模式。计算同一个 warp 中每个线程的 tx 和 ty 的值，并在第 9 行和第 10行中计算 d_M 和 d_N 的索引值时使用 tx 和 ty 的值。说明在每次迭代过程中线程确实访问全局内存中连续的d_M和d_N位置。

假设线程块大小为 $16\times16$ ，第一个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$ 产生了相邻的地址即可。

得证。

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

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];
}

会提升，会受奖赏 在这里插入图片描述