CUDA(八) 周斌 CUDA 编程三

matrix  multiply: CUDA  Kernel

void MatrixMulOnHost(float* M, float* N, float* P, int width)
{
    for(int i = 0; i < width; ++i)
        for(int j = 0; j < width; ++j)//双层循环用来访问结果矩阵P的所有元素
        {
            float sum = 0;
            for(int k = 0; k < width; ++k)//内层循环计算乘累加
            {
                float a = M[i * width + k];//i行*width行宽 + k列 (i, k)
                float b = N[k * width + j];//k行*width行宽 + j列 (k, j)
                sum += a * b;  
            }
            P[i * width + j] = sum;//第i行的所有列，第j列对应列的所有行。
        }
}

//Matrix multiplication kernel - thread specification
__global—— void MatrixMulKernel(float* Md, float* Nd, float* Pd, int Width)
{
    //2D Thread ID
    int tx = threadIdx.x;
    int ty = threadIdx.y;//访问矩阵，因此用二维block（使用了一个block，有x、y方向的线程，每个线程负责一个结果元素）
    
    //Pvalue stores the Pd element that is computed by the thread
    float Pvalue = 0;//每个Kernel线程计算出一个输出结果
    
    for(int k = 0; k < Width; ++k)//CPU版本的2个外层循环被拿掉了
    {
        float Mdelement = Md[ty * Md.width + k];//对global memory的访问，ty行数*每行个数+当前列
        float Ndelement = Nd[k * Nd.width + tx];//读了两个global memory的数据，k行数*每行个数+当前列tx
        Pvalue += Mdelement * Ndelement;//才进行一次乘累加
    }
    
    //Write the matrix to device memory each thread writes one element
    Pd[ty * Width + tx] = Pvalue;//由于计算矩阵的每一个输出结果不依赖于其他的一个结果，不需要锁或者同步。
    //最后还有一次写global memory的操作，所以加一起global memory开销巨大
}

矩阵长度被限制，因为只是用一个block。例如G80和GT200最多512个线程/block，不能计算矩阵大小大于512矩阵。

希望GPU能够处理巨大规模的矩阵运算。

去除长度的限制

       将Pd（结果矩阵）矩阵拆成tile（瓦片）小块，之前是一个block负责一个大的完整矩阵。将一个结果矩阵同个多个block的形式切分。

       把一个tile布置到一个block

      通过threadIdx和blockIdx索引

矩阵相乘中有很多global memory的访问

 

例如

     结果矩阵4*4

   TITLE_WIDTH = 2

   Block尺寸：2*2

   结果矩阵被分成了4块，Block(0,0),Block(1,0),Block(0,1),Block(1,1)

   每个block取出输入矩阵一整行乘上输入矩阵一整列，block线程数有限，但是block数是很多的，

__global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, int Width)
{
    int Row = blockIdx.y * blockDim.y + threadIdx.y;//哪一个block(blockIdx.y)*每一个block有多少个线程(blockDim.y)+在block中的线程号,x和y方向的索引，计算Pd（结果矩阵）和M的行索引
    int Col = blockIdx.x * blockDim.x + threadIdx.x;

    float Pvalue = 0;
    for (int k = 0; k < Width; ++k)
        Pvalue += Md[Row * Width +k] * Nd[k * Width + Col];//做一个乘累加，每一个线程计算块内子矩阵的一个元素
    Pd[Row * Width + Col] = Pvalue;//将计算结果写回结果矩阵
}

调用kernel

dim3 dimGrid(Width / TIE_WIDTH, HEIGHT / TILE_WIDTH);//定义有多少个block，计算宽度和高度方向上需要多少个block
dim3 dimBlock(TILE_WIDTH, TILE_WIDTH);//定义block的维度，dmi3是一个三维的索引类型（例如x, y, z），如果只用前两个，就可以用作二维，处理一个小块，建立的是一个tile的宽度，也就是每个block里面有多少个线程

MatrixMulKernel<<<dimGrid, dimBlock>>>(Md, Nd, Pd, TILE_WIDTH);//dimGrid表示grid里面多少个block,dimBlock表示每个block里面线程数是怎样的。Grid维度和Block维度，参数和tile宽度。

如果不能整除，就补零。设计问题最好设成可以整出的Width / TIE_WIDTH。

global memory 读写没做一次乘累加需要读取两次数据，并且每行，每一列要反复地读。global memory读写很慢

受限于global memory带宽

     G80峰值GFLOPS:346.5GFLOPS，每秒的浮点数运算

    需要1386 GB/s 的带宽来达到：所以每秒钟需要提供的浮点数为1386GB=346.5*4（每一个浮点数的字节数）

    G80存储器实际带宽：86.4GB/s：实际的带宽

               限制代码21.6GFLOPS：86.4/4=21.6GFLOPS  1/15的性能

              实际上，代码运行速度是15GFLOPS因为还有其他额外的开销，峰值性能的1/20

     必须大幅减小对global memory的访问。

每一个输入元素被Width个线程读取，反复读取，例如对每一列元素多读了row number次，每一行也是一样的

使用shared memory（内部共享）来减少global memory带宽需求。

 

把kernel拆分成多个阶段

   每个阶段用Md和Nd的子集累加Pd，每个block负责其中的一个子块，需要M矩阵的一行和N矩阵的一列数据乘累加，才能得到block中P的结果，可以将M矩阵的一行和N矩阵的一列数据全部都加载到block的shared memory里面，在shared memory里面做乘累加，获得block结果。

   每个阶段有很好的数据局部性，可是每一行和每一列都比较大，因此可以把每一行和每一列数据在做一次切分。

每个线程

     读入瓦片tile内的Md和Nd的一个元素      

    存入shared memory

__global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, int Width)
{
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];//为shared memory建立两个临时的存储空间，所有的线程都在执行这段代码，shared memory是对block而言的，
    __shared__ float Ndx[TILE_WIDTH][TILE_WIDTH];//通过shared memory 建立Md和Nd这样的子集

    int bx = blockIdx.x;    int by = blockIdx.y;
    int tx = threadIdx.x;   int ty = threadIdx.y;//获取以前的原始索引，与以前的block是一样的依旧是一个block负责其中一个子块

    int Row = by * TILE_WIDTH + ty;
    int Col = TILE_WIDTH + tx;
    
    float Pvalue = 0;
    for (int m = 0; m < Width/TILE_WIDTH; ++m)//多了一个外层循环，因为乘累加只在shared memory里面的一小块，以shared memory大小作为循环步长。Width/TILE个阶段数目,m是阶段索引
    {
        Mds[ty][tx] = Md[Row*Width + (m*TILE_WIDTH + tx)];
        Nds[ty][tx] = Nd[Col + (m*TILE_WIDTH + ty)*Width];//从M和N里面每一个线程去其中一块数据放到shared memory里面去。目的是为了后面做计算的时候大家都可以用。
        __syncthreads();//等待block内所有线程，即等到整个瓦片存入shared memory

        for (int k = 0; k < TILE_WIDTH; ++k)
            Pvalue += Mds[ty][k] * Nds[k][tx];//计算完成之后其实还是没有算完，只是算了其中一个tile，乘累加子集
        __syncthreads();//同步语句，如果没有同步语句，则Mds和Nds内的部分已经更改
    }//必须在外层循环都做完的时候，才能得到最后的结果
    Pd[Row*Width+Col] = Pvalue;//结果写入global memory
}

TILE_WIDTH大小的选择

    如果太大

           超出一个块允许的最大线程数。

                     G80 and GT200 - 512

                    Fermi - 1024

                    Kerpler - 1025

                   可疑的，依据不同的计算能力，请查表

         超出shared memory 极限

                 G80: 16KB/SM 并且8 blocks / SM

                          2KB / block

                          1 KB 给Nds， 1KB 给Mds（16 * 16 * 4）一个浮点4个字节

                          TILE_WIDTH = 16

                          更大的TILE_WIDTH将导致更少的块数

Shared memory 瓦片化的好处

            global memory 访问次数减少TILE_WIDTH倍

                                 16 * 16 瓦片tile 减少16倍

           G80

                  现在global memory 支持345.6GFLOPS  20*16

                接近峰值346.5GFLOPS

G80线程尺寸的考虑

   每个thread block 有许多个线程

         TILE_WIDTH 为16时， 16*16=256个线程（16k shared memory，8个block，每个block 2k,每个M、N=1k， 1k/4=256个）

需要许多个thread blocks

       一个1024*1024 Pd需要：64*64 = 4K Thread Blocks =4k*16=1024

每个 thread block执行2*256 = 512次global memory的float读入，为了提供256*（2*16）= 8k mul/add 操作。每个线程256个16乘累加M、N

      存储带宽不再是限制因素

 

Atomic Functions原子函数

许多原子操作：GPU为了使用特定操作，特定操作只有一个线程操作完成后，第二个线程才能进行，原子操作对于同步和数据全局共享方法。原子操作尽量少用，很耗时，需要在系统里专门开辟完整的排队

//算术运算

atomicAdd();

atomicSub();

atomicMin();

atomicExch();

atomicMin();

atomicMax();

atomicAdd();

atomicDec();

atomicCAS();

//位运算

atomicAnd();

atomicOr();

atomicXor();

 

不同块里面的线程协作。

尽量少用原子操作。针对同一个内存地址做的原子操作加。意味着所有线程都需要排队，一旦进行内存的排队，访问开销巨大，非必要少使用。