精华内容
下载资源
问答
  • CUDA学习笔记

    2019-02-18 14:34:54
    CUDA学习笔记 目录 CUDA学习笔记 函数类型限定符 __global__ __host__ __device__ 变量类型限定符 __device__ __shared__ __constant__ CUDA线程概念 网格(Grid)、线程块(Block)和线程(Thread)的...

    CUDA学习笔记

    目录

    CUDA学习笔记

    函数类型限定符

    __global__

    __host__

    __device__

    变量类型限定符

    __device__

    __shared__

    __constant__

    CUDA线程概念

    网格(Grid)、线程块(Block)和线程(Thread)的组织关系

    线程索引的计算公式

    dim3结构类型


    函数类型限定符

    __global__

    表示一个内核函数,主机(CPU)上调用,但在设备(GPU)上执行

    __host__

    由CPU调用,且在CPU执行的函数,默认的函数限定符

    __device__

    设备执行的程序,device前缀定义的函数只能在GPU上执行,所以device修饰的函数里面不能调用一般常见的函数


    变量类型限定符

    __device__

     

    __shared__

     

    __constant__

     


    CUDA线程概念

    CUDA线程的概念与CPU的概念基本一致:进程是系统资源分配的基本单位, 而线程是CPU 能调度和独立运行的基本单位(最小单元)


    网格(Grid)、线程块(Block)和线程(Thread)的组织关系

    https://blog.csdn.net/gqixf/article/details/80760701 

    多个线程thread组成一个block,多个block组成一个grid.

    CUDA的软件架构由网格(Grid)、线程块(Block)和线程(Thread)组成,相当于把GPU上的计算单元分为若干(2~3)个网格,每个网格内包含若干(65535)个线程块,每个线程块包含若干(512)个线程,三者的关系如下图:

    Thread,block,grid是CUDA编程上的概念,为了方便程序员软件设计,组织线程。

    • thread:一个CUDA的并行程序会被以许多个threads来执行
    • block:数个threads会被群组成一个block,同一个block中的threads可以同步,也可以通过shared memory通信 
    • grid:多个blocks则会再构成grid

    线程索引的计算公式

    一个Grid可以包含多个Blocks,Blocks的组织方式可以是一维的,二维或者三维的。block包含多个Threads,这些Threads的组织方式也可以是一维,二维或者三维的。

    CUDA中每一个线程都有一个唯一的标识ID—ThreadIdx,这个ID随着Grid和Block的划分方式的不同而变化,这里给出Grid和Block不同划分方式下线程索引ID的计算公式。

    1、 grid划分成1维,block划分为1维
        int threadId = blockIdx.x *blockDim.x + threadIdx.x;  
    2、 grid划分成1维,block划分为2维  
        int threadId = blockIdx.x * blockDim.x * blockDim.y+ threadIdx.y * blockDim.x + threadIdx.x;  
    3、 grid划分成1维,block划分为3维  
        int threadId = blockIdx.x * blockDim.x * blockDim.y * blockDim.z  
                           + threadIdx.z * blockDim.y * blockDim.x  
                           + threadIdx.y * blockDim.x + threadIdx.x;  
    4、 grid划分成2维,block划分为1维  
        int blockId = blockIdx.y * gridDim.x + blockIdx.x;  
        int threadId = blockId * blockDim.x + threadIdx.x;  
    5、 grid划分成2维,block划分为2维 
        int blockId = blockIdx.x + blockIdx.y * gridDim.x;  
        int threadId = blockId * (blockDim.x * blockDim.y)  
                           + (threadIdx.y * blockDim.x) + threadIdx.x;  
    6、 grid划分成2维,block划分为3维
        int blockId = blockIdx.x + blockIdx.y * gridDim.x;  
        int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)  
                           + (threadIdx.z * (blockDim.x * blockDim.y))  
                           + (threadIdx.y * blockDim.x) + threadIdx.x;  
    7、 grid划分成3维,block划分为1维 
        int blockId = blockIdx.x + blockIdx.y * gridDim.x  
                         + gridDim.x * gridDim.y * blockIdx.z;  
        int threadId = blockId * blockDim.x + threadIdx.x;  
    8、 grid划分成3维,block划分为2维  
        int blockId = blockIdx.x + blockIdx.y * gridDim.x  
                         + gridDim.x * gridDim.y * blockIdx.z;  
        int threadId = blockId * (blockDim.x * blockDim.y)  
                           + (threadIdx.y * blockDim.x) + threadIdx.x;  
    9、 grid划分成3维,block划分为3维
        int blockId = blockIdx.x + blockIdx.y * gridDim.x  
                         + gridDim.x * gridDim.y * blockIdx.z;  
        int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z)  
                           + (threadIdx.z * (blockDim.x * blockDim.y))  
                           + (threadIdx.y * blockDim.x) + threadIdx.x;    


    dim3数据结构

    dim3是CUDA 中一个比较特殊的数据结构,我们可以用这个数据结构创建一个二维的线程块与线程网格。例如在长方形布局的方式中,每个线程块的X 轴方向上开启了32 个线程, Y轴方向上开启了4 个线程。在线程网格上, X 轴方向上有1 个线程块, Y 轴方向有4 个线程块。

    dim3 threads_rect(32,4)
    dim3 blocks_rect(1,4)

    dim3    blocks(DIM/16,DIM/16);
    dim3    threads(16,16);
    kernel<<<blocks,threads>>>( d->dev_bitmap, ticks );

     

     

     

     

     

    展开全文
  • cuda学习笔记

    2017-02-19 15:13:29
    cuda学习笔记标签(空格分隔): 学习笔记一、学习平台1.1学习平台搭建与学习平台基本信息1.1.1学习平台搭建此篇文档的整理基于nvidia公司出品的GeForce GTX 950 GPU,在电脑主机当中安装好独立显卡之后,安装cuda...

    cuda学习笔记

    标签(空格分隔): 学习笔记


    一、学习平台

    1.1学习平台搭建与学习平台基本信息

    1.1.1学习平台搭建

    此篇文档的整理基于nvidia公司出品的GeForce GTX 950 GPU,在电脑主机当中安装好独立显卡之后,安装cuda7.0至软件盘(不用再单独安装显卡驱动程序)。在vs下新建cuda工程,就可以编写cuda程序了。

    1.1.2学习平台基本信息

    在编写cuda程序时,程序的头文件应该包括 “cuda_runtime.h”和”device_launch_parameters.h”;以下一段代码用来查看显卡gpu的计算性能和架构

    int main()
    {
        cudaDeviceProp prop;
    
        int count;
        cudaGetDeviceCount(&count);
    
        for (int i = 0; i < count; ++i){
            cudaGetDeviceProperties(&prop, i);
            printf("     --- Genaral Information for Device %d ---\n", i);
            printf("Name : %s\n",prop.name);
            printf("Compute capability : %d.%d\n", prop.major, prop.minor);
            printf("Clock rate:  %d\n",prop.clockRate);
            printf("Device copy overlap: ");
    
            if (prop.deviceOverlap){
                printf("Enabled\n");
            }
            else{
                printf("Disabled\n");
            }
    
            printf("Kernel execition timeout : ");
            if (prop.kernelExecTimeoutEnabled)
                printf("Enabled\n");
            else
                printf("Disabled\n");
    
            printf("  ---Memory Information for Device %d ---\n",i);
            printf("Total global mem: %ld\n",prop.totalGlobalMem);
            printf("Total const mem : %ld\n", prop.totalConstMem);
            printf("Max mem pitch : %ld\n", prop.memPitch);
            printf("Texture Alignment : %ld\n",prop.textureAlignment);
    
    
            printf("   ---MP Information for device %d ---\n", i);
            printf("Multiprocessor count : %d\n", prop.multiProcessorCount);
            printf("shared mem per mp: %d\n", prop.sharedMemPerBlock);
            printf("Register per mp: %d\n",prop.regsPerBlock);
            printf("Threads in warp: %d\n", prop.warpSize);
            printf(" Max threads per block :%d\n", prop.maxThreadsPerBlock);
            printf("Max thread dimentions : (%d, %d, %d)\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
            printf("Max grid dimensions:(%d, %d, %d)\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2] );
    
            printf("\n\n\n");
    
        }
    }

    950显卡运行结果如下图所示:

    此处输入图片的描述

    二、cuda线程模型

    2.1线程索引相关问题

    我们把在GPU上启动的线程块集合称为一个线程格。从名字的含义可以看出,线程格既可以使一维的线程块集合,也可以是二维的线程块集合。核函数的每个副本都可以通过内置变量blockIdx来判断哪个线程块正在执行它。同样,它还可以通过内置变量gridDim来获得线程块的大小。通过这两个变量来计算每个线程块需要的数据索引。
    则当前线程块的索引如下:
    线程块的索引=行索引*线程格的数目+列索引
    blockIdx.y * gridDim.x+blockIdx.x;
    同样的:
    线程索引 = 行索引*线程块的数目+列索引
    threadIdx.y * blockIdx.x + threadIdx.x;
    int offset = x + y * Dim;在这里Dim表示线程块的大小(也就是线程的数量),y为线程块索引,并且x为线程块中的线程索引,所以计算得到如下索引:

    int tid = blockDim.x*blockIdx.x + threadIdx.x;
    tid += blockDim.x*gridDim.x;
    //每个线程块中的数量乘以线程格中线程块的总数量,即为当前线程格中运行的线程总数量。

    对于二维线程的索引,有如下代码:

    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    2.2线程开辟

    一种解决方案是将线程块的大小设置为某个固定数值BLOCKSIZE,然后启动N/BLOCKSIZE个线程块,这样就相当启动了N个线程同时运行。通常我们设置的线程块的个数为(N+BLOCKSIZE-1)/BLOCKSIZE来防止0线程的开辟问题。

    2.3多维数据线程调用

    线性存储器也可以通过cudaMallocPitch()和cudaMalloc3D分配。在分配二维和三维数组的时候,推荐使用,因为上述调用保证了GPU的最佳性能。返回的(pitch,stride)必须用于访问数组元素。下面的代码分配了一个尺寸为weight*height的二维浮点数组,同时演示了怎么在设备代码中遍历数组元素

    //host code
    int width =64,height = 64;
    float *dexPtr;
    int pitch;
    cudaMallocPitch((void **)&devPtr,&pitch,width*sizeof(float),height);
    kernel<<<100,512>>>(devPtr,pitch,widtf,height);
    
    //device code
    __global__void kernel(float* devPtr,int pitch,int width,int height)
    {
    for(int i =0;i<helght;++i){
    float* row =(float*) ((char*)devPtr+i*pitch);
    for(int j = 0;j<width;++j){
    float element = row[i];
            }
        }
    }

    下面的代码演示分配一个尺寸为width*height*depth的三维浮点数组,同时演示了怎么在设备代码中遍历数组元素。

    //host code
    cudaPitchedPtr devPitchedPtr;
    cudaExtent extent = make_cudaExtent(64,64,64);
    cudaMalloc3D(&devPitchedPtr,extent);
    kernel<<<100,512>>>(devPitchedPtr,extent);
    
    //device code
    __global__ void kernel(cudaPitchedPtr devPitchedPtr,cudaExtent extent){
    char *devPtr=devPitchedPtr.ptr;
    size_t pitch =devPitchedPtr.pitch;
    size_t slicePitch = pitch*extent.height;
    
    for(int i=0;i<extent.depth;++i){
    char *slice = devPtr + i*slicePitch;
    for(int j=0;j<extent.height;++j){
    float *row=(float*)(slice + y*pitch);
    for(int x = 0;x<extent.width;++x){
    float element = row[x];
            }
        }
    }

    三、矩阵乘法

    3.1使用全局内存

    #include<iostream>
    #include<fstream>
    using namespace std;
    #define AROWS 16
    #define ACOLS 4
    #define BROWS 4
    #define BCOLS 16
    void randomInit(float* _data, int _size)
    {
        for (int i = 0; i < _size; ++i)
        {
            _data[i] = rand() %5;
            //_data[i] = (float)(rand() / (float)RAND_MAX);
        }
    }
    __global__ void MatMul(float *dev_a, float *dev_b, float *dev_c)
    {
    
        int col = threadIdx.x;
        int row = threadIdx.y;
        float cvalue = 0;
        for (int i = 0; i<ACOLS; i++)
        {
            cvalue += dev_a[row*ACOLS+i] * dev_b[BCOLS*i + col];
        }
        dev_c[row*BCOLS + col] = cvalue;
    }
    
    int main()
    {
        float *a, *b, *c, *dev_a, *dev_b, *dev_c;
        a = (float*)malloc(sizeof(float)*AROWS*ACOLS);
        b = (float*)malloc(sizeof(float)*BROWS*BCOLS);
        c = (float*)malloc(sizeof(float)*AROWS*BCOLS);
        memset(c, 0, sizeof(float)*AROWS*BCOLS);
        cudaMalloc((void **)&dev_a, sizeof(float)*AROWS*ACOLS);
        cudaMalloc((void **)&dev_b, sizeof(float)*BROWS*BCOLS);
        cudaMalloc((void **)&dev_c, sizeof(float)*AROWS*BCOLS);
        cudaMemset(dev_c, 0, sizeof(float)*AROWS*BCOLS);
        //给矩阵AB赋初值
        randomInit(a, AROWS*ACOLS);
        randomInit(b, BROWS*BCOLS);
        for (int i = 0; i < AROWS*ACOLS; i++)
        {
            cout << "a[" << i << "]=" << a[i] << endl;
        }
        for (int i = 0; i < BROWS*BCOLS; i++)
        {
            cout << "b[" << i << "]=" << b[i] << endl;
        }
    
        cudaMemcpy(dev_a, a, sizeof(float)*AROWS*ACOLS, cudaMemcpyHostToDevice);
        cudaMemcpy(dev_b, b, sizeof(float)*BROWS*BCOLS, cudaMemcpyHostToDevice);
        dim3 blocks(16, 16);
        dim3 threads(AROWS+15, BCOLS+15);
        MatMul << <blocks, threads >> >(dev_a, dev_b, dev_c);
        cudaMemcpy(c, dev_c, sizeof(float)*AROWS*BCOLS, cudaMemcpyDeviceToHost);
    
        //打印C矩阵到文件中
        ofstream outFile;
        outFile.open("result.txt");
        outFile << fixed;
        outFile.precision(2);
        outFile.setf(ios_base::showpoint);
        for (int i = 0; i < AROWS*BCOLS; i++)
        {
            outFile << "C[" << i << "]=" << c[i] << endl;
        }
    
        system("pause");
        cudaFree(dev_a);
        cudaFree(dev_b);
        cudaFree(dev_c);
        free(a);
        free(b);
        free(c);
        return 0;
    }

    使用全局内存实现矩阵的一维向量乘法如上程序所示,但这种实现方式并没有充分利用gpu的优势。下面的代码是使用共享内存实现矩阵乘法。

    3.2使用共享内存

    四、关于cuda的基本架构

    在CUDA架构下,线程的最小单元是thread,多个thread组成一个block,多个block再组成一个grid。每一个block中开辟的所有线程共享同一个shared memory。block里面的thread之间的通信和同步所带来的开销是比较大的。SM以 32 个 Thread 为一组的 Warp 来执行 Thread。Warp内的线程是静态的,即在属于同一个warp内的thread之间进行通信,不需要进行栅栏同步(barrier)。

    4.1共享内存

    每个block中开辟的所有线程共享一个shared memory。使用共享内存变量的时候,需要在声明的时候加上shared关键词修饰,使用共享内存的时候应注意同一个线程块中的线程都执行结束才能进行下一步操作,所以需要使用__syncthread关键词使得block中的线程同步。

    4.2常量内存

    常量内存用于保存在核函数执行期间不会发生变化的数据,定义常量内存的时候应该使用关键词constant进行修饰。当从主机内存复制到GPU上的常量内存时,需要使用cudaMemcpyToSymbol()复制数据。常量内存读取数据可以节约内存带宽。注:只要当一个warp的半线程束中的所有16个线程有相同的读取请求时,才值得使用常量内存。

    4.3纹理内存

    纹理内存是专门为那些在内存访问模式中存在大量空间局部性的图形应用程序而设计的。
    以下内容参考博文http://www.cnblogs.com/traceorigin/archive/2013/04/11/3015755.html

    4.3.1、概述

    纹理存储器中的数据以一维、二维或者三维数组的形式存储在显存中,可以通过缓存加速访问,并且可以声明大小比常数存储器要大的多。
    在kernel中访问纹理存储器的操作称为纹理拾取(texture fetching)。将显存中的数据与纹理参照系关联的操作,称为将数据与纹理绑定(texture binding).
    显存中可以绑定到纹理的数据有两种,分别是普通的线性存储器和cuda数组。
    注:线性存储器只能与一维或二维纹理绑定,采用整型纹理拾取坐标,坐标值与数据在存储器中的位置相同;
    CUDA数组可以与一维、二维、三维纹理绑定,纹理拾取坐标为归一化或者非归一化的浮点型,并且支持许多特殊功能。

    4.3.2、纹理缓存:

    (1)、纹理缓存中的数据可以被重复利用
    (2)、纹理缓存一次预取拾取坐标对应位置附近的几个象元,可以实现滤波模式。

    4.3.3纹理存储器的使用:

    使用纹理存储器时,首先要在主机端声明要绑定到纹理的线性存储器或CUDA数组
    (1)声明纹理参考系

    texture<Type, Dim, ReadMode> texRef;
    //Type指定数据类型,特别注意:不支持3元组
    //Dim指定纹理参考系的维度,默认为1
    //ReadMode可以是cudaReadModelNormalizedFloat或cudaReadModelElementType(默认)

    注:纹理参照系必须定义在所有函数体外
    (2) 声明CUDA数组,分配空间
    CUDA数组可以通过cudaMalloc3DArray()或者cudaMallocArray()函数分配。前者可以分配1D、2D、3D的数组,后者一般用于分配2D的CUDA数组。使用完毕,要用cudaFreeArray()函数释放显存。

     //1数组
     cudaMalloc((void**)&dev_A, data_size);
     cudaMemcpy(dev_A, host_A, data_size, cudaMemcpyHostToDevice);
     cudaFree(dev_A);
    
     //2维数组
     cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>()
     cudaArray *cuArray;
     cudaMallocArray(&cuArray, &channelDesc, 64, 32); //64x32
     cudaMemcpyToArray(cuArray, 0, 0, h_data, sizeof(float)*width*height, cudaMemcpyHostToDevice);
     cudaFreeArray(cuArray);
    
     //3维数组 64x32x16
     cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<uchar>();
     cudaArray *d_volumeArray;
     cudaMalloc3DArray(&d_volumeArray, &channelDesc, volumSize);
    
     cudaMemcpy3DParms copyParams = {0};
     copyParams.srcPtr   = make_cudaPitchedPtr((void*)h_volume, volumeSize.width*sizeof(uchar), volumeSize.width, volumeSize.height);
     copyParams.dstArray = d_volumeArray;
     copyParams.extent   = volumeSize;
     copyParams.kind     = cudaMemcpyHostToDevice;
     cudaMemcpy3D(&copyParams);
    
     tex.normalized = true;
     tex.filterMode = cudaFilterModeLinear;
     tex.addressMode[0] = cudaAddressModeWrap;
     tex.addressMode[1] = cudaAddressModeWrap;
     tex.addressMode[2] = cudaAddressModeWrap;

    (3)设置运行时纹理参照系属性

    struct textureReference
    {
        int normalized;
        enum cudaTextureFilterMode filterMode;
        enum cudaTextureAddressMode addressMode[3];
        struct cudaChannelFormatDesc channelDesc;
    }

    normalized设置是否对纹理坐标归一化
    filterMode用于设置纹理的滤波模式
    addressMode说明了寻址方式

    (4)纹理绑定
    通过cudaBindTexture() 或 cudaBindTextureToArray()将数据与纹理绑定。
    通过cudaUnbindTexture()用于解除纹理参照系的绑定
    注:与纹理绑定的数据的类型必须与声明纹理参照系时的参数匹配
    (I).cudaBindTexture() //将1维线性内存绑定到1维纹理

    cudaError_t cudaBindTexture
    (    
        size_t *     offset,
        const struct textureReference *     texref,
        const void *     devPtr,
        const struct cudaChannelFormatDesc *     desc,
        size_t     size = UINT_MAX     
    )

    (II).cudaBindTexture2D //将1维线性内存绑定到2维纹理

    cudaError_t cudaBindTexture2D(    
        size_t *     offset,
        const struct textureReference *     texref,
        const void *     devPtr,
        const struct cudaChannelFormatDesc *     desc,
        size_t     width,
        size_t     height,
        size_t     pitch     
    )

    (III). cudaBindTextureToArray() //将cuda数组绑定到纹理

     cudaError_t cudaBindTextureToArray   
     (    
     const struct textureReference *     texref,
     const struct cudaArray *     array,
     const struct cudaChannelFormatDesc *     desc     
     )

    (5)纹理拾取
      对于线性存储器绑定的纹理,使用tex1Dfetch()访问,采用的纹理坐标是整型。由cudaMallocPitch() 或者 cudaMalloc3D()分配的线性空间实际上仍然是经过填充、对齐的一维线性空   间,因此也用tex1Dfetch()
      对与一维、二维、三维cuda数组绑定的纹理,分别使用tex1D(), tex2D() 和 tex3D()函数访问,并且使用浮点型纹理坐标。

    4.4cuda流

    五、使用事件对cuda程序进行计时

    创建事件

    cudaEvent_t  start,stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start,0);

    ……..此处省略cuda程序

    事件结束

    cudaEventRecord(stop,0);
    cudaEventSynchronize(stop);
    float elapsedTime;
    cudaEventElapsedTime(&elapsedTime,start,stop);
    printf("Time to generate :  %3.1f ms\n",elapsedTime);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);

    六、调试cudaC

    Parallel Nsight, visual profiler

    4.1线程最优配置tips
    最优的cuda线程配置
    1 每个SM上面失少要有192个激活线程,寄存器写后读的数据依赖才能被掩盖

    2 将 寄存器 的bank冲突降到最低,应尽量使每个block含有的线程数是64的倍数

    3 block的数量应设置得令可用的计算资源被充分的利用。由于每个block映射到一个sm上面,所以至少应该让block的数目跟sm的数目一样多。

    4 当Block中的线程被同步时或者等待读取设备存储器时,相应的SM会闲置。通常让block的数目是sm的2倍以上,使其在时间轴上重叠

    5 如果block的数目足够多,则每个Block里的线程数应设置成warp尺寸的整数倍,以免过小的warp浪费计算资源。

    6 给每个block分配越多的线程,能更高效的让他们在时间片上工作。但是相应的每个线程的寄存器也就越少。当寄存器过少,有可能因为访问溢出的寄存器,而导致数据的存储变慢。

    7 当每个线程占用的寄存器较多时,不宜在Block内分配过多的线程,否则也会减少block的数目。从而使SM的工作效率降低

    8 每个block内的线程数应遵循 相应的 计算能力等级中的规定数目。

    9 当线程块的数量为GPU中处理器数量的2倍时,计算性能达到最优。

    展开全文
  • CUDA 学习笔记

    2020-07-12 16:15:23
    书籍:《CUDA高性能并行计算》 windows10 VS2019 CUDA11.0 1.主流的标准操作:创建一些额外的数组并显式地使用cudaMemcpy()函数,在主机和设备间通过PCIe总线进行数据的输入和输出。 #ifndef distance_h #define ...

    windows10 VS2019 CUDA11.0
    1.主流的标准操作:创建一些额外的数组并显式地使用cudaMemcpy()函数,在主机和设备间通过PCIe总线进行数据的输入和输出。

    //程序改编自书籍:《CUDA高性能并行计算》
    #ifndef distance_h
    #define distance_h
    
    void distanceArray(float* out, float* in, float ref, int n);
    float scale(int i, int n);
    
    #endif
    
    #include <stdio.h>
    //也可不加此头文件,因为CUDA内联文件已经包含了math.h,程序可以正常运行
    //只是不加此头文件,VS会因为识别不了sqrt()函数而报错
    #include <cmath>
    //也可不加以下头文件,程序可以正常运行,只是VS会报错
    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    #define TPB 32
    
    //从核函数中调用并在GPU上执行,因此需要添加__device__标识。
    __device__ float distance(float x1, float x2)//左右各两个下划线
    {
    	return sqrt((x2 - x1) * (x2 - x1));
    }
    
    //如所有的核函数一样,是从主机上调用而在设备上执行,因此需要使用修饰符__global__,并且设置返回类型为void。
    __global__ void distanceKernal(float* d_out, float* d_in, float ref)
    {
    	const int i = blockIdx.x * blockDim.x + threadIdx.x;//计算出线程索引号,然后分散到各线程里进行并行计算
    	const float x = d_in[i];
    	d_out[i] = distance(x, ref);
    	printf("i=%2d: dist from %f to %f is %f.\n", i, ref, x, d_out[i]);
    }
    
    //以下为主流的标准操作:创建一些额外的数组并显式地使用cudaMemcpy()函数,在主机和设备间通过PCIe总线进行数据的输入和输出。
    
    //函数distanceArray()调用核函数distanceKernel()。这样的函数被称为封装(warpper)或者启动函数(launcher)。
    __host__ void distanceArray(float* out, float* in, float ref, int n)
    {
    	//声明指向设备端数组的指针
    	float* d_in = NULL;
    	float* d_out = NULL;
    	//在设备端为设备端数组分配内存
    	cudaMalloc(&d_in, n * sizeof(float));
    	cudaMalloc(&d_out, n * sizeof(float));
    	//把主机端上的数组复制到设备端对应的数组中
    	cudaMemcpy(d_in, in, n * sizeof(float), cudaMemcpyHostToDevice);
    	//启动核函数进行相关的计算
    	distanceKernal << <n / TPB, TPB >> > (d_out, d_in, ref);
    	//把设备端的输出数组复制到主机端对应的数组中
    	cudaMemcpy(out, d_out, n * sizeof(float), cudaMemcpyDeviceToHost);
    	//释放分配给设备端数组的内存
    	cudaFree(d_in);
    	cudaFree(d_out);
    }
    
    //在CPU执行的for循环内调用,来生成输入数据。合适的修饰符是__host__,但是我们可以不添加它,因为它是默认的标识符。
    __host__ float scale(int i, int n)
    {
    	return (float)i / (n - 1);
    }
    
    #include <iostream>
    #include "distance.h"
    #define N 64
    int main(int argc, char** argv)
    {
    	float* in = (float*)calloc(N, sizeof(float));
    	float* out = (float*)calloc(N, sizeof(float));
    	const float ref = 0.5f;
    
    	for (int i = 0; i < N; i++)
    	{
    		in[i] = scale(i, N);
    	}
    
    	distanceArray(out, in, ref, N);
    	free(in);
    	free(out);
    	return 0;
    }
    

    2.统一内存:它可以使得一个数组在设备和主机上都能访问统一内存。最关键的部分是函数cudaMallocManaged(),使用这个函数来为托管内存分配空间(在定义了指针以后,与cudaMalloc()函数使用方法类似)。托管内存数组会带来一个开发效率和执行效率的折衷。

    #include <stdio.h>
    #include <cmath>
    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    #include <helper_cuda.h>//需要手动包含此文件,不然VS识别不了
    #define TPB 32
    #define N 1024
    
    //从核函数中调用并在GPU上执行,因此需要添加__device__标识。
    __device__ float distance(float x1, float x2)//左右各两个下划线
    {
    	return sqrt((x2 - x1) * (x2 - x1));
    }
    
    //如所有的核函数一样,是从主机上调用而在设备上执行,因此需要使用修饰符__global__,并且设置返回类型为void。
    __global__ void distanceKernal(float* d_out, float* d_in, float ref)
    {
    	const int i = blockIdx.x * blockDim.x + threadIdx.x;//计算出线程索引号,然后分散到各线程里进行并行计算
    	const float x = d_in[i];
    	d_out[i] = distance(x, ref);
    	printf("i=%2d: dist from %f to %f is %f.\n", i, ref, x, d_out[i]);
    }
    
    //在CPU执行的for循环内调用,来生成输入数据。合适的修饰符是__host__,但是我们可以不添加它,因为它是默认的标识符。
    __host__ float scale(int i, int n)
    {
    	return((float)i) / (n - 1);
    }
    
    
    //统一内存:它可以使得一个数组在设备和主机上都能访问统一内存。
    //最关键的部分是函数cudaMallocManaged(),使用这个函数来为托管内存分配空间(在定义了指针以后,与cudaMalloc()函数使用方法类似)。
    //托管内存数组会带来一个开发效率和执行效率的折衷。
    int main(int argc, char** argv)
    {
    
    	//声明指向在设备和主机上访问统一内存的输入、输出数组的指针
    	float* in = NULL;
    	float* out = NULL;
    
    	//错误检查有一定代价,那就是降低了程序的性能
    	//为输入、输出数组分配统一内存
    	checkCudaErrors(cudaMallocManaged(&in, N * sizeof(float)));
    	checkCudaErrors(cudaMallocManaged(&out, N * sizeof(float)));
    	
    	const float ref = 0.5f;
    	for (int i = 0; i < N; i++)
    	{
    		in[i] = scale(i, N);
    	}
    
    	//启动核函数进行相关的计算
    	distanceKernal << <N / TPB, TPB >> > (out, in, ref);
    	checkCudaErrors(cudaGetLastError());
    	//checkCudaErrors(cudaPeekAtLastError());
    	//调用cudaDeviceSynchronize()保证在返回之前,GPU上的计算已经完成
    	checkCudaErrors(cudaDeviceSynchronize());
    
    	//释放分配给设备端数组的内存
    	checkCudaErrors(cudaFree(in));
    	checkCudaErrors(cudaFree(out));
    
    	return 0;
    }
    

    3.计算运行时间与检查错误

    #include <stdio.h>
    #include <cmath>
    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    #include <helper_cuda.h>
    #define TPB 16
    
    __host__ float scale(int i, int n)
    {
    	return (float)i / (n - 1);
    }
    __device__ float distance(float x1, float x2)//左右各两个下划线
    {
    	return sqrt((x2 - x1) * (x2 - x1));
    }
    
    __global__ void distanceKernal(float* d_out, float* d_in, float ref)
    {
    	const int i = blockIdx.x * blockDim.x + threadIdx.x;//计算出线程索引号,然后分散到各线程里进行并行计算
    	const float x = d_in[i];
    	d_out[i] = distance(x, ref);
    	printf("i=%2d: dist from %f to %f is %f.\n", i, ref, x, d_out[i]);
    }
    
    
    __host__ void distanceArray(float* out, float* in, float ref, int n)
    {
    	//创建用于计时的事件变量
    	cudaEvent_t StartMemcpy, StopMemcpy;
    	cudaEvent_t StartKernal, StopKernal;
    	checkCudaErrors(cudaEventCreate(&StartMemcpy));
    	checkCudaErrors(cudaEventCreate(&StopMemcpy));
    	checkCudaErrors(cudaEventCreate(&StartKernal));
    	checkCudaErrors(cudaEventCreate(&StopKernal));
    
    
    	float* d_in = NULL;
    	float* d_out = NULL;
    
    	checkCudaErrors(cudaMalloc(&d_in, n * sizeof(float)));
    	checkCudaErrors(cudaMalloc(&d_out, n * sizeof(float)));
    
    	//记录事件
    	checkCudaErrors(cudaEventRecord(StartMemcpy));
    
    	checkCudaErrors(cudaMemcpy(d_in, in, n * sizeof(float), cudaMemcpyHostToDevice));
    	checkCudaErrors(cudaEventRecord(StopMemcpy));
    
    	checkCudaErrors(cudaEventRecord(StartKernal));
    
    	distanceKernal << <n / TPB, TPB >> > (d_out, d_in, ref);
    	checkCudaErrors(cudaEventRecord(StopKernal));
    
    
    	checkCudaErrors(cudaMemcpy(out, d_out, n * sizeof(float), cudaMemcpyDeviceToHost));
    
    	//确保函数全部执行完毕
    	checkCudaErrors(cudaEventSynchronize(StopMemcpy));
    	checkCudaErrors(cudaEventSynchronize(StopKernal));
    
    	//计算并输出运行时间
    	float MemcpyTime = 0;
    	checkCudaErrors(cudaEventElapsedTime(&MemcpyTime, StartMemcpy, StopMemcpy));
    	float KernalTime = 0;
    	checkCudaErrors(cudaEventElapsedTime(&KernalTime, StartKernal,StopKernal));
    	printf("Kernal time is %f ms.\n", KernalTime);
    	printf("Data transfer time is %f ms.\n", MemcpyTime);
    
    	checkCudaErrors(cudaFree(d_in));
    	checkCudaErrors(cudaFree(d_out));
    }
    
    展开全文
  • Cuda学习笔记

    千次阅读 2017-04-18 19:03:07
    CUDA C简介 基本操作 读取GPU的信息 CUDA C并行编程 向量和 Julia集 线程协作 点积的计算 申请共享内存 每个线程单独工作 多个线程协同工作 保存归约结果 总的代码CUDA C简介基本操作以下是调用GPU的基本操作代码。...

    CUDA C简介

    基本操作

    以下是调用GPU的基本操作代码。代码作用是将两个数相加。
    其中要注意的是:
    1. cudaMemcpy() 函数前两个参数传递的是地址。
    2. cudaMalloc() 函数原型为:

    cudaError_t cudaMalloc (void **devPtr, size_t  size );   

    所以调用时首先将第一个参数强制转换为 (void **) 类型,再取地址 & 得到之前定义的那个一维指针地址。

    代码如下:

    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    
    #include <stdio.h>
    
    __global__ void add(int a, int b, int* c){
        *c = a + b;
    }
    
    int main()
    {
        int c[10];
        int* dev_c;
    
        cudaMalloc( (void**)&dev_c, sizeof(int) ) ;
    
        add<<<1, 1>>>(2, 7, dev_c);
        // 第一个参数和第二个参数都传的是地址.
        // int c;
        // cudaMemcpy( &c, dev_c, sizeof(int),cudaMemcpyDeviceToHost ) ;
        // 以下代码传址传的是数组c的首地址
        cudaMemcpy( c, dev_c, sizeof(int),cudaMemcpyDeviceToHost ) ;
    
        printf("%d\n", c[0]);
        cudaFree(dev_c);
        return 0;
    }

    读取GPU的信息

    用如下方式读取GPU信息:

    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    
    #include <stdio.h>
    
    int main()
    {
        cudaDeviceProp prop;
    
        int count;
        cudaGetDeviceCount(&count);
    
        for (int i = 0; i < count; ++i)
        {
            cudaGetDeviceProperties(&prop, i);
            printf("  --- General Information for device %d --- \n", i);
            printf("Name:  %s\n", prop.name);
            printf("Compute capability: %d.%d\n", prop.major, prop.minor);
            printf("Clock rate: %d\n", prop.clockRate);
            printf("Device copy overlap:");
            if(prop.deviceOverlap)
                printf("Enable\n");
            else
                printf("Disable\n");
    
            printf("Kernel execition timeout:");
            if(prop.kernelExecTimeoutEnabled)
                puts("Enable");
            else
                puts("Disable");
    
            printf("  --- Memory Information for device %d ---\n", i);
            printf("Total global mem: %ld\n", prop.totalGlobalMem);
            printf("Total constant mem: %ld\n", prop.totalConstMem);
            printf("Max mem pitch: %ld\n", prop.memPitch);
            printf("Texture Alignment: %ld\n", prop.textureAlignment);
            printf("  --- MP Information for device %d --- \n", i);
            printf("Multiprocessor count %d\n", prop.multiProcessorCount);
            printf("Shared mem per mp : %ld\n", prop.sharedMemPerBlock);
            printf("Registers per mp: %d\n", prop.regsPerBlock);
            printf("Threads in warp: %d\n", prop.warpSize);
            printf("Max threads per block: %d\n", prop.maxThreadsPerBlock);
            printf("Max therad dimensions: (%d, %d, %d)\n",  prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
    
            printf("Max grid dimensions: (%d, %d, %d)\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
            puts("");
    
        }
    
        return 0;
    }

    CUDA C并行编程

    向量和

    使用CUDA C并行编程来计算向量与向量的和:

    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    
    #include <stdio.h>
    
    const int N = 10;
    
    __global__ void add(int* a, int* b, int* c)
    {
        // blockIdx - 内置变量,保存当前执行设备代码的线程块的索引
        int tid = blockIdx.x;
        if (tid < N)
        {
            c[tid] = a[tid] + b[tid];
        }
    }
    
    int main()
    {
        int a[N], b[N], c[N];
        int *dev_a, *dev_b, *dev_c;
    
        // GPU内存分配
        cudaMalloc((void**)&dev_a, N * sizeof(int));
        cudaMalloc((void**)&dev_b, N * sizeof(int));
        cudaMalloc((void**)&dev_c, N * sizeof(int));
    
        // 在CPU上给a, b 赋值
        for (int i = 0; i < N; ++i)
        {
            a[i] = -i;
            b[i] = i * i;
        }
    
        // 将数组a, b复制到GPU
        cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice);
        cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice);
    
        // <<<a, b>>> 第一个参数a表示设备在执行核函数时使用的并行线程块的数量
        // 以下代码创建N个线程块在GPU上运行
        add<<<N, 1>>>(dev_a, dev_b, dev_c);
    
        // 将结果dev_c复制到CPU
        cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost);
    
        for (int i = 0; i < N; ++i)
        {
            printf("%d + %d = %d\n", a[i], b[i], c[i]);
        }
    
        cudaFree(dev_a);
        cudaFree(dev_b);
        cudaFree(dev_c);
    
        return 0;
    }

    Julia集

    首先要在VS和系统路径中加入glut,之后头文件导入就行了。
    CPU版本:

    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    //#include "cuda_gl_interop.h"
    #include <gl/glut.h>
    #include <stdio.h>
    
    #define DIM 1000
    
    struct CPUBitmap
    {
        unsigned char *pixels;    /*像素点的总个数*/
        int x, y;                 /*图像的长宽*/
        void *dataBlock;          /*  */
    
        void (*bitmapExit)(void*);  /*这是一个函数 */
    
        CPUBitmap( int width, int height, void *d = NULL )
        {
            pixels = new unsigned char[width * height * 4];   /*计算总的像素点数,并分配新的空间*/
            x = width;                                        /*图像的宽*/
            y = height;                                       /*图像的高*/
            dataBlock = d;                                    /* */
        }
    
        /*析构函数*/
        ~CPUBitmap()
        {
            /*删除像素点*/
            delete [] pixels;       
        }
        /*取得所有像素点*/       
        unsigned char* get_ptr( void ) const   { return pixels; }
        /*取得图片总大小*/
        long image_size( void ) const { return x * y * 4; }
    
        void display_and_exit( void(*e)(void*) = NULL )
        {
            CPUBitmap**   bitmap = get_bitmap_ptr();
            *bitmap = this;
            bitmapExit = e;
    
            // a bug in the Windows GLUT implementation prevents us from
            // passing zero arguments to glutInit()
            int c=1;
            char* dummy = "";
    
            /*glutInit,对 GLUT (OpenGl 里面的一个工具包,包含很多函数)进行初始化,这个函数必须在其它的 GLUT使用之前调用一次。其格式比较死板,一般照抄这句glutInit(&argc, argv)就可以了*/
    
            glutInit( &c, &dummy );        
    
            /*设置显示方式,其中 GLUT_RGBA 表示使用 RGBA 颜色,与之对应的还有GLUT_INDEX(表示使用索引颜色) ;GLUT_SINGLE 表示使用单缓冲,。与之对应的还有 GLUT_DOUBLE(使用双缓冲)。*/    
            glutInitDisplayMode( GLUT_SINGLE | GLUT_RGBA );
    
            /*这个也简单,设置窗口的大小*/
    
            glutInitWindowSize( x, y );
    
            /*根据前面设置的信息创建窗口。参数将被作为窗口的标题。注意:窗口被创建后,并不立即显示到屏幕上。需要调用 glutMainLoop 才能看到窗口。*/
    
            glutCreateWindow( "bitmap" );
    
            glutKeyboardFunc(Key);
    
            /* 设置一个函数,当需要进行画图时,这个函数就会被调用。*/
            glutDisplayFunc(Draw);
    
            /*显示窗口*/
    
            glutMainLoop();
    
        }
    
            // static method used for glut callbacks
        static CPUBitmap** get_bitmap_ptr( void )
        {
            static CPUBitmap   *gBitmap;
            return &gBitmap;
        }
    
            // static method used for glut callbacks
        static void Key(unsigned char key, int x, int y)
        {
            /* 如果按键按的是Esc按键,则退出程序。*/
            switch (key)
            {
                case 27:
                CPUBitmap*   bitmap = *(get_bitmap_ptr());
                if (bitmap->dataBlock != NULL && bitmap->bitmapExit != NULL)
                    bitmap->bitmapExit( bitmap->dataBlock );
            }
         }
    
        // static method used for glut callbacks
    
        /* 画图 */
        static void Draw( void )
        {
            CPUBitmap*   bitmap = *(get_bitmap_ptr());
    
            /*设置背景颜色*/
            glClearColor( 0.0, 0.0, 0.0, 1.0 );
    
            /*清除。GL_COLOR_BUFFER_BIT 表示清除颜色*/
            glClear( GL_COLOR_BUFFER_BIT );
    
            glDrawPixels( bitmap->x, bitmap->y, GL_RGBA, GL_UNSIGNED_BYTE, bitmap->pixels );
    
            /*保证前面的 OpenGL 命令立即执行(而不是让它们在缓冲区中等待)。其作用跟 fflush(stdout)类似。*/
            glFlush();
        }
    };
    
    struct cuComplex
    {
        float r;
        float i;
        cuComplex(float a, float b) : r(a), i(b) {}
        float magnitude2(void) {return r * r + i * i;}
        cuComplex operator * (const cuComplex& a){
            return cuComplex(r * a.r - i * a.i, i * a.r + r*a.i);
        }
        cuComplex operator+(const cuComplex& a)
        {
            return cuComplex(r + a.r, i + a.i);
        }
    };
    
    int julia(int x, int y)
    {
        const float scale = 1.5;
        float jx = scale * (float)((DIM >> 1) - x) / (DIM >> 1);
        float jy = scale * (float)((DIM >> 1) - y) / (DIM >> 1);
    
        cuComplex c(-0.8, 0.156);
        cuComplex a(jx, jy);
    
        int i = 0;
        for (int i = 0; i < 200; ++i)
        {
            a = a * a + c;
            if (a.magnitude2() > 1000)
            {
                return 0;
            }
        }
        return 1;
    }
    
    void kernel(unsigned char *ptr)
    {
        for (int y = 0; y < DIM; y++)
        {
            for (int x = 0; x < DIM; x++)
            {
                int offset = x + y * DIM;
    
                int juliaValue = julia(x, y);
                ptr[(offset << 2) + 0] = 255 * juliaValue;
                ptr[(offset << 2) + 1] = 0;
                ptr[(offset << 2) + 2] = 0;
                ptr[(offset << 2) + 3] = 255;
            }
        }
    }
    
    int main()
    {
        CPUBitmap bitmap(DIM, DIM);
        unsigned char *ptr = bitmap.get_ptr();
        kernel(ptr);
    
        bitmap.display_and_exit();
    
        return 0;
    }

    改写为GPU版本:

    #include "cuda_runtime.h"
    #include "device_launch_parameters.h"
    //#include "cuda_gl_interop.h"
    #include <gl/glut.h>
    #include <stdio.h>
    #include <time.h>
    
    typedef long clock_t;
    #define CLOCKS_PER_SEC ((clock_t)1000)
    #define DIM 1000
    
    struct CPUBitmap
    {
        unsigned char *pixels;    /*像素点的总个数*/
        int x, y;                 /*图像的长宽*/
        void *dataBlock;          /*  */
    
        void (*bitmapExit)(void*);  /*这是一个函数 */
    
        CPUBitmap( int width, int height, void *d = NULL )
        {
            pixels = new unsigned char[width * height * 4];   /*计算总的像素点数,并分配新的空间*/
            x = width;                                        /*图像的宽*/
            y = height;                                       /*图像的高*/
            dataBlock = d;                                    /* */
        }
    
        /*析构函数*/
        ~CPUBitmap()
        {
            /*删除像素点*/
            delete [] pixels;       
        }
        /*取得所有像素点*/       
        unsigned char* get_ptr( void ) const   { return pixels; }
        /*取得图片总大小*/
        long image_size( void ) const { return x * y * 4; }
    
        void display_and_exit( void(*e)(void*) = NULL )
        {
            CPUBitmap**   bitmap = get_bitmap_ptr();
            *bitmap = this;
            bitmapExit = e;
    
            // a bug in the Windows GLUT implementation prevents us from
            // passing zero arguments to glutInit()
            int c=1;
            char* dummy = "";
    
            /*glutInit,对 GLUT (OpenGl 里面的一个工具包,包含很多函数)进行初始化,这个函数必须在其它的 GLUT使用之前调用一次。其格式比较死板,一般照抄这句glutInit(&argc, argv)就可以了*/
    
            glutInit( &c, &dummy );        
    
            /*设置显示方式,其中 GLUT_RGBA 表示使用 RGBA 颜色,与之对应的还有GLUT_INDEX(表示使用索引颜色) ;GLUT_SINGLE 表示使用单缓冲,。与之对应的还有 GLUT_DOUBLE(使用双缓冲)。*/    
            glutInitDisplayMode( GLUT_SINGLE | GLUT_RGBA );
    
            /*这个也简单,设置窗口的大小*/
    
            glutInitWindowSize( x, y );
    
            /*根据前面设置的信息创建窗口。参数将被作为窗口的标题。注意:窗口被创建后,并不立即显示到屏幕上。需要调用 glutMainLoop 才能看到窗口。*/
    
            glutCreateWindow( "bitmap" );
    
            glutKeyboardFunc(Key);
    
            /* 设置一个函数,当需要进行画图时,这个函数就会被调用。*/
            glutDisplayFunc(Draw);
    
            /*显示窗口*/
    
            glutMainLoop();
    
        }
    
            // static method used for glut callbacks
        static CPUBitmap** get_bitmap_ptr( void )
        {
            static CPUBitmap   *gBitmap;
            return &gBitmap;
        }
    
            // static method used for glut callbacks
        static void Key(unsigned char key, int x, int y)
        {
            /* 如果按键按的是Esc按键,则退出程序。*/
            switch (key)
            {
                case 27:
                CPUBitmap*   bitmap = *(get_bitmap_ptr());
                if (bitmap->dataBlock != NULL && bitmap->bitmapExit != NULL)
                    bitmap->bitmapExit( bitmap->dataBlock );
            }
         }
    
        // static method used for glut callbacks
    
        /* 画图 */
        static void Draw( void )
        {
            CPUBitmap*   bitmap = *(get_bitmap_ptr());
    
            /*设置背景颜色*/
            glClearColor( 0.0, 0.0, 0.0, 1.0 );
    
            /*清除。GL_COLOR_BUFFER_BIT 表示清除颜色*/
            glClear( GL_COLOR_BUFFER_BIT );
    
            glDrawPixels( bitmap->x, bitmap->y, GL_RGBA, GL_UNSIGNED_BYTE, bitmap->pixels );
    
            /*保证前面的 OpenGL 命令立即执行(而不是让它们在缓冲区中等待)。其作用跟 fflush(stdout)类似。*/
            glFlush();
        }
    };
    
    struct cuComplex
    {
        float r;
        float i;
        __device__ cuComplex(float a, float b) : r(a), i(b) {}
        __device__ float magnitude2(void) {return r * r + i * i;}
        __device__ cuComplex operator * (const cuComplex& a)
        {
            return cuComplex(r * a.r - i * a.i, i * a.r + r*a.i);
        }
        __device__ cuComplex operator + (const cuComplex& a)
        { 
            return cuComplex(r + a.r, i + a.i);
        }
    };
    
    __device__ int julia(int x, int y)
    {
        const float scale = 1.5;
        float jx = scale * (float)((DIM >> 1) - x) / (DIM >> 1);
        float jy = scale * (float)((DIM >> 1) - y) / (DIM >> 1);
    
        cuComplex c(-0.8, 0.156);
        cuComplex a(jx, jy);
    
        int i = 0;
        for (int i = 0; i < 200; ++i)
        {
            a = a * a + c;
            if (a.magnitude2() > 1000)
            {
                return 0;
            }
        }
        return 1;
    }
    
    __global__ void kernel(unsigned char *ptr)
    {
        int x = blockIdx.x;
        int y = blockIdx.y;
        int offset = x + y * gridDim.x;
    
        int juliaValue = julia(x, y);
        ptr[(offset << 2) + 0] = 255 * juliaValue;
        ptr[(offset << 2) + 1] = 0;
        ptr[(offset << 2) + 2] = 0;
        ptr[(offset << 2) + 3] = 255;
    }
    
    int main()
    {
        CPUBitmap bitmap(DIM, DIM);
        unsigned char *dev_bitmap;
    
        cudaMalloc( (void**)&dev_bitmap, bitmap.image_size() );
    
        dim3 grid(DIM, DIM);
        //clock_t start = clock();
        kernel<<<grid, 1>>>(dev_bitmap);
        //clock_t finish = clock();
    
        //printf( "%f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );
    
        cudaMemcpy(bitmap.get_ptr(), dev_bitmap, bitmap.image_size(), cudaMemcpyDeviceToHost);
    
        bitmap.display_and_exit();
    
        cudaFree(dev_bitmap);
        return 0;
    }

    线程协作

    点积的计算

    申请共享内存

    首先我们需要申请共享内存,在这个例子中声明的是数组cache:

    __shared__ float cache[threadsPerBlock];

    这里我们需要明白的是,一旦这样声明数组,就会创建与线程块的数量相同的数组cahce,即每个线程块都会对应一个这样的数组cache。我们都知道,共享内存是用于同一个线程块内的线程之间交流的,不同线程块之间是无法通过共享内存进行交流的。另外,数组cache的大小是每个线程块中线程的个数,即线程块的大小。

    每个线程单独工作

    现在来看每个线程完成的是什么工作!
    如果向量长度不是特别长(假设大小等于总线程个数)的话,每个线程只需要工作一次,即计算两个元素的积并保存在中间变量 temp 里。但是实际计算过程中由于向量长度过长,一次计算可能会计算不完,每个线程需要多次计算才能完成所有工作,因此 temp 保存的值可能为多个元素乘积之和,如下图所示:
    这里写图片描述
    假设数组大小为16,线程总数为4。此时一次并行是无法完成工作的,所以需要多次并行,即每个线程需要做四次工作才可完成计算。
    相应的代码如下:

    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    double tmp = 0;
    while (tid < N)
    {
        tmp += a[tid] * b[tid];
        tid += blockDim.x * gridDim.x;
    }

    多个线程协同工作

    线程之间通过共享内存进行协作。每个线程将temp的值保存到每个线程块的共享内存(shared memory)中,即数组cache中,相应的代码如下:

    cache[cacheIndex] = temp;
    __syncthreads();

    这样每个线程块中对应的数组cache保存的就是每个线程的计算结果。为了节省带宽,这里又采用了并行计算中常用的归约算法,来计算数组中所有值之和,并保存在第一个元素(cache[0])内。这样每个线程就通过共享内存(shared memory)进行数据交流了。具体代码如下所示:

    //归约算法将每个线程块上的cache数组归约为一个值cache[0],最终保存在数组c里
    int i = blockDim.x /2;
    while (i != 0)
    {
        if (cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();        //确保每个线程已经执行完前面的语句
        i /= 2;
    }

    保存归约结果

    现在每个线程块的计算结果已经保存到每个共享数组cache的第一个元素cache[0]中,这样可以大大节省带宽。下面就需要将这些归约结果保存到全局内存(global memory)中。

    观察核函数你会发现有一个传入参数——数组c。这个数组是位于全局内存中,每次使用线程块中线程ID为0的线程来将每个线程块的归约结果保存到该数组中,注意这里每个线程块中的结果保存到数组c中与之相对应的位置,即c[blockIdx.x]。

    总的代码

    #include <stdlib.h>
    #include <stdio.h>
    #include "cuda_runtime.h"
    #include "device_functions.h"
    #include "device_launch_parameters.h"
    
    #define Min(a,b) (a<b?a:b)
    
    const int N = 33 * 1024;
    // 线程块里的线程256个,线程格一共有32个线程,这就意味着,每个线程将会计算4次,因为数组元素很大
    const int threadsPerBlock = 256;
    // 
    const int blocksPerGrid = Min( 32, (N + threadsPerBlock - 1) / threadsPerBlock );
    
    __global__ void dot(float *a, float *b, float *c)
    {
        __shared__ float cache[threadsPerBlock];
        int tid = threadIdx.x + blockIdx.x * blockDim.x;
        int cacheIndex = threadIdx.x;
    
        float tmp = 0;
        while (tid < N)
        {
            tmp += a[tid] * b[tid];
            tid += blockDim.x * gridDim.x;
        }
    
        // 设置cache中对应位置上的值
        cache[cacheIndex] = tmp;
    
        // 对线程块中的线程进行同步
        __syncthreads();
    
        // 归约,以下代码要求threadPerBlock是2的指数
        int i = blockDim.x /2;
        while (i != 0)
        {
            if (cacheIndex < i)
            {
                cache[cacheIndex] += cache[cacheIndex + i];
            }
            __syncthreads();
            i /= 2;
        }
        if (cacheIndex == 0)
        {
            c[blockIdx.x] = cache[0];
        }
    
    }
    
    int main()
    {
        float *a, *b, c, *partial_c;
        float *dev_a, *dev_b, *dev_partial_c;
    
        // 在CPU上分配内存
        a = (float*)malloc(N * sizeof(float));
        b = (float*)malloc(N * sizeof(float));
        partial_c = (float*)malloc(blocksPerGrid * sizeof(float));
    
        // 在GPU上分配内存
        cudaMalloc((void**)&dev_a, N * sizeof(float));
        cudaMalloc((void**)&dev_b, N * sizeof(float));
        cudaMalloc((void**)&dev_partial_c, blocksPerGrid * sizeof(float));
    
        // 初始化
        for (int i = 0; i < N; ++i)
        {
            a[i] = i;
            b[i] = i << 1;
        }
    
        // 将数组a, b从CPU上复制到GPU
        cudaMemcpy(dev_a, a, N * sizeof(float), cudaMemcpyHostToDevice);
        cudaMemcpy(dev_b, b, N * sizeof(float), cudaMemcpyHostToDevice);
    
        // 计算点积
        dot<<<blocksPerGrid, threadsPerBlock>>>(dev_a, dev_b, dev_partial_c);
    
        // 将partial_c从GPU上复制到CPU
        cudaMemcpy(partial_c, dev_partial_c, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost);
    
        // 在CPU上完成最终计算
        c = 0;
        for (int i = 0; i < blocksPerGrid; ++i)
        {
            c += partial_c[i];
        }
        #define sum_squares(x) (x*(x+1)*(x*2+1)/6)
        printf( "Does GPU value %.6g == %.6g?\n", c, 2 * sum_squares( (N - 1.0) ) );
        cudaFree(dev_a);
        cudaFree(dev_b);
        cudaFree(dev_partial_c);
    
        free(a);
        free(b);
        free(partial_c);
    
        return 0;
    }
    

    常量内存与事件

    常量内存

    常量内存用

    __constant__

    来声明,将把变量的访问限制为只读。与全局内存中读取数据相比,从常量内存中读取相同的数据可以节约内存带宽。
    但使用常量内存是否可以使性能变好,可以由事件来计算运行时间判断。

    事件

    CUDA中事件的本质上是一个GPU时间戳,这个时间戳是在用户制定的时间点上记录的。由于GPU本身支持记录时间戳,因此避免了当使用CPU定时器来统计GPU执行的时间时可能遇到的诸多问题。
    其使用过程如下:
    首先创建一个起始事件,结束事件,然后记录一个事件,最后告诉CPU在某个事件上需要同步。

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start, 0);
    
    // 在GPU上执行一些工作
    
    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    
    float elapsedTime;
    cudaEventElapsedTime(&elapsedTime, start, stop);
    
    printf("Time: %.3lf ms\n", elapsedTime);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    

    其中cudaEventRecord函数的第二个参数,书中会在讨论流(stream)的时候再介绍。

    原子性

    简介

    从C/C++的递增运算符入手:

    x++;

    这条语句的操作包括:

    1. 读取x中的值。
    2. 将步骤1中读到的值增加1。
    3. 将递增后的结果写回到x。

    现在考虑:如果有两个线程需要对x的值进行递增。会有非常多的调度方式,如果调度方式不正确,将会得到错误的结果。因此我们需要通过某种方式一次性执行完读取-修改-写入三个操作,并且在执行过程中不能被其他线程中断,除非已经完成了这三个操作,否则其他的线程都不能读取或写入x的值。
    由于这些操作的执行过程不能分解为更小的部分,因此我们将满足这种条件限制的操作称为原子操作
    CUDA C支持多种原子操作,当有数千个线程在内存访问上发生竞争时,这些操作能够确保在内存上实现安全的操作。

    直方图的计算

    在CPU上计算直方图

    这是一个在CPU上计算直方图的程序,非常简单。

    #include <stdlib.h>
    #include <stdio.h>
    #include <time.h>
    #include "cuda_runtime.h"
    #include "device_functions.h"
    #include "device_launch_parameters.h"
    
    #define Min(a,b) (a<b?a:b)
    #define SIZE (100*1024*1024)
    
    // 生成随机数据
    void* big_random_block( int size )
    {
        unsigned char *data = (unsigned char*)malloc( size );
        for (int i=0; i<size; i++)
            data[i] = rand();
        return data;
    }
    
    int main()
    {
        unsigned char *buffer = (unsigned char*)big_random_block( SIZE );
    
        // capture the start time
        clock_t start, stop;
        start = clock();
    
        unsigned int histo[256];
        for (int i = 0; i < 256; i++)
            histo[i] = 0;
    
        for (int i = 0; i < SIZE; i++)
            histo[buffer[i]]++;
    
        stop = clock();
        float elapsedTime = (1.0 * stop - start) / CLOCKS_PER_SEC * 1000.0f;
        printf( "Time to generate:  %3.1f ms\n", elapsedTime );
    
        long histoCount = 0;
        for (int i = 0; i < 256; i++)
        {
            histoCount += histo[i];
        }
        printf( "Histogram Sum:  %ld\n", histoCount );
    
        free( buffer );
        return 0;
    }

    在GPU上计算直方图

    如果输入数组足够大,通过多个线程处理缓冲区的不同部分,将会节约大量的计算时间。
    不同的线程来读取不同部分的输入数据非常容易。但在计算输入数据的直方图时,多个线程可能同时对输出直方图的同一个元素进行递增。在这种情况下,需要通过原子递增操作来避免问题。
    main函数与CPU版本的基本差不多。思路是先分配内存,然后调用GPU计算,然后检验结果是否正确,最后释放内存:

    int main()
    {
        unsigned char *buffer = (unsigned char*)big_random_block( SIZE );
    
        // 初始化计时事件
        cudaEvent_t start, stop;
        cudaEventCreate( &start );
        cudaEventCreate( &stop );
        cudaEventRecord( start, 0 );
    
        // 在GPU上为数据分配内存
        unsigned char *dev_buffer;
        unsigned int *dev_histo;
        cudaMalloc( (void**)&dev_buffer, SIZE );
        cudaMalloc( (void**)&dev_histo, 256 * sizeof(int) );
    
        cudaMemcpy( dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice );
        cudaMemset( dev_histo, 0, 256 * sizeof(int) );
    
        // 计算直方图
        /// .......................
        /// .......................
    
        // 将GPU上运行后的数据复制到CPU
        unsigned int histo[256];
        cudaMemcpy( histo, dev_histo, 256 * sizeof(int), cudaMemcpyDeviceToHost );
    
        // 计时结束
        cudaEventRecord( stop, 0 );
        cudaEventSynchronize( stop );
        float elapsedTime;
        cudaEventElapsedTime( &elapsedTime, start, stop );
        printf( "Time to generate: %3.lf ms\n", elapsedTime );
    
        // 在CPU上检验计算结果是否正确
        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);
            }
        }
    
        // 释放内存
        cudaEventDestroy( start );
        cudaEventDestroy( stop );
        cudaFree( dev_histo );
        cudaFree( dev_buffer );
        free( buffer );
    
        return 0;
    }

    出于性能的考虑,这个示例中的核函数调用比通常的核函数调用复杂一点。由于直方图包含256个元素,因此可以在每个线程块中包含256个线程,这种方式不仅方便而且高效。
    但是在线程块的数量上还可以有更多选择。比如在100MB数据中共有104 857 600个字节。我们可以启动一个线程块,让每个线程处理409 600个数据元素。同样,还可以启动409 600个线程块,让每个线程处理一个数据元素。
    通过一些实验,当线程块的数量为GPU处理器数量的2倍时,具有最优性能。
    通过以下代码来实现这个操作。

    cudaDeviceProp prop;
    cudaGetDeviceProperties( &prop, 0 );
    int blocks = prop.multiProcessorCount << 1;
    histo_kernel<<<blocks, 256>>>( dev_buffer, SIZE, dev_histo );

    因此完整的main函数如下:

    int main()
    {
        unsigned char *buffer = (unsigned char*)big_random_block( SIZE );
    
        // 初始化计时事件
        cudaEvent_t start, stop;
        cudaEventCreate( &start );
        cudaEventCreate( &stop );
        cudaEventRecord( start, 0 );
    
        // 在GPU上为数据分配内存
        unsigned char *dev_buffer;
        unsigned int *dev_histo;
        cudaMalloc( (void**)&dev_buffer, SIZE );
        cudaMalloc( (void**)&dev_histo, 256 * sizeof(int) );
    
        cudaMemcpy( dev_buffer, buffer, SIZE, cudaMemcpyHostToDevice );
        cudaMemset( dev_histo, 0, 256 * sizeof(int) );
    
        // GPU计算直方图
        cudaDeviceProp prop;
        cudaGetDeviceProperties( &prop, 0 );
        int blocks = prop.multiProcessorCount << 1;
        histo_kernel<<<blocks, 256>>>( dev_buffer, SIZE, dev_histo );
    
    
        // 将GPU上运行后的数据复制到CPU
        unsigned int histo[256];
        cudaMemcpy( histo, dev_histo, 256 * sizeof(int), cudaMemcpyDeviceToHost );
    
        // 计时结束
        cudaEventRecord( stop, 0 );
        cudaEventSynchronize( stop );
        float elapsedTime;
        cudaEventElapsedTime( &elapsedTime, start, stop );
        printf( "Time to generate: %3.lf ms\n", elapsedTime );
    
        // 在CPU上检验计算结果是否正确
        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);
            }
        }
    
        // 释放内存
        cudaEventDestroy( start );
        cudaEventDestroy( stop );
        cudaFree( dev_histo );
        cudaFree( dev_buffer );
        free( buffer );
    
        return 0;
    }

    全局内存原子操作

    在GPU上计算直方图的代码,首先先采用全局内存的原子操作。
    代码如下:

    __global__ void histo_kernel( unsigned char *buffer, long size, unsigned int *histo )
    {
        int i = threadIdx.x + blockIdx.x * blockDim.x;
        int stride = blockDim.x * gridDim.x;
    
        // 每个线程知道它的起始偏移i以及递增数量
        // 遍历输入数组,递增直方图中的元素
        while ( i < size )
        {
                // 原子操作
            atomicAdd( &(histo[buffer[i]]), 1 );
            i += stride;
        }
    
    }

    然而,发现在运行这个代码比原来CPU版本运行的还慢大概4倍。
    设置基性能标准很重要。
    由于在核函数里面只包含了非常少的计算工作,所以可能是全局内存上的原子操作导致性能降低。当数千个线程尝试访问少量的内存位置时,将会发生大量的竞争。为了确保递增操作的原子性,对相同内存位置的操作都将被硬件串行化。这可能导致保存未完成操作的队列非常长,会抵消通过并行运行的线程获得的性能提升。
    因此考虑共享内存的操作。

    共享内存原子操作

    虽然原子操作是导致上面性能降低的原因,但是解决这个问题的方法确实使用更多的原子操作。因为问题出在有数千个线程在少量的内存地址上发生竞争。解决这个问题分两步。
    首先,对每个并行线程块计算它所处理数据的直方图。由于每个线程块在执行这个操作时是相互独立的,所以可以在共享内存中计算这些直方图。但这种方式依然需要原子操作,因为在线程块中的多个线程之间还是会处理相同值的数据元素。但现在只有256个线程在256个地址上发生竞争,将大大减少在使用全局内存时数千个线程之间发生竞争的情况。
    然后,在上个阶段中分配一个共享内存缓冲区进行初始化,用来保存每个线程块的临时直方图。
    代码如下:

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

    因为使用256个线程,且直方图中刚好包含256个元素,因此每个线程将自动把它计算得到的元素只增加到最终直方图的元素上。如果线程数量不等于元素数量,那么这个阶段将更为复杂。

    总结

    在一些情况中,比如成千上万个线程同时修改一个地址的元素,在这些情况中,大规模并行机器反而会带来负担。而硬件中支持的原子操作可以帮助减轻这种痛苦。
    然而,在上面例子中可以看到,有时候依赖原子操作会带来性能问题,并且这些问题只能通过对算法的某些部分重构来解决。在上面例子中用了一种两阶段的算法,降低了在全集内存访问上竞争的程度。
    通常,这种降低内存竞争程度的策略总能带来不错的效果。因此如果遇到在程序中使用原子操作的时候,要记住这种策略

    散列表的实现

    散列表

    散列表实际上就是一个hash table,但我目前没搞懂hash table和inverted list的关系,感觉是同一个东西。
    散列表是一种保存 键-值 二元组的数据结构。散列表根据与值相应的键,把值放入 桶(bucket) 中。这种将键值映射到值的方法叫做散列函数。好的散列函数可以把键均匀地分布到所有桶中。这种情况下,有可能发生哈希冲突,哈希冲突的解决方法是表头键后面跟上一个链表,来保存被散列函数映射到同一个键值的桶中。

    CPU实现

    散列表主要包含两个部分:一个散列函数,一个表示桶的数据结构。
    桶的实现用分配一个长度为N的数组来表示,数组中每个元素都表示一个 键-值 的二元组链表。
    以下为这个数据结构:

    struct Entry
    {
        // 键
        unsigned int key;
        // 值,任意数据类型
        void* value;
        // 冲突时指向下一结点的指针
        Entry* next;
    }
    
    struct Table
    {
        // 哈希表的长度
        size_t cnt;
        // 哈希表,第一维指针连接表头,第二维指针指向entry
        Entry **entries;
        // 每添加一个Entry结点时需要重新分配新的内存,对程序性能产生影响
        // 用pool来维持一个可用Entry节点的数组,避免这种情况
        Entry *pool;
        // 指向下一个可用的Entry节点
        // 需要将一个结点添加到哈希表中时,只需使用firstFree指向的Entry
        Entry *firstFree;
    }

    初始化代码:

    void initialize(Table &table, int entries, int elements)
    {
        table.cnt = entries;
        // calloc分配空间并初始化为0,malloc分配空间不初始化
        table.entries = (Entry**)calloc( entries, sizeof(Entry*) );
        table.pool = (Entry*)malloc( elements * sizeof(Entry) );
        table.firstFree = table.pool;
    }

    在初始化过程中,主要操作有为哈希表entries分配内存,为结点的池分配内存,将指针firstFree初始化为指向结点池中的第一个结点。

    程序结束之后,需要释放内存:

    void free_table(Table &table)
    {
        free( table.entries );
        free( table.pool );
    }

    直接使用键值作为索引,也就是说,将结点e保存在table.entries[e.key]中。因此散列表函数如下,下面这个函数并不能保证生成数据的均匀的,这里假设生成的键是随机并且均匀的。

    size_t hash(unsigned int key, size_t cnt)
    {
        return key % cnt;
    }

    接下来是插入操作:
    1. 首先将键放入散列表函数中计算出新的结点所属于的桶。
    2. 从结点池中取出一个预先分配的Entry结点,赋值。
    3. 将这个结点插入到得到的桶的首部。
    代码如下:

    void add_to_table(Table &table, unsigned int key, void* value)
    {
        // 计算要插入的新结点的表头
        size_t hashValue = hash(key, table.cnt);
    
        // 从结点池中取出一个预先分配Entry结点
        Entry* location = table.firstFree++;
        location -> key = key;
        location -> value = value;
    
        // 插入当前表的链表首部
        location -> next = table.entries[hashValue];
        table.entries[hashValue] = location;
    }

    用如下代码来检验上面代码能否工作。首先遍历这张哈希表,然后查看每个结点,将结点放入散列表函数计算,确认这个结点被保存到了正确的桶中,检查完每个结点之后,验证散列表中的结点数量确实等于添加到散列表的元素数量。如果这些数值不相等,不是无意中将一个结点添加到了多个桶,就是没有正确的插入结点。

    void verigy_table(const Table &table)
    {
        int cnt = 0;
        for (size_t i = 0; i < table.cnt; ++i)
        {
            Entry *current = table.entries[i];
            while (current != NULL)
            {
                cnt++;
                if (hash( current->key, table.cnt) != i)
                {
                    printf("%d hashed to %ld, but was located at %ld\n", current->key, hash(current->key, table.cnt), i);
                }
                current = current -> next;
            }
        }
        if (cnt != ELEMENTS)
            printf("%d elements found in hash table. Should be %ld\n", cnt, ELEMENTS);
        else
            printf("All %d elements found in hash table.\n", cnt);
    }
    
    #define HASH_ENTRIES 1024
    
    int main()
    {
        unsigned int *buffer = (unsigned int*)big_random_block( SIZE );
        clock_t start, stop;
        start = clock();
    
        Table table;
        initialize( table, HASH_ENTRIES, ELEMENTS );
    
        for (int i = 0; i < ELEMENTS; i++)
        {
            add_to_table( table, buffer[i], (void*)NULL);
        }
    
        stop = clock();
        double elaspsedTime = (stop - start) / CLOCKS_PER_SEC * 1000.0;
        printf("Time to hash: %3.lf ms\n", elaspsedTime);
    
        verigy_table(table);
        free_table(table);
        free(buffer);
    
        return 0;
    }

    GPU多线程下的散列表

    当两个线程,同时对同一个表头插入结点的时候,就会出现两个指针同时指向原表头的情况。因此,每次只有一个线程可以安全地对表头进行插入结点。如果每个表头都有一个相应的原子锁,那么我们可以确保每次只有一个线程对指定的桶进行修改。

    首先,我们先需要一个原子锁结构,其定义如下:

    struct Lock
    {
        int *mutex;
        Lock( void )
        {
            cudaMalloc( (void**)&mutex, sizeof(int) );
            cudaMemset( mutex, 0, sizeof(int) );
        }
    
        ~Lock( void )
        {
            cudaFree( mutex );
        }
    
        __device__ void lock( void )
        {
            while( atomicCAS( mutex, 0, 1 ) != 0 );
        }
    
        __device__ void unlock( void )
        {
            atomicExch( mutex, 0 );
        }
    };

    其他数据结构的定义相同,只需要把散列表函数的声明改为 _ _ device _ _ , _ _ host _ _,当这两个关键字一起使用时,会告诉NVIDIA编译器,同时生成函数在设备上和主机上的版本。设备版本将在设备上运行,并且只能从设备代码中调用。同样,主机版本的函数将在主机上运行,并且只能从主机代码中调用。

    struct Entry
    {
        // 键
        unsigned int key;
        // 值,任意数据类型
        void* value;
        // 冲突时指向下一结点的指针
        Entry* next;
    };
    
    struct Table
    {
        // 哈希表的长度
        size_t cnt;
        // 哈希表,第一维指针连接表头,第二维指针指向entry
        Entry **entries;
        // 每添加一个Entry结点时需要重新分配新的内存,对程序性能产生影响
        // 用pool来维持一个可用Entry节点的数组,避免这种情况
        Entry *pool;
        // 指向下一个可用的Entry节点
        // 需要将一个结点添加到哈希表中时,只需使用firstFree指向的Entry
        Entry *firstFree;
    };
    
    __device__ __host__ size_t hash(unsigned int key, size_t cnt)
    {
        return key % cnt;
    }

    初始化和释放内存的函数大多与CPU版本中相同,但使用的是CUDA开辟内存的函数:

    void initialize(Table &table, int cnt, int elements)
    {
        table.cnt = cnt;
        cudaMalloc( (void**)&table.entries, cnt * sizeof(Entry*) );
        cudaMemset( table.entries, 0, cnt * sizeof(Entry*) );
        cudaMalloc( (void**)&table.pool, elements * sizeof(Entry) );
    }
    
    void free_table(Table &table)
    {
        cudaFree( table.entries );
        cudaFree( table.pool );
    }

    检查散列表的函数,可以编写一个在GPU上运行的函数,也可以使用原来CPU上的检查函数。第二种方法比较好,可以函数复用,节约开发时间。
    这里的verify_table()函数与CPU中的完全相同。由于选择了重用CPU版本的函数,因此需要把散列表从GPU内存复制到主机内存。这个函数将包括三个步骤。首先为散列表数据分配主机内存,通过cudaMemcpy()函数将GPU上的数据复制到这块内存里,这部分代码并不困难。

    复杂的地方在于,有一部分的数据是指针。不能简单地将这些指针复制到主机上,因为这些指针指向的地址存在与GPU上,他们在主机上并不是有效的指针。但这些指针的相对偏移还是有效的,每个指向Entry结点的GPU指针都指向数据table.pool[]中的某个位置,但是为了在主机上使用散列表,需要他们指向数组hostTable.pool[]中相同的Entry。因此给定一个GPU指针X,需要给目前的CPU指针加上偏移:

    hostTable.pool + (X - table.pool)

    对每个被复制的Entry指针,都要执行这个更新操作:包括hostTable.entries中的Entry指针,以及散列表的结点池中每个Entry的next指针。

    void copy_table_to_host(const Table &table, Table &hostTable)
    {
        // 创建CPU的空间,并将GPU上的数据复制到CPU上
        hostTable.cnt = table.cnt;
        hostTable.entries = (Entry**)calloc( table.cnt, sizeof(Entry*) );
        hostTable.pool = (Entry*)malloc( ELEMENTS * sizeof(Entry) );
    
        cudaMemcpy( hostTable.entries, table.entries, table.cnt * sizeof(Entry*), cudaMemcpyDeviceToHost );
        cudaMemcpy( hostTable.pool, table.pool, ELEMENTS * sizeof(Entry), cudaMemcpyDeviceToHost );
    
        // 原来复制到CPU上的指针,所指向的地址仍旧是GPU上的地址,但其偏移是不变的
        // 因此,计算复制到CPU上的GPU指针的偏移,用来重新定位在GPU当前这个内存上的元素
    
        // 重新定位在GPU上元素的表头
        for (int i = 0; i < table.cnt; ++i)
        {
            if (hostTable.entries[i] != NULL)
            {
                // hostTable.entries[i]在GPU中指向table.pool[X],减去GPU中的位置首元素table.pool,得到偏移
                // 用GPU地址的hostTable.pool + 偏移,就得到GPU上的指针指向的元素了
                hostTable.entries[i] =(Entry*)( (size_t)hostTable.pool + ((size_t)hostTable.entries[i] - (size_t)table.pool) );
            }
        }
        // 重新定位每个元素的next指针
        for (int i = 0; i < ELEMENTS; ++i)
        {
            if (hostTable.pool[i].next != NULL)
            {
                // 与上面类似
                hostTable.pool[i].next = (Entry*)( (size_t)hostTable.pool + ((size_t)hostTable.pool[i].next - (size_t)table.pool) );
            }
        }
    }

    接下来就是CUDA C原子锁语句的使用了。核函数add_to_table()的参数包括一个键的数组,一个值的数组,一个散列表和原子锁数组。原子锁数组用于锁定散列表中的每个桶。
    由于输入的数据是两个数组,并且在线程中需要对这两个数组进行索引,因此还需要将索引线性化。

    之后遍历输入数组。对于数据key[]中的每个键,线程将通过散列表函数计算出这个 键-值 二元组属于哪个桶。计算出目标桶之后,线程会锁定这个桶,添加它的 键-值 二元组,然后解锁这个桶:

    __global__ void add_to_table(unsigned int *keys, void **values, Table table, Lock *lock)
    {
        // 计算当前所在线程索引
        int tid = threadIdx.x + blockIdx.x * blockDim.x;
        int stride = blockDim.x * gridDim.x;
    
        //遍历输入数组
        while (tid < ELEMENTS)
        {
            unsigned int key = keys[tid];
            size_t hashValue = hash( key, table.cnt );
            for (int i = 0; i < 32; i++)
            {
                if ( (tid % 32) == i )
                {
                    Entry* location = &(table.pool[tid]);
                    location->key = key;
                    location->value = values[tid];
    
                    // 原子锁,锁定哈希表头的内存,这块内存只能当前线程操作
                    lock[hashValue].lock();
    
                    location->next = table.entries[hashValue];
                    table.entries[hashValue] = location;
    
                    // 解锁原子锁
                    lock[hashValue].unlock();
                }
            }
            tid += stride;
        }
    
    }

    然而这段代码中有一个非常特别的地方,for()循环和后面的if()语句似乎没必要。在之前有引入线程束的概念,线程束是一个包含32个线程的集合,并且这些线程以不掉已知的方式执行。这本书中并没有讨论如何在GPU上实现这种步调一致的执行方式,但每次在线程束中只有一个线程可以获得这个锁,如果让线程束中所有的32个线程都同时竞争这个锁,将会发生严重的问题。在这种情况下,最好的方式是在软件中执行一部分的工作,遍历线程束中的线程,并给每个线程依此机会来获取数据结构的锁,执行工作,然后解锁。

    main函数与CPU版本大致相同。首先分配一大块随机数据作为散列表的键。然后用CUDA事件来计时。接下来为随机数组分配GPU内存,将数组幅值到GPU上,并初始化散列表。
    之后的步骤就是为散列表的桶准备好原子锁,为散列表中每一个桶都分配一个锁。
    总的代码如下:

    #include <stdlib.h>
    #include <stdio.h>
    #include <time.h>
    #include "cuda_runtime.h"
    #include "device_functions.h"
    #include "device_launch_parameters.h"
    
    #define Min(a,b) (a<b?a:b)
    #define SIZE (100*1024*1024)
    #define ELEMENTS (SIZE / sizeof(unsigned int))
    
    struct Lock
    {
        int *mutex;
        Lock( void )
        {
            cudaMalloc( (void**)&mutex, sizeof(int) );
            cudaMemset( mutex, 0, sizeof(int) );
        }
    
        ~Lock( void )
        {
            cudaFree( mutex );
        }
    
        __device__ void lock( void )
        {
            while( atomicCAS( mutex, 0, 1 ) != 0 );
        }
    
        __device__ void unlock( void )
        {
            atomicExch( mutex, 0 );
        }
    };
    
    // 生成随机数据
    void* big_random_block( int size )
    {
        unsigned char *data = (unsigned char*)malloc( size );
        for (int i=0; i<size; i++)
            data[i] = rand();
        return data;
    }
    
    struct Entry
    {
        // 键
        unsigned int key;
        // 值,任意数据类型
        void* value;
        // 冲突时指向下一结点的指针
        Entry* next;
    };
    
    struct Table
    {
        // 哈希表的长度
        size_t cnt;
        // 哈希表,第一维指针连接表头,第二维指针指向entry
        Entry **entries;
        // 每添加一个Entry结点时需要重新分配新的内存,对程序性能产生影响
        // 用pool来维持一个可用Entry节点的数组,避免这种情况
        Entry *pool;
        // 指向下一个可用的Entry节点
        // 需要将一个结点添加到哈希表中时,只需使用firstFree指向的Entry
        Entry *firstFree;
    };
    
    void initialize(Table &table, int cnt, int elements)
    {
        table.cnt = cnt;
        cudaMalloc( (void**)&table.entries, cnt * sizeof(Entry*) );
        cudaMemset( table.entries, 0, cnt * sizeof(Entry*) );
        cudaMalloc( (void**)&table.pool, elements * sizeof(Entry) );
    }
    
    void free_table(Table &table)
    {
        cudaFree( table.entries );
        cudaFree( table.pool );
    }
    
    __device__ __host__ size_t hash(unsigned int key, size_t cnt)
    {
        return key % cnt;
    }
    
    void copy_table_to_host(const Table &table, Table &hostTable)
    {
        // 创建CPU的空间,并将GPU上的数据复制到CPU上
        hostTable.cnt = table.cnt;
        hostTable.entries = (Entry**)calloc( table.cnt, sizeof(Entry*) );
        hostTable.pool = (Entry*)malloc( ELEMENTS * sizeof(Entry) );
    
        cudaMemcpy( hostTable.entries, table.entries, table.cnt * sizeof(Entry*), cudaMemcpyDeviceToHost );
        cudaMemcpy( hostTable.pool, table.pool, ELEMENTS * sizeof(Entry), cudaMemcpyDeviceToHost );
    
        // 原来复制到CPU上的指针,所指向的地址仍旧是GPU上的地址,但其偏移是不变的
        // 因此,计算复制到CPU上的GPU指针的偏移,用来重新定位在GPU当前这个内存上的元素
    
        // 重新定位在GPU上元素的表头
        for (int i = 0; i < table.cnt; ++i)
        {
            if (hostTable.entries[i] != NULL)
            {
                // hostTable.entries[i]在GPU中指向table.pool[X],减去GPU中的位置首元素table.pool,得到偏移
                // 用GPU地址的hostTable.pool + 偏移,就得到GPU上的指针指向的元素了
                hostTable.entries[i] =(Entry*)( (size_t)hostTable.pool + ((size_t)hostTable.entries[i] - (size_t)table.pool) );
            }
        }
        // 重新定位每个元素的next指针
        for (int i = 0; i < ELEMENTS; ++i)
        {
            if (hostTable.pool[i].next != NULL)
            {
                // 与上面类似
                hostTable.pool[i].next = (Entry*)( (size_t)hostTable.pool + ((size_t)hostTable.pool[i].next - (size_t)table.pool) );
            }
        }
    }
    
    void verify_table(const Table &dev_table)
    {
        int cnt = 0;
        Table table;
        copy_table_to_host(dev_table, table);
        for (size_t i = 0; i < table.cnt; ++i)
        {
            Entry *current = table.entries[i];
            while (current != NULL)
            {
                cnt++;
                if (hash( current->key, table.cnt) != i)
                {
                    printf("%d hashed to %ld, but was located at %ld\n", current->key, hash(current->key, table.cnt), i);
                }
                current = current -> next;
            }
        }
        if (cnt != ELEMENTS)
            printf("%d elements found in hash table. Should be %ld\n", cnt, ELEMENTS);
        else
            printf("All %d elements found in hash table.\n", cnt);
    }
    
    /// CPU version
    //void add_to_table(Table &table, unsigned int key, void* value)
    //{
    //  // 计算要插入的新结点的表头
    //  size_t hashValue = hash(key, table.cnt);
    //
    //  // 从结点池中取出一个预先分配Entry结点
    //  Entry* location = table.firstFree++;
    //  location -> key = key;
    //  location -> value = value;
    //
    //  // 插入当前表的链表首部
    //  location -> next = table.entries[hashValue];
    //  table.entries[hashValue] = location;
    //}
    
    /// GPU version
    __global__ void add_to_table(unsigned int *keys, void **values, Table table, Lock *lock)
    {
        // 计算当前所在线程索引
        int tid = threadIdx.x + blockIdx.x * blockDim.x;
        int stride = blockDim.x * gridDim.x;
    
        //遍历输入数组
        while (tid < ELEMENTS)
        {
            unsigned int key = keys[tid];
            size_t hashValue = hash( key, table.cnt );
            for (int i = 0; i < 32; i++)
            {
                if ( (tid % 32) == i )
                {
                    Entry* location = &(table.pool[tid]);
                    location->key = key;
                    location->value = values[tid];
    
                    // 原子锁,锁定哈希表头的内存,这块内存只能当前线程操作
                    lock[hashValue].lock();
    
                    location->next = table.entries[hashValue];
                    table.entries[hashValue] = location;
    
                    // 解锁原子锁
                    lock[hashValue].unlock();
                }
            }
            tid += stride;
        }
    
    }
    
    #define HASH_ENTRIES 1024
    
    int main()
    {
        unsigned int *buffer = (unsigned int*)big_random_block( SIZE );
    
        cudaEvent_t start, stop;
        cudaEventCreate( &start );
        cudaEventCreate( &stop );
        cudaEventRecord( start, 0 );
    
        unsigned int *dev_keys;
        void **dev_values;
    
        cudaMalloc( (void**)&dev_keys, SIZE );
        cudaMalloc( (void**)&dev_values, SIZE );
        cudaMemcpy( dev_keys, buffer, SIZE, cudaMemcpyHostToDevice );
    
        // 分配锁
        Table table;
        initialize( table, HASH_ENTRIES, ELEMENTS );
    
        Lock lock[HASH_ENTRIES];
        Lock* dev_lock;
    
        cudaMalloc( (void**)&dev_lock, HASH_ENTRIES * sizeof(Lock) );
        cudaMemcpy( dev_lock, lock, HASH_ENTRIES * sizeof(Lock), cudaMemcpyHostToDevice );
    
        // kernel
        add_to_table<<<60, 256>>>(dev_keys, dev_values, table, dev_lock);
    
        cudaEventRecord( stop, 0 );
        cudaEventSynchronize( stop );
        float elapsedTime;
        cudaEventElapsedTime( &elapsedTime, start, stop );
        printf("Time to hash: %3.lf ms\n", elapsedTime);
    
        verify_table(table);
    
        cudaEventDestroy( start );
        cudaEventDestroy( stop );
        free_table( table );
        cudaFree( dev_lock );
        cudaFree( dev_keys );
        cudaFree( dev_values );
        free( buffer );
    
    
        return 0;
    }
    
    展开全文
  • cuda学习笔记[part 1]

    2021-02-03 18:50:18
    cuda学习笔记[part 1] cuda学习笔记[part 1]cuda学习笔记[part 1]环境程序配置编写测试程序-----打印GPU参数。参考文章 环境 我这里的环境默认为CUDA10.1,vs2019 1.右键项目属性->包含目录->编辑,加入你的...
  • Pytorch-CUDA学习笔记

    2020-12-09 16:29:51
    Pytorch-CUDA学习笔记 C++ FRONTEND C++的高层封装,是Pytorch的C++版本。PyTorch利用CPython在front end的基础上添加一个胶水层,使使用者能够用Python调用。 #include <torch/torch.h> //引入包 auto model...
  • CUDA学习笔记(LESSON1/2)——架构、通信模式与GPU硬件 CUDA学习笔记(LESSON3)——GPU基本算法(Part I) CUDA学习笔记(LESSON4)——GPU基本算法(Part II) CUDA学习笔记(LESSON5)——GPU优化 CUDA学习...
  • CUDA学习笔记之 CUDA存储器模型 标签: cuda存储bindingcache编程api 2010-12-14 01:33 1223人阅读 评论(0) 收藏 举报 分类: CUDA(26) GPU片内:register,shared memory; 板载...
  • Cuda 学习笔记 (一)

    2019-03-06 21:03:10
    Cuda 学习笔记 (一) 1.Hello World #include&lt;stdio.h&gt; __global__ void helloFromGPU(void){ printf("Hello World from GPU!\n"); } int main(void){ printf("Hello World from ...
  • CUDA学习笔记(LESSON1/2)——架构、通信模式与GPU硬件 CUDA学习笔记(LESSON3)——GPU基本算法(Part I) CUDA学习笔记(LESSON4)——GPU基本算法(Part II) CUDA学习笔记(LESSON5)——GPU优化 CUDA学习...
  • CUDA学习笔记(LESSON1/2)——架构、通信模式与GPU硬件 CUDA学习笔记(LESSON3)——GPU基本算法(Part I) CUDA学习笔记(LESSON4)——GPU基本算法(Part II) CUDA学习笔记(LESSON5)——GPU优化 CUDA学习...
  • CUDA学习笔记(LESSON1/2)——架构、通信模式与GPU硬件 CUDA学习笔记(LESSON3)——GPU基本算法(Part I) CUDA学习笔记(LESSON4)——GPU基本算法(Part II) CUDA学习笔记(LESSON5)——GPU优化 CUDA学习...
  • CUDA学习笔记(1) 必备工具 NVIDIA的显卡 CUDA编译器NVCC windwos可去官网下载(https://developer.nvidia.com/cuda-downloads) linux(本机为manjiaro) sudo pacman -Ss //查找 sudo pacman -S //安装 CUDA,...
  • CUDA学习笔记之程序优化 标签: cuda优化conflict存储算法数学计算 2010-01-05 17:18 5035人阅读 评论(4) 收藏 举报 分类: CUDA(6) 版权声明:本文为博主原创文章,未经博主允许...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,433
精华内容 573
关键字:

cuda学习笔记