CUDA(六) 周斌 CUDA 编程模型

GPU架构概览

1. GPU特别适用于：

•              密集计算，高度可并行计算
•              图形学

       2.晶体管主要被用于：

•                       执行计算
•                        而不是

                                         缓存数据

                                          控制指令流

CPU中ALU中占用比例不大，主要是做缓存和控制逻辑

GPU大量的面积用于做数据处理，对缓存和控制逻辑的部分不多。

 

GPU计算历史

2001/2002-GPU数据并行处理器

         GPGPU，通用GPU

2007 NVIDIA 发表CUDA

          CUDA Compute Uniform Device Architecture统一计算设备架构

         GPGPU从协处理器发展成完整的 进行并行计算的GPU Computing 

 

2008 opencL并行设备

 

CUDA的一些信息

1.   层次化线程集合  A hierarchy of thread groups
2.   共享存储Shared memories
3.  同步Barrier synchronization

CUDA术语

Host -即主机端，通常指CPU

         采用ANSI标准的C语言编程

Device - 即设备端 通常指GPU（数据可并行）

           采用ANSI标准C的扩展语言编程

Host 和 Device拥有各自的存储器，数据管理和数据存储分布设计

CUDA编程

            包括主机端和设备端两部分代码

Kernel-数据并行处理函数

               通过调用kernel函数在设备端创建轻量级线程

              线程由硬件负责创建并调度

类似于OpenGL的shader

 

CUDA 的和函数（kernels）

在N个不同的CUDA线程上并行执行

  通过编译器系统进行翻译，将整个代码翻译到GPU可以执行。

 

_global_ void VecAdd(float* A, float* B, float* C)
{
    int i = threadIdx.x;
    C[i] = A[i] + B[i];
}
int main()
{
    VecAdd<<1,N>>(A,B,C);
}

 

_global_ : Declaration Specifier修饰词

<<1,N>>:Execution Configuration 告诉指令系统需要多少个线程，多少个线程块的配置来启动核函数。

threadIdx.x: Thread ID线程ID标识在整个CUDA系统中线程在那个部分。

CUDA程序的执行

                                  CPU Serial Code    在CPU上串行执行的代码

                               GPU Parallel Kernel 调用到GPU核函数

                    KernelA<<<nBlk,nTid>>>(args);核函数在GPU上进行硬件的并行，返回结果到CPU端

    

                              CPU Serial Code

                               GPU Parallel Kernel

                    KernelA<<<nBlk,nTid>>>(args);

线程层次Thread Hierarchies

                   Grid - 一维 或 多维线程块（block）。多个线程块，便于索引

                              一维 或 二维

                  Block -  一组线程，线程块，便于索引

                             一维、二维或三维

                              一个Grid里面的每一个Block的线程数是一样的。

                             block内部的每个线程可以：

                                            同步访问 synchronize

                                            访问共享存储器 shared memory

grid--》block，有自己的编号-->thread，有自己的编号

A thread block is a batch of threads that can cooperate with each other by:

                                            Sychronizing their exectution

                                                             For hazard-free shared memory accesses

                                                 Efficiently sharing data through a low latency shared memory

Two threads from two different blocks cannot cooperate

编号有主意组织和处理特定维度的数据

            Block  - 一维，二维 或 三维。索引数组，矩阵，体。

                   Thread ID : Scalar thread identifier

                   线程索引：threadIdx

                 一维Block: Thread ID =Thread Index，线程属于Block，线程有编号

                二维Block(Dx,Dy): Thread ID of index(x,y)== x + yDy

                三维Block(Dx,Dy,Dz): Thread ID of index (x,y,z) == x+yDy+zDxDy

_global_ void MatAdd(float* A[N][N], float* B[N][N], float* C[N][N])
{
    int i = threadIdx.x;
    int j = threadIdx.y;
    C[i][j] = A[i][j] + B[i][j];
}
int main()
{
    int numBlocks = 1;
    dim3 threadsPerBlock(N,N); 
    MatAdd<<numBlocks ,threadsPerBlock>>(A,B,C);
}

numBlocks =1只使用一个block，1 Thread Block

threadsPerBlock block是二维索引 2D Block

 threadIdx.x  threadIdx.y 线程二维索引

 线程块Thread Block 

      线程的集合，线程的限制，每个线程块最多不超过

         G80和GT200：多达512个线程

         Fremi:多达1024个线程

      位于相同的处理器核心（SM）

      共享所在核的存储器

 

块索引 blockIdx

维度：blockDim

  一维、二维、三维

_global_ void MatAdd(float* A[N][N], float* B[N][N], float* C[N][N])
{
    int i = blockIdx.x*blockDim.x+threadIdx.x;
    int j = blockIdx.y*blockDim.y+threadIdx.y;
    if(i<N&&j<N)
        C[i][j] = A[i][j] + B[i][j];
}
int main()
{
    dim3 threadsPerBlock(16,16); 
    dim3 numBlocks(N/threadsPerBlock.x,N/threadsPerBlock.y);
    MatAdd<<numBlocks ,threadsPerBlock>>(A,B,C);
}

下标在全局的位置

16*16 threads per block

numBlocks :2D thread Block

对于每一个线程块的编号，都有下标的索引，

          例如：N=32

          每个块有16*16个线程（跟N无关）

           threadIdx([0,15],[0,15])

          Grid里面有2*2 个线程块Block

            blockIdx([0,1],[0,1])线程块的下标

            blockDim=16，线程块大小

                  int i = blockIdx.x*blockDim.x+threadIdx.x;
                  int j = blockIdx.y*blockDim.y+threadIdx.y;

             i=[0,1] *16 +[0,15]二维数据转换成1维索引。

线程索引对应组织架构理解

线程块之间彼此独立执行

       任意顺序：并行或串行（不是按照线程块的编号）

       被任意数量的处理器（SM）以任意顺序调度

                  处理器的数量具有可扩展性

  例如一个CUDA有8个Block放在某一个SM上执行，例如一个GPU有两个完整的SM，可以通过4次调度完成Block，如果有4个SM（完整的核，CUDA的核心是ALU计算单元），相同的程序只需要两次调度完成，每次执行4个。如果CUDA的Block数足够多，则可以利用SM资源。

一个块内部的线程

      共享容量有限的低延迟存储器（shared memory）

      同步执行

      合并访存

             __syncThreads()

            -barrier - 块内线程一起等待所有线程执行到某处语句

            -轻量级

 

CUDA 内存传输

   Device code can :GPU设备，不同部分存储器

          -R/W per-thread registers  读写每一个线程的私有寄存器

          -R/W per-thread local memory  

          -R/W per-block shared memory每一个线程块有一个公共的共享存储

         -R/W per-grid global memory  读写所有线程共享的显存上的global memory

         -Read only per-grid constant memory独立的存储空间，固定值的存储器，能够在多个线程在使用一个不太变化的内存，只能读

Host code can： CPU主机端代码

        -R/W per grid global and constant memories 读写global memory核constant memory

主机端Host 可以从device往返传输数据，通过PCIe总线

          Global memory 全局存储器（GPU中）

          Constant memory常量存储器（GPU中）

                      cudaMalloc()在设备端分配global memory 

                      cudaFree()释放存储空间

float *Md
int size = Width * Width * sizeof(float);
cudaMalloc((void**) &Md, size);
cudaFree(Md);

Md指针，返回设备端地址。指向了已经分配好的地址。指向设备端的指针 Pointer to device memory，只能在设备端使用，不能在CPU端使用。

size：size in bytes

cudaMemcpy()

        内存传输

                Host to host：Cpu内存上的数据可以通过GPU的方式拷贝

               Host to device

               Device to host

               Device to device

cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice);
cudaMemcpy(P, Pd, size, cudaMemcpyDeviceToHost);

cudaMemcpyHostToDevice:Host to device拷贝方向

Md:目的地址，一定要是指向GPU设备的指针 Destination(device)

M：原地址，在CPU主机端（Host）端地址.Source(host)

size：大小

 

Matrix Multiply 矩阵相乘算法提示

     向量

     点乘

      行优先或列优先

     每次点成结果输出一个元素

P=M*N假定M和N都是方阵

1000*1000矩阵

        1000000点乘

        each 1000 multiples and 1000 adds

void MatrixMulOnHost(float* M, float* N, float* P, int width)
{
    for (int i = 0; i < width; ++i)
        for (int j = 0; j < width; ++j)
        {
            float sum = 0;
            for (int k = 0; k < width; ++k)
            {
                folat a = M[i * width + k];
                float b = N[k * width + j];
                sum += a * b;
            }
            P[i * width + j] = sum;
        }
}

获取行和列，做累加，放到对应位置

int main(void){
1、    //Allocate and initialize the matrices M, N, P
       //I/O to read the input matrices M and N
2、    //M * N on the device
       MatrixMulOnDevice(M, N, P, width);

3、    //I/O to write the output matrix P
       //Free matrices M, N, P
return 0;
}

1、管理整个内存，为输入数据和输出结果分配空间。

2、数据放到GPU上之后，GPU并行计算。

3、将结果返回到CPU，释放内存。

Matrix Multiply 数据传输

第一步

       在算法框架中添加CUDA memory transfers  内存拷贝

void MatrixMulOnDevice(float* M, float* N, float* P, int Width)
{
    int size = Width * Width * sizeof(float);
1、//Load M and N to device memory
    cudaMalloc(Md, size);
    cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice);
    cudaMalloc(Nd, size);
    cudaMemcpy(Nd, N, size, cudaMemcpyHostToDevice);
  //Allocate P on the device
    cudaMalloc(Pd, size);

2、//Kernel invocation code - to be shown later
3、//Read P from the device
    cudaMemcpy(P, Pd, size, cudaMemcpyDeviceToHost);
  //Free device matrices
    cudaFree(Md); cudaFree(Nd); cudaFree(Pd);
}

核函数，处理的核心。

第二部

     CUDA C编程实现kernel

Matrix Multiply CUDA Kernel

//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;//访问一个matrix，所以采用二维block，每一个线程处理一个结果数据
    //Pvalue stores the Pd element that is computed by the thread
    float Pvalue = 0;//每一个Kernel线程计算1个输出
    for(int k = 0; k < Width; ++k)
    {    
        float Mdelement = Md[ty * Md.width + k];//行索引
        float Ndelement = Nd[k * Nd.width + tx];//方阵
        Pvalue += Mdelement * Ndelement;
    }
    // Write the matrix to device memory each thread writes one element
    Pd[ty * Width + tx] = Pvalue;
}

双层循环i,j使用并行处理替代了。

不需要锁和同步，因为没有依赖关系。每次计算结果都是独立的，不会用到其他计算结果。

Matrix Multiply ：调用Kernel

第三部

     CUDA C编程调用Kernel

    

//Setup the execution configuration
dim3.dimBlock(WIDTH, WIDTH);//一个block里面包含width*width个线程
dim3.dimGrid(1, 1);//每个Grid里面有几个block
//Launch the device computation threads
MatrixMulKernel<<dimGrid, dimBlock>>(Md, Nd, Pd);

Matrix Multiply ：样例

  一个线程block 计算Pd

     每个线程计算Pd的一个元素

每个线程

    读入矩阵Md的一行

    读入矩阵Nd的一行

    为每对Md和Nd元素执行一次乘法和加法

  （not very high）计算次数和片外访存次数比率接近1：1（不是很高）

 矩阵长度受限于一个线程快允许的线程数目

 

算法中最主要的性能问题是访存。