风辰的 CUDA 入门教程


1 风辰的 CUDA 入门教程 作者:风辰 二零一零年七月二十四日 于中国科学院研究生院青年公寓 基于共同进步、分享的原则,任何个人都可使用此文档,但是本人保留所有权利。 2 目录 风辰的 CUDA 入门教程 ........................................................................................................................ 1 第一章、CUDA 的基本内容 ........................................................................................................ 3 第一节、CUDA 及 GPU 简介 ............................................................................................. 3 第二节、Linux 下 CUDA 开发环境安装 .......................................................................... 3 第三节、CUDA 与 fork/join 模式 .................................................................................... 4 第四节、CUDA C 语言 ........................................................................................................ 5 第五节、计算π.................................................................................................................. 6 第六节、编程模式.............................................................................................................. 8 第七节、线程层次.............................................................................................................. 9 第八节、存储器组织 ....................................................................................................... 10 第九节、执行模式............................................................................................................ 16 第十节、NVIDIA GPU 结构 .............................................................................................. 17 第二章、CUDA 程序优化 ......................................................................................................... 20 每一节、CUDA 总体优化策略 .......................................................................................... 20 第二节、计时器的设计 ................................................................................................... 20 第三节、错误处理............................................................................................................ 22 第四节、串行 C 程序的优化 ........................................................................................... 25 第五节、CUDA程序的优化 ....................................................................................... 27 第三章、一些例子.................................................................................................................... 30 第一节、 两向量的距离 .................................................................................................. 30 第二节、 矩阵与向量乘积 .............................................................................................. 32 第三节、 线性方程组的求解 .......................................................................................... 35 3 第一章、CUDA 的基本内容 第一节、CUDA 及 GPU 简介 GPU 是图形处理单元 (Graphic Processing Unit) 的简称,最初主要用于图形渲染。自 九十年代开始,GPU 的发展产生了较大的变化,NVIDIA、AMD(ATI)等 GPU 生产商敏锐的观察 到 GPU 天生的并行性,经过他们对硬件和软件的改进,GPU 的可编程能力不断提高,GPU 通 用计算应运而生。由于 GPU 具有比 CPU 强大的计算能力(见图 1-1),为科学计算的应用提 供了新的选择。 图 1-1 GPU 与 CPU 峰值计算能力对比图 由于 GPU 拥有比 CPU 更强的计算能力,很早就有人想到将 GPU 应用到通用计算上,这就 是 GPGPU,所谓 GPGPU 是指直接使用了图形学的 API,将任务映射成纹理的渲染过程,使用 汇编或者高级着色器语言 Cg,HLSL 等编写程序,然后通过图形学 API 执行(Direct3D 和 OpenGL),这样的开发不仅难度较大,而且难以优化,对开发人员的要求非常高,因此,传 统的 GPGPU 计算并没有广泛应用。 2007 年 6 月,NVIDIA 公司推出了 CUDA (Compute Unified Device Architecture) , CUDA 不需要借助图形学 API,而是采用了类 C 语言进行开发。同时,CUDA 采用了统一处理 架构,降低了编程的难度,使得 NVIDIA 相比 AMD/ATI 后来居上。相比 AMD 的 GPU,NVIDIA GPU 引入了片内共享存储器,提高了效率。这两项改进使 CUDA 架构更加适合进行 GPU 通用 计算。由于这些特性,CUDA 推出后迅速发展,被应用于石油勘测、天文计算、流体力学模 拟、分子动力学仿真、生物计算、图像处理、音视频编解码等领域。 CUDA 直接采用 C/C++编译器作为前端(如 Linux 下的 gcc,windows 下的 vs),以 C/C++ 语法为基础设计,因此对熟悉 C 系列语言的程序员来说,CUDA 的语法比较容易掌握。但是 这并不意味着 CUDA 容易,整体来说 CUDA 是相当难的,难在优化,难在开发出健壮、可扩展 性的程序;难在没有成熟的库和算法以借用。从语言的角度说,CUDA 只对 ANSI C 进行了最 小的必要扩展,以实现其关键特性--线程按照两个层次进行组织、共享存储器和栅栏同步。 这些关键特性使得 CUDA 拥有了两个层次的并行:线程级并行实现的细粒度数据并行,和任 务级并行实现的粗粒度并行。 第二节、Linux 下 CUDA 开发环境安装 前一节已经简单了说了一下 CUDA,为了能够使用 CUDA 开发,这一节将说明怎样构建 CUDA 开发环境。 4 本节介绍在 ubuntu9.04 操作系统和 gcc 前端的基础上安装 CUDA 开发环境。 首先,要保证自己机器上的 gcc 能够使用,因为 Ubuntu 缺少 gcc 的一些包和 g++,所以 这些得自己安装。安装命令:sudo apt-get install g++,待其完成后,弄个 C 代码试试 看:);当然你得保证你的显卡支持 CUDA。 其次,到 nVidia 官方网站(http://www.nvidia.cn/object/cuda_get_cn.html)上下载 对应操作系统的驱动(driver)和工具包(toolkit)。 再次,转换到控制台,命令为 Ctrl+Alt+F1/F2/F3/F4,关掉 gdm,命令为:sudo /etc/init.d/gdm stop,要确定已经关闭,否则在安装时会提示你有x server 程序在运 行。 再次,进入 driver 和 toolkit 目录,执行安装命令,为了方便,请一定按照默认安装。 命令 sudo sh ./dev……………… sudo sh ./nv………… 然后,打开个人目录下的.bashrc 文件或者/etc/profile 文件,在其中加入命令: PATH=${PATH}:/usr/local/cuda/bin/ LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/usr/local/cuda/lib,执行 source .bashrc 或者 /etc/profile,依据你添加 PATH 和 LD_LIBRARY_PATH 时修改了那个文件确定。 最后执行 nvcc 命令,看看,如果提示你没有输入文件,就安装完成了。 如果你要安装 SDK 的话,还得安装一些包,依据 make 时的提示,Google 和新力得应该 可以搞定一切,现在你可以享受 CUDA 了! 第三节、CUDA 与 fork/join 模式 说到 CUDA 就不能不说一下异构并行计算,所谓异构就是指多个计算平台,如 CPU,GPU 和其它一些加速部件。一般而言,各个计算平台有其独立的硬件资源。CPU 和 GPU 拥有各 自独立的 ALU、存储器。本文中的内存特指 CPU 内存,而涉及到 GPU 的存储器,则用全称。 CUDA 认为可以用于计算的硬件系统包含两个部分,一个是 CPU,一个是 GPU,CPU 控制/ 指挥 GPU 工作,GPU 只是 CPU 的协处理器。 熟悉 Linux、java、pthread 或 OpenMP 的人对于 fork/join 一定不会陌生,下图给出了 fork/join 模式的示意图。在 Linux 中,fork 函数会产生子进程,子进程和父进程共同工 作,父进程调用 wait/waitpid 函数来等待子进程;在 pthread 中,pthread_create 函数产 生子线程,子线程执行指定的工作,在某个位置,父线程调用 pthread_join 来等待子线程 完成工作;java Thread 中有一个函数 join 也是等待某个运行至此,在 java 7 中,则直接 引入 fork/join 框架;当然,这些语文/库相应的函数功能都远超过这里所说的。在这一点 上和 CUDA 最相似的是 OpenMP,其#pragma omp parallel 指令表示下面的一个代码块是由 多个线程执行的,到块的结束,所有的线程又都自动消失,只留下主线程。CUDA 也是这 样,CUDA 中多线程执行的块称为内核,内核有其特殊的声明和调用函数,这在下一节说 明。 5 图 1-3.fork/join 模式 第四节、CUDA C 语言 对一个语言来说,最简单的,也是用得最多的当然是它的语法了,下面简单的介绍一下 CUDA 的语法。 CUDA C 不是 C 语言,而是对 C 语言进行扩展。CUDA 对 C 的扩展主要包括以下四个方 面:  函数类型限定符,用来确定函数是在 CPU 还是在 GPU 上执行,以及这个函数是从 CPU 调用还是从 GPU 调用。  __device__ , 表示从 GPU 上调用,在 GPU 上执行;也就是说其可以被 __global__ 或者__device__修饰的函数调用。此限定符修饰的函数使用有限 制,比如在 G80/GT200 架构上不能使用递归,不能使用函数指针等,具体可参 见我翻译的编程指南。  __global__,表示在 CPU 上调用,在 GPU 上执行,也就是所谓的内核(kernel) 函数;内核只能够被主机调用,内核并不是一个完整的程序,它只是一个数据 并行步骤,其指令流由多个线程执行。  __host__,表示在 CPU 上调用,在 CPU 上执行,这是默认时的情况,也就是传 统的 C 函数。CUDA 支持__host__和__device__的联用,表示同时为主机和设备 编译。  变量类型限定符,用来规定变量存储什么位置上。在传统的 CPU 程序上,这个任务 由编译器承担。在 CUDA 中,不仅要使用主机端的内存,还要使用设备端的显存和 GPU 片上的寄存器、共享存储器和缓存。在 CUDA 存储器模型中,一共抽象出来了 8 种不同的存储器。复杂的存储器模型使得必须要使用限定符要说明变量的存储位 置。  __device__,__device__表明声明的数据存放在显存中,所有的线程都可以访 问,而且主机也可以通过运行时库访问;  __shared__,__shared__表示数据存放在共享存储器在,只有在所在的块内的 线程可以访问,其它块内的线程不能访问;  __constant__,__constant__表明数据存放在常量存储器中,可以被所有的线 程访问,也可以被主机通过运行时库访问;  Texture,texture 表明其绑定的数据可以被纹理缓存加速存取,其实数据本身 的存放位置并没有改变,纹理是来源于图形学的一介概念,CUDA 使用它的原因 一部分在于支持图形处理,另一方面也可以利用它的一些特殊功能。  如果在 GPU 上执行的函数内部的变量没有限定符,那表示它存放在寄存器或者 本地存储器中,在寄存器中的数据只归线程所有,其它线程不可见。  如果 SM 的寄存器用完,那么编译器就会将本应放到寄存器中的变量放到本地存 储器中。  执行配置运算符<<< >>>,用来传递内核函数的执行参数。执行配置有四个参数,第 一个参数声明网格的大小,第二个参数声明块的大小,第三个参数声明动态分配的 共享存储器大小,默认为 0,最后一个参数声明执行的流,默认为 0。  五个内建变量,用于在运行时获得网格和块的尺寸及线程索引等信息  gridDim, gridDim 是一个包含三个元素 x,y,z 的结构体,分别表示网格在 x,y,z 三个方向上的尺寸,虽然其有三维,但是目前只能使用二维;  blockDim, blockDim 也是一个包含三个元素 x,y,z 的结构体,分别表示块在 6 x,y,z 三个方向上的尺寸,对应于执行配置中的第一个参数,对应于执行配置 的第二个参数;  blockIdx, blockIdx 也是一个包含三个元素 x,y,z 的结构体,分别表示当前线 程所在块在网格中 x,y,z 三个方向上的索引;  threadIdx, threadIdx 也是一个包含三个元素 x,y,z 的结构体,分别表示当前 线程在其所在块中 x,y,z 三个方向上的索引;  warpSize,warpSize 表明 warp 的尺寸,在计算能力为 1.0 的设备中,这个值是 24,在 1.0 以上的设备中,这个值是 32。 其它的还有数学函数,原子函数,纹理读取、绑定函数,内建栅栏,内存 fence 函数 等。一般而言,知道这些就应该能够写出 CUDA 程序了,当然要写好的话,必须知道很多其 它的细节。 第五节、计算π 上一节,已经简单的说了一下 CUDA C 的基本语法;因而在本节,兄弟决定以一个例子 为基础说明 CUDA 程序的基本组成部分,不过说实话兄弟选择的例子并不太好,这个例子就 是采用积分法计算圆周率/π的值。其计算原理是:在[0,1]范围内积分 1/(1+x*x) 首先,让我们看一下在 CPU 上的计算流程,其计算流程如下 /* 串行计算 PI 的程序,基本思想为:将积分区间均分为 num 小块,将每小块的面积加起来。 */ float cpuPI(int num){ float sum=0.0f; float temp; for(int i=0;i>1);i>0;i>>=1){ if(threadIdx.x>1);i>0;i>>=1){ if(id vi; Int len = vi.size(); For(int I = 0; I < len; i++) … 代码二 vector vi; … for(int i = 0; I < vi.size(); i++) …. 27 最后要说的是,很多时候编译器能够帮助我们优化,因此,这里说的优化方法在某些编译器下,或者某 些编译器开启了某些选项的情况下,并不有效。 第五节、CUDA程序的优化 对于一个问题来说,往往有很多办法可以解决它。选择那个方法,这是一个问题,而且 是一个与政策、技术,甚至阴谋相关的问题,具体怎么做,往往考量很多。或者说条件不一 样,评判的标准不一样,选择的方法也就不一样。编程也类似,一般来说,不是为了特地使 用 GPU,就没有必要没有原因就将原来的CPU代码修改为GPU。 说这段话,是因为发现有很多人在不适合 GPU 的情况下,却使用它,结果当然很坏具 了。 使用CUDA可分两种情况,一种是从头到尾从新设计;一种是已经有了串行的版本, 但不能满足速度要求,想通过使用CUDA来加快速度,此时可以考虑将代码中最耗时的部 分使用 CUDA 重写。 1 选择合适的工具 由于历史的原因,我们目前所使用的串行编程模式是完全一致的,说白了就是以控制为 中心,虽然面向对象说是以数据为中心,但是我认为其本质上还是以控制为中心,大家想想 看你写面向对象的代码时,其整个中心是不是处理逻辑?是不是系统要做什么?而目前的并 行编程模式却分为数据并行和任务并行,或者分布式存储并行和共享存储并行,CUDA是 基于数据并行和共享存储器的并行,如果你有 OpenMP 或者 MPI 的编程经验,相信你也许也 会有这种看法。 所以如果你要解决的问题并不是数据并行的或者共享存储器的,使用 CUDA 就绝不是一 个好的选择。 2 优化最应该优化的地方-优化准则 相信很多人都知道 80%-20%定律,也就是说程序的 80%时间花在执行 20%的代码上,如果 我们将那 20%的代码加速了10倍,最终也才加速一点点,如果我们将那80%加速了一 倍,最终的结果将加速近一倍。 根据公式 s = 1/(f + (1-f)/p),其中f是程序中串行代码的比重,p是处理器数目。可 知f越大,s就越小。这也许解释了nv的说法(尽量让更多的代码在 GPU 上执行,就算没 有加速比)。也就是尽量提高并行代码的比重。这潜在的能够提升程序的并行效率,Fermi 的某些功能显然是向前迈进了一些,但是nv是不是真的能够达到他们的设想,那就是另一 回事了。 3 全局存储器优化 在 CUDA 的存储器结构中,全局存储器是最慢的,其延迟在几百个时钟,但是其容量是 最大的,我们使用 cudaMalloc 函数分配的指针都是指向全局存储器空间的,我们不可能避 开它们的使用。这两个因素使得全局存储器的使用对 CUDA 程序的性能至关重要。 3.1 合并访问 为了提高访问全局存储器的效率,CUDA 提供了一种称为合并访问的机制来加速全局存 储器的访问。在计算能力 1.0 和 1.1 的设备上,合并访问能够在一次存储器访问中访问最多 达 16 个数据;在计算能力 1.2 和 1.3 的设备上,合并访问一次最多能够访问 32 个数据。另 外线程间的切换也能够隐藏全局存储器的访问延迟。 所谓合并访问就是尽量要求相邻的线程访问相邻的地址空间,当然在不同计算能力的设 28 备上,其具体要求有所差别,在1.0和1.1的设备上,合并访问要求第k个线程访问对 齐的第k个存储地址,必须对齐,不能交叉,而1.2以上的设备没有此要求。 对于不同的数据类型,合并访问的带宽不相同,一般而言使用四字节的数据类型比较 好,但是并不绝对,在某些情况下,使用8字节或16字节的数据性能可能会更好,尤其是 在满足1.0和1.1设备的合并访问条件下。 在设计的时候要考虑到全局存储器的合并访问,并以此为要求来设计,比如在必要的时 候对矩阵进行转置。对数据的组织方式加以改变。 4 共享存储器优化 共享存储器比较靠近SM的,因此其速度相当的快,一般而言可以在一到二个时钟内读 写,因此使用共享存储器取代全局存储器会极大的节约带宽,但是共享存储器的使用是有要 求,它要求数据有局部性(一个block内共享)重用。当然有时也可用共享存储器取代 寄存器,这会在后面说到。 一般的使用共享存储器的方式是:s_x[tid] = x[id];其中 tid 是块内线程的块内索 引,也就是类似 threadIdx,id 是线程的网格内索引。然后可以访问s_x取代访问x。这 样一般可以提升程序的效率。 要注意的是,使用共享存储器往往伴随着存储器的一致性问题,此时可求助于 __syncthreads()和 memory fence 大神。 4.1 存储体冲突 类似于合并访问,共享存储器一次也能满足 半束线程的要求,条件是这半束线程访问 的数据在不同的存储体中。CUDA 将共享存储器组织成 16 列,每一列称为一个存储体。如果 有两个或以上的线程访问的数据在同一个存储体中,此时不能在一次存储体访问中,满足半 束线程的要求,这称为存储体冲突。一般而言,存储体冲突对性能的影响还不是太大,当然 具体的影响要看半束内最多有多少线程落入同一存储体内。要提到的 编程时要注意的是对于一个SM来说,共享存储器的数量是有限的,如果超出其使用限 制,其结果未知,我的经验告诉我,有时没有问题,有时会有大问题。另外,使用共享存储 器有时可能会影响程序的可扩展性,这个在编写库代码的时候要特别注意。 5 常量存储器优化 CUDA 允许分配最多 64KB的常量存储器,常量存储器顾名思义内容是不变的,所以也 有人称其为不变存储器。每个SM有6-8KB的常量缓存,一般而言一到两个周期可读取 常量存储器。如果半束内的线程访问的不是同一个地址,那么各个线程的访问将会串行化。 常量存储器的设计个人认为不是太好,有点像鸡肋。因为要求常量存储器是全局的,因 此对程序的可读性和可扩展性都有影响。 对于常量存储器大小的限制问题,有些情况下,可以使用多次导入并多次执行内核的方 式解决。 6 纹理存储器 纹理存储器是来自图形学的一个概念。由于硬件的支持,提供了很多额外处理功能,如 边界处理、滤波等。访问纹理存储器要通过纹理参考和纹理获取。 CUDA对于随机存取的性能极其的悲剧,而纹理存储器可减弱这种悲剧的效果,尤其 是访问的数据之间具有极大的局部性的时候。 另外在某些情况下可以利用纹理的插值功能来加速计算。此时要注意的纹理执行的是低 精度的线性插值。 29 7 寄存器 寄存器是 NVIDIA GPU 存储器空间中,最快的硬件,读取它几乎不用耗费时间。因此合 理的使用它是至关重要的。 7.1 寄存器溢出和本地存储器 由于GPU上寄存器的数目是有限的,在计算能力 1.0,1.1 的设备上,其数目是 8 K, 在计算能力 1.2,1.3 的设备上,其值是 16 K。如果我们使用的寄存器数目超过了系统的最 大数目,编译器便会将其转入设备存储器上的一个区域内,这个区域称为本地存储器。本地 存储器的速度和全局存储器一样。 使用 2.3 及以前的 toolkit,可以使用-keep 指令输出中间文件,然后直接查看 cubin 文 件就可以知道内核使用的寄存器数目,同时也能看到共享存储器、常量存储器和本地存储器 的使用量。如果使用的是 3.0 的 toolkit,此时 cubin 文件已经不是文本格式,因此不能直 接查看,此时可以在用 nvcc 编译的时候,使用—ptxas-options=-v 选项,编译器会报告存 储器使用信息。 8 执行配置和占用率 使用<<>>语法指定执行线程配置的时候,grid 和 block 大小也影响程序 的效率。一般而言,grid 要大于多处理器的数目,这样才能让多处理器不至于空闲,但是 这样也会导致一些问题,比如负载均衡,如果 grid 大小不能比 sm 数目整除的话,就会有 SM计算的时候,另外一些SM空闲,如果 grid 远大于sm数目的话,可忽略,但是如果 SM数目与grid大小相差不大的话,性能损耗就很可观了。一般而言,grid大小至 少是sm数目的三倍,如果有块内同步的话,grid要大于sm的四倍以上。一般而言, block大小要是束大小的四倍以上,此时基本上可隐藏访存延迟。如果数据量比较小的 话,grid大小和block大小可能相互牵制,此时要综合考虑。在数据量比较大的时 候,只要考虑block大小就行了。 另外,CUDA 还有一个占用率问题,所谓占用率就是SM上活跃块数目和最大允许数目的 比例。使用cuda visual profiler的时候会有这个值。这个值的影响 因素有内核内寄存器、共享存储器的用量。由于 CUDA 线程在切换的时候并不转储线程的状 态,这是CUDA线程极轻的原因,但是这也使得为了保证线程能够切换,内核内使用的资 源最多只能是最大资源的一半,此时占用率是 0.25,一般而言占用率要在 0.5 左右。 9 指令优化 各种指令有不同的执行时间,这样要更快的速度可尽量使用占用时间少的指令,比如 x *= 2; 最好写成 x += x;因为乘法的占用的时间比加法多,但是这个并不明显,因为在现代 处理器上乘法的效率已经和加法差不多了。要注意的是除法和模余,它们占用的时间远大于 加法,要尽量少使用。另外一些超越函数占用的时间更长,为此 CUDA 提供了超越函数的更 快版本,更快的代价是精度损失。所以在速度远比精度重要的情况下,可以使用,但是如果 精度比较重要的话,就要注意了。 9.1 寄存器依赖 在程序的一些部分,我们经常要存取数据,如果在后面的比较近的指令中要读取这个数 据,此时就要等到这个数据存取完成,这就是寄存器依赖,其后果就是由于要等待导致的性 能损耗。这种损耗往往可以通过重新安排程序的语句减弱。 30 10 分支优化 在 CUDA 中,分支会极大的减弱性能,因为SM没有分支预测,因此只能让束内线程在 每个分支上都执行一遍,当然如果某个分支没有线程执行,就可以忽略,因此要减少分支的 数目。但是在实践中很多时候分支是没有办法减少的。此时目前还没有好的解决方案。 第三章、一些例子 第一节、 两向量的距离 本程序中,我以 a[N]和 b[N]代表两个向量,其欧氏距离计算的串行C代码如下: 其中 dis 表示两个向量的距离,sqrt 表示求平方根函数。从代码可以看出,我们可以使用N 个线程,每个线程计算 a,b 一维上数据差的平方,但是由于这个计算量太小,为了加大计算 量,我采用了一个线程计算多个数据差的平方和,在这步后,我们要将所有线程数据加起来 以得到最终的结果,考虑到CUDA不支持全局同步而支持块内同步,因此我们可以先得到 每个块的和,然后再将每个块的和保存到全局存储器 d_temp 中。最后再重新建立一个内核 来将d_temp的数据求和,由于此时数据量小,因此只要一个块就够了。 考虑到分支会极大的影响到 CUDA 程序的性能,因此我们将块内各线程数据加和的方式 是折半相加,此时可以利用高速的共享存储器。折半相加的思想是:假设有 n = 2^k 个线 程,那么第一次相加在 0 和 n/2, 1 和 n/2+1,……,n/2-1 和 n-1 之间进行。此后线程 0 至 n/2- 1 拥有的数据和等于之前块内所有数据和,然后可以进行下一次折半相加。在再次折半相加 之间为了保证数据一致性,必须使用__syncthreads()函数在块内所有线程间同步。具体代码 如下: dis = 0; for(int i = 0; i < N; i++){ dis += (a[i] - b[i]) * (a[i - b[i]); } dis = sqrt(dis); 31 为了使代码能够同时应用于整形和浮点型,使用了模板机制。下面来具体分析: 1.取得当前线程的索引,CUDA 线程有两种组织,一是网格,一是块。 2.取得整个网格内的线程数。 3.动态声明共享存储器,其大小由内核调用的第三个参数决定,单位是字节。 4.初始化共享存储器为 0。 5.-6.声明两个寄存器,并初始化一个为 0。 7.-10.遍历两个向量的所有维并将向量对应维差的平方加到对应线程的寄存器上,为了保 证全局存储器的合并访问,采用了循环读取模式,循环读取是指,假设一个有5个线程,1 0个数,那个线程0就读第0和第5个数,线程1读第1个和第六个数,以此类推。 11.-12.将块内线程取得的对应维差的平方和存入共享存储器并同步,同步的目的是防止 有些线程可能先执行 15.语句而取得了不正确结果。 13.-17.将块内共享存储器的数据全加到共享存储器数组的第0个元素上。 18. 将整个块的数据和存入全局存储器。为什么要做这一步的原因在于:共享存储器并 不持久,当内核执行完后,其就会被回收,而全局存储器是持久的。这类似于CPU中内存 和栈。 在获得每个块的数据和并存入全局存储器后,下一步要做的就是将上一内核存入全局存 储器的数据加起来得到最终的结果。由于内容和上一内核函数差不多,因此就不详细说了。 template static __global__ void reduceMultiBlock(const T* __restrict__ d_a, const T* __restrict__ d_b, const int dim, T* __restrict__ d_temp){ //global thread id 1. unsigned int id = blockDim.x*blockIdx.x + threadIdx.x; //iterate length 2. unsigned int iter = blockDim.x*gridDim.x; //for temp distance, length = blockDim.x 3. extern __shared__ T s_dis[]; //use for when dim is shorter than iter 4. s_dis[threadIdx.x] = 0; 5. T temp; 6. T tempDis = (T)0; 7. for(int i = id; i < dim; i += iter){ 8. temp = d_a[i] - d_b[i]; 9. tempDis += temp*temp; 10. } 11. s_dis[threadIdx.x] = tempDis; 12. __syncthreads(); 13. for(int i = (blockDim.x>>1); i > 0; i >>= 1){ 14. if(threadIdx.x < i) 15. s_dis[threadIdx.x] += s_dis[i + threadIdx.x]; 16. __syncthreads(); 17. } 18. d_temp[blockIdx.x] = s_dis[0]; } 32 代码如下: 经过这两个内核的处理,d_temp[0]就是向量所有维差的平方和了。我们只要将其从显存 上复制到内存并求其平方根就得到了这两个向量的距离。 第二节、 矩阵与向量乘积 本文使用一个线程块计算矩阵的一行和向量的乘积的模式。主要代码如下: template static __global__ void reduceOneBlock(T* __restrict__ d_temp, const int len){ unsigned int tid = threadIdx.x; T temp = (T)0; //length == blockDim.x extern __shared__ T s_dis[]; //use for when len is shorter than d_temp's length s_dis[tid] = (T)0; for(int i = tid; i < len; i += blockDim.x){ temp += d_temp[i]; } s_dis[tid] = temp; __syncthreads(); for(int i = (blockDim.x>>1); i > 0; i >>= 1){ if(tid < i) s_dis[tid] += s_dis[i + tid]; __syncthreads(); } if(0 == tid) *d_temp = s_dis[0]; } 33 这种做法的主要思想还是 reduction,在矩阵的列数和行数都比较大的时候,这种方法 的性能很好,但是如果矩阵是那种行比较多,列比较少的话,应该怎么做?我们可以使用一 个 warp 计算一行(有人称之为 warp 模式),代码如下所示: void __global__ matrixMultiVectorBlock(const int rowSize, const int columnSize, const int pitchItem, const float* __restrict__ d_matrix, const float* __restrict__ d_vec, float* __restrict__ d_r){ unsigned int tid = threadIdx.x; extern __shared__ float s_r[]; float temp = 0.0f; for(int i = tid; i < columnSize; i += blockDim.x){ temp += d_matrix[blockIdx.x*pitchItem+i]*d_vec[i]; } s_r[tid] = temp; __syncthreads(); for(int i = (blockDim.x>>1); i > 32; i >>= 1){ if(tid < i){ s_r[tid] += s_r[tid+i]; } __syncthreads(); } if(tid < 32){ s_r[tid] += s_r[tid+32]; } if(tid < 16){ s_r[tid] += s_r[tid+16]; } if(tid < 8){ s_r[tid] += s_r[tid+8]; } if(tid < 4){ s_r[tid] += s_r[tid+4]; } if(tid < 2){ s_r[tid] += s_r[tid+2]; } if(tid < 1){ s_r[tid] += s_r[tid+1]; d_r[blockIdx.x] = s_r[0]; } } 34 void __global__ matrixMultiVectorWarp(const int rowSize, const int columnSize, const int pitchItem, const float* __restrict__ d_matrix, const float* __restrict__ d_vec, float* __restrict__ d_r ){ unsigned int warpid = threadIdx.y + blockIdx.x*blockDim.y; //length equal to two block size extern __shared__ float s_s[]; float *s_v = s_s; unsigned int blockSize = blockDim.y*blockDim.x; float *s_r = s_s + blockSize; unsigned int tid = threadIdx.x + threadIdx.y*blockDim.x; s_r[tid] = 0.0f; int iterTimes = (columnSize+blockSize-1)/blockSize; for(int i = 0; i < iterTimes ; i++){ if(i*blockSize+tid < columnSize){ s_v[tid] = d_vec[i*blockSize+tid]; } __syncthreads(); int temp = min(blockSize, columnSize-i*blockSize); for(int j = threadIdx.x; j < temp; j += blockDim.x){ s_r[tid] += d_matrix[i*blockSize+j+pitchItem*warpid]*s_v[j]; } __syncthreads(); } if(warpid >= rowSize){ return; } if(threadIdx.x < 16){ s_r[tid] += s_r[tid+16]; } if(threadIdx.x < 8){ s_r[tid] += s_r[tid+8]; } if(threadIdx.x < 4){ s_r[tid] += s_r[tid+4]; } if(threadIdx.x < 2){ s_r[tid] += s_r[tid+2]; } if(threadIdx.x < 1){ s_r[tid] += s_r[tid+1]; } if(threadIdx.x == 0){ d_r[warpid] = s_r[threadIdx.y<<5]; } } 35 由于如果没有 host 端代码的话,上面的代码非常难以理解,因此本文提供 host 端的代码以 供参考。 不知道这些代码给了你灵感没有,我第一次设计这些的代码的时候,非常的兴奋,在 GTX295 上的测试表明在每行元素比较多的话,warp 模式没有优势,而在每行元素比较少的 时候,warp 可比 block 模式快一倍以上。如果你知道更好的计算方式,请你联系我,我的 qq:304128534 第三节、 线性方程组的求解 本文的线性方程组求解主要采用主元素高斯消元法。其分为三个步骤:一、选主元; 二、交换主元素行与当前行;三、消元,消元时也消除对角线上面的元素。 一、选主元:选主元的做法就是对于当前列,先用一个内核将每一个线程块中的绝对值 最大的数据求出来并保存其下标,将这两数据存到全局存储器。然后再另开一内核求得最大 值所在的位置并保存其下标。 void runMatrixMultiVectorWarpGPU(const int rowSize, const int columnSize, const float *matrix, const float *v, float *r){ float *d_matrix; size_t pitch; cutilSafeCall(cudaMallocPitch((void**)&d_matrix, &pitch, columnSize*sizeof(float), rowSize)); cutilSafeCall(cudaMemcpy2D(d_matrix, pitch, matrix, columnSize*sizeof(float), columnSize*sizeof(float), rowSize, cudaMemcpyHostToDevice)); float *d_v; cutilSafeCall(cudaMalloc((void**)&d_v, columnSize*sizeof(float))); cutilSafeCall(cudaMemcpy(d_v, v, columnSize*sizeof(float), cudaMemcpyHostToDevice)); float *d_r; cutilSafeCall(cudaMalloc((void**)&d_r, rowSize*sizeof(float))); for(int i = 0; i < 1000; i++){ dim3 threads(32, 8); matrixMultiVectorWarp<<<(rowSize+threads.y-1)/threads.y, threads, 2*sizeof(float)*threads.x*threads.y>>>(rowSize, columnSize, pitch/sizeof(float), d_matrix, d_v, d_r); cutilSafeCall(cudaThreadSynchronize()); } cutilSafeCall(cudaMemcpy(r, d_r, rowSize*sizeof(float), cudaMemcpyDeviceToHost)); cutilSafeCall(cudaFree(d_r)); cutilSafeCall(cudaFree(d_v)); cutilSafeCall(cudaFree(d_matrix)); } 36 //find the main row static __global__ void mainRowMultiBlock(const int dim, const int currentRow, float *d_matrix, float *d_matrixTemp, int *d_mainRowTemp){ const unsigned int id = blockDim.x*blockIdx.x + threadIdx.x; extern __shared__ float s_s[];//len 2*blockDim float *s_m = s_s; int *s_id = (int*)(s_s + blockDim.x); float r = 0.0f; int index = id + currentRow; int rId = 0; for(int i = index; i < dim; i += gridDim.x * blockDim.x){ if( fabsf(d_matrix[i*(dim+1) + currentRow]) > r){ r = fabsf(d_matrix[i*(dim+1) + currentRow]); rId = i; } } s_m[threadIdx.x] = r; s_id[threadIdx.x] = rId; __syncthreads(); for(int i = (blockDim.x/2); i > 0; i /= 2){ if((threadIdx.x < i) && (s_m[threadIdx.x + i] > s_m[threadIdx.x])){ s_m[threadIdx.x] = s_m[threadIdx.x + i]; s_id[threadIdx.x] = s_id[i + threadIdx.x]; } __syncthreads(); } if(0 == threadIdx.x){ d_matrixTemp[blockIdx.x] = s_m[0]; d_mainRowTemp[blockIdx.x] = s_id[0]; } } static __global__ void mainRowOneBlock(int len, float *d_matrixTemp, int *d_mainRowTemp){ const unsigned int tid = threadIdx.x; extern __shared__ float s_s[]; float *s_m = s_s; int *s_id = (int*)(s_s + blockDim.x); float r = 0.0f; int rId = 0; for(int i = tid; i < len; i += blockDim.x){ if(r < d_matrixTemp[i]){ r = d_matrixTemp[i]; rId = d_mainRowTemp[i]; } } s_m[tid] = r; s_id[tid] = rId; // __syncthreads(); for(int i = (blockDim.x>>1); i > 0; i >>= 1){ if(tid < i && s_m[tid + i] > s_m[tid]){ s_m[tid] = s_m[tid + i]; s_id[tid] = s_id[tid + i]; } // __syncthreads(); } if(0 == tid){ d_mainRowTemp[0] = s_id[0]; } } 37 二、交换主元素行,代码如下: 三、消元,代码如下: 主机端代码如下: //swap the main row and current row static __global__ void swapMain(const int dim, const int currentRow, const int *mainRow, float *d_matrix){ const unsigned int id = blockDim.x*blockIdx.x + threadIdx.x; int row = *mainRow; for(int i = id + currentRow; i < dim+1; i += blockDim.x*gridDim.x){ float temp = d_matrix[currentRow*(dim + 1) + i]; d_matrix[currentRow*(dim + 1) + i] = d_matrix[row*(dim + 1) + i]; d_matrix[row*(dim + 1) + i] = temp; } } static __global__ void reduceRow(const int dim, const int currentRow, float *d_matrix){ const unsigned bid = blockIdx.x; const unsigned tid = threadIdx.x; if(bid == currentRow) return; __shared__ float fac; fac = d_matrix[bid*(dim+1) + currentRow]/d_matrix[currentRow*(dim+1) + currentRow]; __syncthreads(); for(int k = tid + currentRow; k < dim+1; k += blockDim.x){ d_matrix[bid*(dim+1) + k] -= d_matrix[currentRow*(dim+1) + k]*fac; } } void solveFuncGPU(const int dim, float *d_matrix){ unsigned int grid = 30*4; unsigned int block; int *d_mainRowTemp; cutilSafeCall(cudaMalloc((void**)&d_mainRowTemp, grid*sizeof(int))); float *d_matrixTemp; cutilSafeCall(cudaMalloc((void**)&d_matrixTemp, grid*sizeof(float))); for(int i = 0; i < dim; i++){ block = 256; if(dim - i < 256) block = 32; mainRowMultiBlock<<>>(dim, i, d_matrix, d_matrixTemp, d_mainRowTemp); mainRowOneBlock<<<1, 32, 64*sizeof(int)>>>(grid, d_matrixTemp, d_mainRowTemp); swapMain<<>>(dim, i, d_mainRowTemp, d_matrix); reduceRow<<>>(dim , i, d_matrix); } cutilSafeCall(cudaFree(d_matrixTemp)); cutilSafeCall(cudaFree(d_mainRowTemp)); } 38 很明显的可以看出来,这个程序中最消耗时间的应该是消元,而消元的内核消耗依旧有 优化的余地,原因在于对全局存储器的访问是合并的,但是并不完美(每次访问的首地址可 能没有对齐),想到这一点,你应该能够想到优化的办法了,优化后性能提高了差不多一 倍。
还剩37页未读

继续阅读

下载pdf到电脑,查找使用更方便

pdf的实际排版效果,会与网站的显示效果略有不同!!

需要 8 金币 [ 分享pdf获得金币 ] 1 人已下载

下载pdf

pdf贡献者

goat

贡献于2014-09-24

下载需要 8 金币 [金币充值 ]
亲,您也可以通过 分享原创pdf 来获得金币奖励!
下载pdf