深入浅出谈CUDA
发表时间：2008-11-21
“CUDA是NVIDIA的GPGPU模型，它使用C语言为基础，可以直接以大多数人熟悉的C语言，写出在显示芯片上执行的程序，而不需要去学习特定的显示芯片的指令或是特殊的结构。”
“CUDA是NVIDIA的GPGPU模型，它使用C语言为基础，可以直接以大多数人熟悉的C语言，写出在显示芯片上执行的程序，而不需要去学习特定的显示芯片的指令或是特殊的结构。”
CUDA是什么？能吃吗？
编者注：NVIDIA的GeFoce8800GTX发布后，它的通用计算架构CUDA经过一年多的推广后，现在已经在有相当多的论文发表，在商业应用软件等方面也初步出现了视频编解码、金融、地质勘探、科学计算等领域的产品，是时候让我们对其作更深一步的了解。为了让大家更容易了解CUDA，我们征得Hotball的本人同意，发表他最近亲自撰写的本文。这篇文章的特点是深入浅出，也包含了hotball本人编写一些简单CUDA程序的亲身体验，对于希望了解CUDA的读者来说是非常不错的入门文章，PCINLIFE对本文的发表没有作任何的删减，主要是把一些台湾的词汇转换成大陆的词汇以及作了若干"编者注"的注释。
现代的显示芯片已经具有高度的可程序化能力，由于显示芯片通常具有相当高的内存带宽，以及大量的执行单元，因此开始有利用显示芯片来帮助进行一些计算工作的想法，即GPGPU。CUDA即是NVIDIA的GPGPU模型。
NVIDIA的新一代显示芯片，包括GeForce8系列及更新的显示芯片都支持CUDA。NVIDIA免费提供CUDA的开发工具（包括Windows版本和Linux版本）、程序范例、文件等等，可以在CUDAZone下载。
GPGPU的优缺点 使用显示芯片来进行运算工作，和使用CPU相比，主要有几个好处：
显示芯片通常具有更大的内存带宽。例如，NVIDIA的GeForce8800GTX具有超过50GB/s的内存带宽，而目前高阶CPU的内存带宽则在10GB/s左右。
显示芯片具有更大量的执行单元。例如GeForce8800GTX具有128个"streamprocessors"，频率为1.35GHz。CPU频率通常较高，但是执行单元的数目则要少得多。
和高阶CPU相比，显卡的价格较为低廉。例如目前一张GeForce8800GT包括512MB内存的价格，和一颗2.4GHz四核心CPU的价格相若。
当然，使用显示芯片也有它的一些缺点：
显示芯片的运算单元数量很多，因此对于不能高度并行化的工作，所能带来的帮助就不大。
显示芯片目前通常只支持32bits浮点数，且多半不能完全支持IEEE754规格，有些运算的精确度可能较低。目前许多显示芯片并没有分开的整数运算单元，因此整数运算的效率较差。
显示芯片通常不具有分支预测等复杂的流程控制单元，因此对于具有高度分支的程序，效率会比较差。
目前GPGPU的程序模型仍不成熟，也还没有公认的标准。例如NVIDIA和AMD/ATI就有各自不同的程序模型。
整体来说，显示芯片的性质类似streamprocessor，适合一次进行大量相同的工作。CPU则比较有弹性，能同时进行变化较多的工作。
CUDA架构 CUDA是NVIDIA的GPGPU模型，它使用C语言为基础，可以直接以大多数人熟悉的C语言，写出在显示芯片上执行的程序，而不需要去学习特定的显示芯片的指令或是特殊的结构。
在CUDA的架构下，一个程序分为两个部份：host端和device端。Host端是指在CPU上执行的部份，而device端则是在显示芯片上执行的部份。Device端的程序又称为"kernel"。通常host端程序会将数据准备好后，复制到显卡的内存中，再由显示芯片执行device端程序，完成后再由host端程序将结果从显卡的内存中取回。
由于CPU存取显卡内存时只能透过PCIExpress接口，因此速度较慢（PCIExpressx16的理论带宽是双向各4GB/s），因此不能太常进行这类动作，以免降低效率。
在CUDA架构下，显示芯片执行时的最小单位是thread。数个thread可以组成一个block。一个block中的thread能存取同一块共享的内存，而且可以快速进行同步的动作。
每一个block所能包含的thread数目是有限的。不过，执行相同程序的block，可以组成grid。不同block中的thread无法存取同一个共享的内存，因此无法直接互通或进行同步。因此，不同block中的thread能合作的程度是比较低的。不过，利用这个模式，可以让程序不用担心显示芯片实际上能同时执行的thread数目限制。例如，一个具有很少量执行单元的显示芯片，可能会把各个block中的thread顺序执行，而非同时执行。不同的grid则可以执行不同的程序（即kernel）。
Grid、block和thread的关系，如下图所示：
每个thread都有自己的一份register和localmemory的空间。同一个block中的每个thread则有共享的一份sharememory。此外，所有的thread（包括不同block的thread）都共享一份globalmemory、constantmemory、和texturememory。不同的grid则有各自的globalmemory、constantmemory和texturememory。这些不同的内存的差别，会在之后讨论。
执行模式 由于显示芯片大量并行计算的特性，它处理一些问题的方式，和一般CPU是不同的。主要的特点包括：
内存存取latency的问题：CPU通常使用cache来减少存取主内存的次数，以避免内存latency影响到执行效率。显示芯片则多半没有cache（或很小），而利用并行化执行的方式来隐藏内存的latency（即，当第一个thread需要等待内存读取结果时，则开始执行第二个thread，依此类推）。
分支指令的问题：CPU通常利用分支预测等方式来减少分支指令造成的pipelinebubble。显示芯片则多半使用类似处理内存latency的方式。不过，通常显示芯片处理分支的效率会比较差。
因此，最适合利用CUDA处理的问题，是可以大量并行化的问题，才能有效隐藏内存的latency，并有效利用显示芯片上的大量执行单元。使用CUDA时，同时有上千个thread在执行是很正常的。因此，如果不能大量并行化的问题，使用CUDA就没办法达到最好的效率了。
CUDAToolkit的安装
目前NVIDIA提供的CUDAToolkit（可从这里下载）支持Windows（32bits及64bits版本）及许多不同的Linux版本。
CUDAToolkit需要配合C/C++compiler。在Windows下，目前只支持VisualStudio7.x及VisualStudio8（包括免费的VisualStudioC++2005Express）。VisualStudio6和gcc在Windows下是不支援的。在Linux下则只支援gcc。
这里简单介绍一下在Windows下设定并使用CUDA的方式。
下载及安装 在Windows下，CUDAToolkit和CUDASDK都是由安装程序的形式安装的。CUDAToolkit包括CUDA的基本工具，而CUDASDK则包括许多范例程序以及链接库。基本上要写CUDA的程序，只需要安装CUDAToolkit即可。不过CUDASDK仍值得安装，因为里面的许多范例程序和链接库都相当有用。
CUDAToolkit安装完后，预设会安装在C:\CUDA目录里。其中包括几个目录：
bin--工具程序及动态链接库
doc--文件
include--header檔
lib--链接库档案
open64--基于Open64的CUDAcompiler
src--一些原始码
安装程序也会设定一些环境变量，包括：
CUDA_BIN_PATH--工具程序的目录，默认为C:\CUDA\bin
CUDA_INC_PATH--header文件的目录，默认为C:\CUDA\inc
CUDA_LIB_PATH--链接库文件的目录，默认为C:\CUDA\lib
在VisualStudio中使用CUDA CUDA的主要工具是nvcc，它会执行所需要的程序，将CUDA程序代码编译成执行档(或object檔)。在VisualStudio下，我们透过设定custombuildtool的方式，让VisualStudio会自动执行nvcc。
这里以VisualStudio2005为例：
首先，建立一个Win32Console模式的project（在ApplicationSettings中记得勾选Emptyproject），并新增一个档案，例如main.cu。
在main.cu上右键单击，并选择Properties。点选General，确定Tool的部份是选择CustomBuildTool。
选择CustomBuildStep，在CommandLine使用以下设定：
Release模式："$(CUDA_BIN_PATH)\nvcc.exe"-ccbin"$(VCInstallDir)bin"-c-DWIN32-D_CONSOLE-D_MBCS-Xcompiler/EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT-I"$(CUDA_INC_PATH)"-o$(ConfigurationName)\$(InputName).obj$(InputFileName)
Debug模式："$(CUDA_BIN_PATH)\nvcc.exe"-ccbin"$(VCInstallDir)bin"-c-D_DEBUG-DWIN32-D_CONSOLE-D_MBCS-Xcompiler/EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd-I"$(CUDA_INC_PATH)"-o$(ConfigurationName)\$(InputName).obj$(InputFileName)
如果想要使用软件仿真的模式，可以新增两个额外的设定：
EmuRelease模式："$(CUDA_BIN_PATH)\nvcc.exe"-ccbin"$(VCInstallDir)bin"-deviceemu-c-DWIN32-D_CONSOLE-D_MBCS-Xcompiler/EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT-I"$(CUDA_INC_PATH)"-o$(ConfigurationName)\$(InputName).obj$(InputFileName)
EmuDebug模式："$(CUDA_BIN_PATH)\nvcc.exe"-ccbin"$(VCInstallDir)bin"-deviceemu-c-D_DEBUG-DWIN32-D_CONSOLE-D_MBCS-Xcompiler/EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd-I"$(CUDA_INC_PATH)"-o$(ConfigurationName)\$(InputName).obj$(InputFileName)
对所有的配置文件，在CustomBuildStep的Outputs中加入$(ConfigurationName)\$(InputName).obj。
选择project，右键单击选择Properties，再点选Linker。对所有的配置文件修改以下设定：
General/EnableIncrementalLinking：No
General/AdditionalLibraryDirectories：$(CUDA_LIB_PATH)
Input/AdditionalDependencies：cudart.lib
这样应该就可以直接在VisualStudio的IDE中，编辑CUDA程序后，直接build以及执行程序了。
第一个CUDA程序
CUDA目前有两种不同的API：RuntimeAPI和DriverAPI，两种API各有其适用的范围。由于runtimeAPI较容易使用，一开始我们会以runetimeAPI为主。
CUDA的初始化 首先，先建立一个档案first_cuda.cu。如果是使用VisualStudio的话，则请先按照这里的设定方式设定project。
要使用runtimeAPI的时候，需要includecuda_runtime.h。所以，在程序的最前面，加上
#include#include
接下来是一个InitCUDA函式，会呼叫runtimeAPI中，有关初始化CUDA的功能：
boolInitCUDA(){???intcount;???cudaGetDeviceCount(&count);???if(count==0){??????fprintf(stderr,"Thereisnodevice.\n");??????returnfalse;???}???inti;???for(i=0;i=1){???????????break;?????????}??????}???}???if(i==count){??????fprintf(stderr,"ThereisnodevicesupportingCUDA1.x.\n");??????returnfalse;???}???cudaSetDevice(i);???returntrue;}
这个函式会先呼叫cudaGetDeviceCount函式，取得支持CUDA的装置的数目。如果系统上没有支持CUDA的装置，则它会传回1，而device0会是一个仿真的装置，但不支持CUDA1.0以上的功能。所以，要确定系统上是否有支持CUDA的装置，需要对每个device呼叫cudaGetDeviceProperties函式，取得装置的各项数据，并判断装置支持的CUDA版本（prop.major和prop.minor分别代表装置支持的版本号码，例如1.0则prop.major为1而prop.minor为0）。
透过cudaGetDeviceProperties函式可以取得许多数据，除了装置支持的CUDA版本之外，还有装置的名称、内存的大小、最大的thread数目、执行单元的频率等等。详情可参考NVIDIA的CUDAProgrammingGuide。
在找到支持CUDA1.0以上的装置之后，就可以呼叫cudaSetDevice函式，把它设为目前要使用的装置。
最后是main函式。在main函式中我们直接呼叫刚才的InitCUDA函式，并显示适当的讯息：
intmain(){???if(!InitCUDA()){??????return0;???}???printf("CUDAinitialized.\n");???return0;}
这样就可以利用nvcc来compile这个程序了。使用VisualStudio的话，若按照先前的设定方式，可以直接BuildProject并执行。
nvcc是CUDA的compile工具，它会将.cu檔拆解出在GPU上执行的部份，及在host上执行的部份，并呼叫适当的程序进行compile动作。在GPU执行的部份会透过NVIDIA提供的compiler编译成中介码，而host执行的部份则会透过系统上的C++compiler编译（在Windows上使用VisualC++而在Linux上使用gcc）。
编译后的程序，执行时如果系统上有支持CUDA的装置，应该会显示CUDAinitialized.的讯息，否则会显示相关的错误讯息。
?
利用CUDA进行运算
到目前为止，我们的程序并没有做什么有用的工作。所以，现在我们加入一个简单的动作，就是把一大堆数字，计算出它的平方和。
首先，把程序最前面的include部份改成：
#include#include#include#defineDATA_SIZE1048576intdata[DATA_SIZE];
并加入一个新函式GenerateNumbers：
voidGenerateNumbers(intnumber,intsize){???for(inti=0;i
这个函式会产生一大堆0~9之间的随机数。
要利用CUDA进行计算之前，要先把数据复制到显卡内存中，才能让显示芯片使用。因此，需要取得一块适当大小的显卡内存，再把产生好的数据复制进去。在main函式中加入：
???GenerateNumbers(data,DATA_SIZE);???intgpudata,result;???cudaMalloc((void)&gpudata,sizeof(int)DATA_SIZE);???cudaMalloc((void)&result,sizeof(int));???cudaMemcpy(gpudata,data,sizeof(int)DATA_SIZE,??????cudaMemcpyHostToDevice);
上面这段程序会先呼叫GenerateNumbers产生随机数，并呼叫cudaMalloc取得一块显卡内存（result则是用来存取计算结果，在稍后会用到），并透过cudaMemcpy将产生的随机数复制到显卡内存中。cudaMalloc和cudaMemcpy的用法和一般的malloc及memcpy类似，不过cudaMemcpy则多出一个参数，指示复制内存的方向。在这里因为是从主内存复制到显卡内存，所以使用cudaMemcpyHostToDevice。如果是从显卡内存到主内存，则使用cudaMemcpyDeviceToHost。这在之后会用到。
接下来是要写在显示芯片上执行的程序。在CUDA中，在函式前面加上__global__表示这个函式是要在显示芯片上执行的。因此，加入以下的函式：
__global__staticvoidsumOfSquares(intnum,intresult){???intsum=0;???inti;???for(i=0;i
在显示芯片上执行的程序有一些限制，例如它不能有传回值。其它的限制会在之后提到。
接下来是要让CUDA执行这个函式。在CUDA中，要执行一个函式，使用以下的语法：
???函式名称<<>>(参数...);
呼叫完后，还要把结果从显示芯片复制回主内存上。在main函式中加入以下的程序：
???sumOfSquares<<<1,1,0>>>(gpudata,result);???intsum;???cudaMemcpy(&sum,result,sizeof(int),cudaMemcpyDeviceToHost);???cudaFree(gpudata);???cudaFree(result);???printf("sum:%d\n",sum);
因为这个程序只使用一个thread，所以block数目、thread数目都是1。我们也没有使用到任何sharedmemory，所以设为0。编译后执行，应该可以看到执行的结果。
为了确定执行的结果正确，我们可以加上一段以CPU执行的程序代码，来验证结果：
???sum=0;???for(inti=0;i
编译后执行，确认两个结果相同。
?
计算运行时间 ?
CUDA提供了一个clock函式，可以取得目前的timestamp，很适合用来判断一段程序执行所花费的时间（单位为GPU执行单元的频率）。这对程序的优化也相当有用。要在我们的程序中记录时间，把sumOfSquares函式改成：
?
__global__staticvoidsumOfSquares(intnum,intresult,???clock_ttime){???intsum=0;???inti;???clock_tstart=clock();???for(i=0;i
把main函式中间部份改成：
???intgpudata,result;???clock_ttime;???cudaMalloc((void)&gpudata,sizeof(int)DATA_SIZE);???cudaMalloc((void)&result,sizeof(int));???cudaMalloc((void)&time,sizeof(clock_t));???cudaMemcpy(gpudata,data,sizeof(int)DATA_SIZE,??????cudaMemcpyHostToDevice);
???sumOfSquares<<<1,1,0>>>(gpudata,result,time);???intsum;???clock_ttime_used;???cudaMemcpy(&sum,result,sizeof(int),cudaMemcpyDeviceToHost);???cudaMemcpy(&time_used,time,sizeof(clock_t),??????cudaMemcpyDeviceToHost);???cudaFree(gpudata);???cudaFree(result);???printf("sum:%dtime:%d\n",sum,time_used);
编译后执行，就可以看到执行所花费的时间了。
如果计算实际运行时间的话，可能会注意到它的执行效率并不好。这是因为我们的程序并没有利用到CUDA的主要的优势，即并行化执行。在下一段文章中，会讨论如何进行优化的动作。
改良第一个CUDA程序
在上一篇文章中，我们做了一个计算一大堆数字的平方和的程序。不过，我们也提到这个程序的执行效率并不理想。当然，实际上来说，如果只是要做计算平方和的动作，用CPU做会比用GPU快得多。这是因为平方和的计算并不需要太多运算能力，所以几乎都是被内存带宽所限制。因此，光是把数据复制到显卡内存上的这个动作，所需要的时间，可能已经和直接在CPU上进行计算差不多了。
不过，如果进行平方和的计算，只是一个更复杂的计算过程的一部份的话，那么当然在GPU上计算还是有它的好处的。而且，如果数据已经在显卡内存上（例如在GPU上透过某种算法产生），那么，使用GPU进行这样的运算，还是会比较快的。
刚才也提到了，由于这个计算的主要瓶颈是内存带宽，所以，理论上显卡的内存带宽是相当大的。这里我们就来看看，倒底我们的第一个程序，能利用到多少内存带宽。
?
程序的并行化
我们的第一个程序，并没有利用到任何并行化的功能。整个程序只有一个thread。在GeForce8800GT上面，在GPU上执行的部份（称为"kernel"）大约花费640M个频率。GeForce8800GT的执行单元的频率是1.5GHz，因此这表示它花费了约0.43秒的时间。1M个32bits数字的数据量是4MB，因此，这个程序实际上使用的内存带宽，只有9.3MB/s左右！这是非常糟糕的表现。
为什么会有这样差的表现呢？这是因为GPU的架构特性所造成的。在CUDA中，一般的数据复制到的显卡内存的部份，称为globalmemory。这些内存是没有cache的，而且，存取globalmemory所需要的时间（即latency）是非常长的，通常是数百个cycles。由于我们的程序只有一个thread，所以每次它读取globalmemory的内容，就要等到实际读取到数据、累加到sum之后，才能进行下一步。这就是为什么它的表现会这么的差。
由于globalmemory并没有cache，所以要避开巨大的latency的方法，就是要利用大量的threads。假设现在有大量的threads在同时执行，那么当一个thread读取内存，开始等待结果的时候，GPU就可以立刻切换到下一个thread，并读取下一个内存位置。因此，理想上当thread的数目够多的时候，就可以完全把globalmemory的巨大latency隐藏起来了。
要怎么把计算平方和的程序并行化呢？最简单的方法，似乎就是把数字分成若干组，把各组数字分别计算平方和后，最后再把每组的和加总起来就可以了。一开始，我们可以把最后加总的动作，由CPU来进行。
首先，在first_cuda.cu中，在#defineDATA_SIZE的后面增加一个#define，设定thread的数目：
#defineDATA_SIZE???1048576#defineTHREAD_NUM??256
接着，把kernel程序改成：
__global__staticvoidsumOfSquares(intnum,intresult,???clock_ttime){???constinttid=threadIdx.x;???constintsize=DATA_SIZE/THREAD_NUM;???intsum=0;???inti;???clock_tstart;???if(tid==0)start=clock();???for(i=tidsize;i<(tid+1)size;i++){??????sum+=num[i]num[i];???}???result[tid]=sum;???if(tid==0)time=clock()-start;}
程序里的threadIdx是CUDA的一个内建的变量，表示目前的thread是第几个thread（由0开始计算）。以我们的例子来说，会有256个threads，所以同时会有256个sumOfSquares函式在执行，但每一个的threadIdx.x则分别会是0~255。利用这个变量，我们就可以让每一份函式执行时，对整个数据不同的部份计算平方和。另外，我们也让计算时间的动作，只在thread0（即threadIdx.x=0的时候）进行。
同样的，由于会有256个计算结果，所以原来存放result的内存位置也要扩大。把main函式中的中间部份改成：
???intgpudata,result;???clock_ttime;???cudaMalloc((void)&gpudata,sizeof(int)DATA_SIZE);???cudaMalloc((void)&result,sizeof(int)THREAD_NUM);???cudaMalloc((void)&time,sizeof(clock_t));???cudaMemcpy(gpudata,data,sizeof(int)DATA_SIZE,??????cudaMemcpyHostToDevice);???sumOfSquares<<<1,THREAD_NUM,0>>>(gpudata,result,time);???intsum[THREAD_NUM];???clock_ttime_used;???cudaMemcpy(&sum,result,sizeof(int)THREAD_NUM,??????cudaMemcpyDeviceToHost);???cudaMemcpy(&time_used,time,sizeof(clock_t),??????cudaMemcpyDeviceToHost);???cudaFree(gpudata);???cudaFree(result);???cudaFree(time);
可以注意到我们在呼叫sumOfSquares函式时，指定THREAD_NUM为thread的数目。最后，在CPU端把计算好的各组数据的平方和进行加总：
???intfinal_sum=0;???for(inti=0;i
编译后执行，确认结果和原来相同。
这个版本的程序，在GeForce8800GT上执行，只需要约8.3Mcycles，比前一版程序快了77倍！这就是透过大量thread来隐藏latency所带来的效果。
不过，如果计算一下它使用的内存带宽，就会发现其实仍不是很理想，大约只有723MB/s而已。这和GeForce8800GT所具有的内存带宽是很大的差距。为什么会这样呢？
内存的存取模式 显卡上的内存是DRAM，因此最有效率的存取方式，是以连续的方式存取。前面的程序，虽然看起来是连续存取内存位置（每个thread对一块连续的数字计算平方和），但是我们要考虑到实际上thread的执行方式。前面提过，当一个thread在等待内存的数据时，GPU会切换到下一个thread。也就是说，实际上执行的顺序是类似
???thread0->thread1->thread2->...
因此，在同一个thread中连续存取内存，在实际执行时反而不是连续了。要让实际执行结果是连续的存取，我们应该要让thread0读取第一个数字，thread1读取第二个数字…依此类推。所以，我们可以把kernel程序改成如下：
__global__staticvoidsumOfSquares(intnum,intresult,???clock_ttime){???constinttid=threadIdx.x;???intsum=0;???inti;???clock_tstart;???if(tid==0)start=clock();???for(i=tid;i
编译后执行，确认结果相同。
仅仅是这样简单的修改，实际执行的效率就有很大的差别。在GeForce8800GT上，上面的程序执行需要的频率是2.6Mcycles，又比前一版程序快了三倍。不过，这样仍只有2.3GB/s的带宽而已。
这是因为我们使用的thread数目还是不够多的原因。理论上256个threads最多只能隐藏256cycles的latency。但是GPU存取globalmemory时的latency可能高达500cycles以上。如果增加thread数目，就可以看到更好的效率。例如，可以把THREAD_NUM改成512。在GeForce8800GT上，这可以让执行花费的时间减少到1.95Mcycles。有些改进，但是仍不够大。不幸的是，目前GeForce8800GT一个block最多只能有512个threads，所以不能再增加了，而且，如果thread数目增加太多，那么在CPU端要做的最后加总工作也会变多。
更多的并行化 前面提到了block。在之前介绍呼叫CUDA函式时，也有提到"block数目"这个参数。到目前为止，我们都只使用一个block。究竟block是什么呢？
在CUDA中，thread是可以分组的，也就是block。一个block中的thread，具有一个共享的sharedmemory，也可以进行同步工作。不同block之间的thread则不行。在我们的程序中，其实不太需要进行thread的同步动作，因此我们可以使用多个block来进一步增加thread的数目。
首先，在#defineDATA_SIZE的地方，改成如下：
#defineDATA_SIZE??1048576#defineBLOCK_NUM??32#defineTHREAD_NUM??256
这表示我们会建立32个blocks，每个blocks有256个threads，总共有32256=8192个threads。
接着，我们把kernel部份改成：
__global__staticvoidsumOfSquares(intnum,intresult,???clock_ttime){???constinttid=threadIdx.x;???constintbid=blockIdx.x;???intsum=0;???inti;???if(tid==0)time[bid]=clock();???for(i=bidTHREAD_NUM+tid;i
blockIdx.x和threadIdx.x一样是CUDA内建的变量，它表示的是目前的block编号。另外，注意到我们把计算时间的方式改成每个block都会记录开始时间及结束时间。
main函式部份，修改成：
???intgpudata,result;???clock_ttime;???cudaMalloc((void)&gpudata,sizeof(int)DATA_SIZE);???cudaMalloc((void)&result,??????sizeof(int)THREAD_NUMBLOCK_NUM);???cudaMalloc((void)&time,sizeof(clock_t)BLOCK_NUM2);???cudaMemcpy(gpudata,data,sizeof(int)DATA_SIZE,??????cudaMemcpyHostToDevice);???sumOfSquares<<>>(gpudata,result,??????time);???intsum[THREAD_NUMBLOCK_NUM];???clock_ttime_used[BLOCK_NUM2];???cudaMemcpy(&sum,result,sizeof(int)THREAD_NUMBLOCK_NUM,??????cudaMemcpyDeviceToHost);???cudaMemcpy(&time_used,time,sizeof(clock_t)BLOCK_NUM2,??????cudaMemcpyDeviceToHost);???cudaFree(gpudata);???cudaFree(result);???cudaFree(time);???intfinal_sum=0;???for(inti=0;itime_used[i])?????????min_start=time_used[i];??????if(max_end
基本上我们只是把result的大小变大，并修改计算时间的方式，把每个block最早的开始时间，和最晚的结束时间相减，取得总运行时间。
这个版本的程序，执行的时间减少很多，在GeForce8800GT上只需要约150Kcycles，相当于40GB/s左右的带宽。不过，它在CPU上执行的部份，需要的时间加长了（因为CPU现在需要加总8192个数字）。为了避免这个问题，我们可以让每个block把自己的每个thread的计算结果进行加总。
Thread的同步 前面提过，一个block内的thread可以有共享的内存，也可以进行同步。我们可以利用这一点，让每个block内的所有thread把自己计算的结果加总起来。把kernel改成如下：
__global__staticvoidsumOfSquares(intnum,intresult,???clock_ttime){???extern__shared__intshared[];???constinttid=threadIdx.x;???constintbid=blockIdx.x;???inti;???if(tid==0)time[bid]=clock();???shared[tid]=0;???for(i=bidTHREAD_NUM+tid;i
利用__shared__声明的变量表示这是sharedmemory，是一个block中每个thread都共享的内存。它会使用在GPU上的内存，所以存取的速度相当快，不需要担心latency的问题。
__syncthreads()是一个CUDA的内部函数，表示block中所有的thread都要同步到这个点，才能继续执行。在我们的例子中，由于之后要把所有thread计算的结果进行加总，所以我们需要确定每个thread都已经把结果写到shared[tid]里面了。
接下来，把main函式的一部份改成：
???intgpudata,result;???clock_ttime;???cudaMalloc((void)&gpudata,sizeof(int)DATA_SIZE);???cudaMalloc((void)&result,sizeof(int)BLOCK_NUM);???cudaMalloc((void)&time,sizeof(clock_t)BLOCK_NUM2);???cudaMemcpy(gpudata,data,sizeof(int)DATA_SIZE,??????cudaMemcpyHostToDevice);???sumOfSquares<<>>(gpudata,result,time);???intsum[BLOCK_NUM];???clock_ttime_used[BLOCK_NUM2];???cudaMemcpy(&sum,result,sizeof(int)BLOCK_NUM,??????cudaMemcpyDeviceToHost);???cudaMemcpy(&time_used,time,sizeof(clock_t)BLOCK_NUM2,??????cudaMemcpyDeviceToHost);???cudaFree(gpudata);???cudaFree(result);???cudaFree(time);???intfinal_sum=0;???for(inti=0;i
可以注意到，现在CPU只需要加总BLOCK_NUM也就是32个数字就可以了。
不过，这个程序由于在GPU上多做了一些动作，所以它的效率会比较差一些。在GeForce8800GT上，它需要约164Kcycles。
当然，效率会变差的一个原因是，在这一版的程序中，最后加总的工作，只由每个block的thread0来进行，但这并不是最有效率的方法。理论上，把256个数字加总的动作，是可以并行化的。最常见的方法，是透过树状的加法：
把kernel改成如下：
__global__staticvoidsumOfSquares(intnum,intresult,???clock_ttime){???extern__shared__intshared[];???constinttid=threadIdx.x;???constintbid=blockIdx.x;???inti;???intoffset=1,mask=1;???if(tid==0)time[bid]=clock();???shared[tid]=0;???for(i=bidTHREAD_NUM+tid;i
后面的while循环就是进行树状加法。main函式则不需要修改。
这一版的程序，在GeForce8800GT上执行需要的时间，大约是140Kcycles（相当于约43GB/s），比完全不在GPU上进行加总的版本还快！这是因为，在完全不在GPU上进行加总的版本，写入到globalmemory的数据数量很大（8192个数字），也对效率会有影响。所以，这一版程序不但在CPU上的运算需求降低，在GPU上也能跑的更快。
进一步改善 上一个版本的树状加法是一般的写法，但是它在GPU上执行的时候，会有sharememory的bankconflict的问题（详情在后面介绍GPU架构时会提到）。采用下面的方法，可以避免这个问题：
???offset=THREAD_NUM/2;???while(offset>0){??????if(tid>=1;??????__syncthreads();???}
这样同时也省去了mask变数。因此，这个版本的执行的效率就可以再提高一些。在GeForce8800GT上，这个版本执行的时间是约137Kcycles。当然，这时差别已经很小了。如果还要再提高效率，可以把树状加法整个展开：
???if(tid<128){shared[tid]+=shared[tid+128];}???__syncthreads();???if(tid<64){shared[tid]+=shared[tid+64];}???__syncthreads();???if(tid<32){shared[tid]+=shared[tid+32];}???__syncthreads();???if(tid<16){shared[tid]+=shared[tid+16];}???__syncthreads();???if(tid<8){shared[tid]+=shared[tid+8];}???__syncthreads();???if(tid<4){shared[tid]+=shared[tid+4];}???__syncthreads();???if(tid<2){shared[tid]+=shared[tid+2];}???__syncthreads();???if(tid<1){shared[tid]+=shared[tid+1];}???__syncthreads();
当然这只适用于THREAD_NUM是256的情形。这样可以再省下约1000cycles左右（约44GB/s）。最后完整的程序文件可以从这里下载。
第二个CUDA程序
前面介绍的计算平方和的程序，似乎没有什么实用价值。所以我们的第二个CUDA程序，要做一个确实有（某些）实用价值的程序，也就是进行矩阵乘法。而且，这次我们会使用浮点数。
虽然矩阵乘法有点老套，不过因为它相当简单，而且也可以用来介绍一些有关CUDA的有趣性质。
矩阵乘法 为了单纯起见，我们这里以方形的矩阵为例子。基本上，假设有两个矩阵A和B，则计算AB=C的方法如下：
???for(i=0;i
一开始，我们先准备好产生数据、设定CUDA等等的工作：
???intmain()???{??????floata,b,c,d;??????intn=1000;??????if(!InitCUDA())return0;??????a=(float)malloc(sizeof(float)nn);??????b=(float)malloc(sizeof(float)nn);??????c=(float)malloc(sizeof(float)nn);??????d=(float)malloc(sizeof(float)nn);??????srand(0);??????matgen(a,n,n);??????matgen(b,n,n);??????clock_ttime=matmultCUDA(a,n,b,n,c,n,n);??????matmult(a,n,b,n,d,n,n);??????compare_mat(c,n,d,n,n);??????doublesec=(double)time/CLOCKS_PER_SEC;??????printf("Timeused:%.2f(%.2lfGFLOPS)\n",sec,?????????2.0nnn/(sec1E9));??????return0;???}
InitCUDA函式和第一个CUDA程序一样，可以直接参考前面的文章。以下是上面用到的一些其它的函式：
产生矩阵：
???voidmatgen(floata,intlda,intn)???{??????inti,j;??????for(i=0;i
这个函式只是利用随机数生成器把矩阵填满0~1之间的数字。特别注意到因为C语言中无法声明变动大小的二维矩阵，所以我们使用ilda+j的方式。
进行矩阵乘法：
???voidmatmult(constfloata,intlda,constfloatb,intldb,??????floatc,intldc,intn)???{??????inti,j,k;??????for(i=0;i
这是以CPU进行矩阵乘法、用来进行验证答案正确与否的程序。特别注意到它用double来储存暂时的计算结果，以提高精确度。
验证结果：
???voidcompare_mat(constfloata,intlda,??????constfloatb,intldb,intn)???{??????floatmax_err=0;?????floataverage_err=0;??????inti,j;??????for(i=0;i
这个函式计算两个矩阵的最大相对误差和平均相对误差，并把结果印出来。
最后是CUDA的矩阵乘法的部份：
???#defineNUM_THREADS256???clock_tmatmultCUDA(constfloata,intlda,??????constfloatb,intldb,floatc,intldc,intn)???{??????floatac,bc,cc;??????clock_tstart,end;??????start=clock();??????cudaMalloc((void)&ac,sizeof(float)nn);??????cudaMalloc((void)&bc,sizeof(float)nn);??????cudaMalloc((void)&cc,sizeof(float)nn);??????cudaMemcpy2D(ac,sizeof(float)n,a,sizeof(float)lda,?????????sizeof(float)n,n,cudaMemcpyHostToDevice);??????cudaMemcpy2D(bc,sizeof(float)n,b,sizeof(float)ldb,?????????sizeof(float)n,n,cudaMemcpyHostToDevice);??????intblocks=(n+NUM_THREADS-1)/NUM_THREADS;??????matMultCUDA<<>>?????????(ac,n,bc,n,cc,n,n);??????cudaMemcpy2D(c,sizeof(float)ldc,cc,sizeof(float)n,??????sizeof(float)n,n,cudaMemcpyDeviceToHost);??????cudaFree(ac);??????cudaFree(bc);??????cudaFree(cc);??????end=clock();??????returnend-start;???}
这个函式相当单纯，就是在显卡内存中配置存放矩阵的内存，然后把主内存中的矩阵数据复制到显卡内存上。不过，因为我们的矩阵乘法函式可以指定pitch（即lda、ldb、和ldc），所以如果用一般的cudaMemcpy函式来复制内存的话，会需要每个row都分开复制，那会需要呼叫很多次cudaMemcpy函式，会使效率变得很差。因此，在这里我们用了一个新的cudaMemcpy2D函式，它是用来复制二维数组，可以指定数组的pitch。这样就可以透过一次函数调用就可以了。
进行计算的kernel如下：
???__global__staticvoidmatMultCUDA(constfloata,size_tlda,??????constfloatb,size_tldb,floatc,size_tldc,intn)???{??????constinttid=threadIdx.x;??????constintbid=blockIdx.x;??????constintidx=bidblockDim.x+tid;??????constintrow=idx/n;??????constintcolumn=idx%n;??????inti;??????if(row
这个函式一开始先从bid和tid计算出这个thread应该计算的row和column，在判断row和column在范围内之后，就直接进行计算，并把结果写到c矩阵中，是非常单纯的函式。
在GeForce8800GT上实际执行的结果如下：
???Maxerror:2.01484e-006Averageerror:3.36637e-007???Timeused:1.1560(1.73GFLOPS)
可以看到两个问题：
很明显的，执行效率相当低落。
最大相对误差偏高。理想上应该要低于1e-6。
计算结果的误差偏高的原因是，在CPU上进行计算时，我们使用double（即64bits浮点数）来累进计算过程，而在GPU上则只能用float（32bits浮点数）。在累加大量数字的时候，由于累加结果很快会变大，因此后面的数字很容易被舍去过多的位数。
由于CUDA的浮点数运算，在进行加、减、乘法时是符合IEEE754规定的精确度的，因此，我们可以利用Kahan''sSummationFormula来提高精确度。把程序改成：
???if(row
修改后的程序的执行结果是：
???Maxerror:1.19209e-007Averageerror:4.22751e-008???Timeused:1.1560(1.73GFLOPS)
可以看到相对误差有很大的改善，效率则没什么变化。
由于Kahan''sSummationFormula需要的运算量提高，但是效率却没有什么改变，可以看出这个kernel主要的瓶颈应该是在内存的存取动作上。这是因为有大量的内存读取是重复的。例如，矩阵a的一个row在每次进行计算时都被重复读入，但这是相当浪费的。这样的计算方式，总共需要读取2n3次内存。如果让一个row只需要读入一次的话，就可以减到为n3+n2次。
?
第一个改良 ?
和我们的第一个CUDA程序一样，我们可以利用sharedmemory来储存每个row的数据。不过，因为只有同一个block的threads可以共享sharedmemory，因此现在一个row只能由同一个block的threads来进行计算。另外我们也需要能存放一整个row的sharedmemory。因此，把先把呼叫kernel的部份改成：
??????matMultCUDA<<>>?????????(ac,n,bc,n,cc,n,n);
kernel的部份则改成：
???__global__staticvoidmatMultCUDA(constfloata,size_tlda,??????constfloatb,size_tldb,floatc,size_tldc,intn)???{??????extern__shared__floatdata[];??????constinttid=threadIdx.x;??????constintrow=blockIdx.x;??????inti,j;??????for(i=tid;i
第一个部份先把整个row读到sharedmemory中，而第二个部份则进行计算，并没有太大的变化。主要的差别是现在一个row只由一个block进行计算。
在GeForce8800GT上，执行的结果是：
???Maxerror:1.19209e-007?Averageerror:4.22751e-008???Timeused:0.4220??(4.74GFLOPS)
很明显的，计算的结果并没有改变，不过速度则提高了超过一倍。虽然如此，但是这样的效率仍不尽理想，因为理论上GeForce8800GT有超过300GFLOPS的运算性能。即使是把Kahan''sSummationFormula所需要的额外运算考虑进去，这样的效率仍然连理论最大值的十分之一都不到。
会有这样的结果，原因其实还是同样的：对内存的存取次数太多了。虽然现在A矩阵的row的数据已经不再需要重复读取，但是B矩阵的column的数据仍然一直被重复读取。
另一个问题比较不是那么明显：对B矩阵的读取，虽然看起来不连续，但实际上它是连续的。这是因为不同的thread会读取不同的column，因此同时间每个thread读取的各个column加起来，就是一个连续的内存区块。那么，为什么效率还是不佳呢？这是因为，GPU上的内存控制器，从某个固定的倍数地址开始读取，才会有最高的效率（例如16bytes的倍数）。由于矩阵大小并不是16的倍数（这里使用的是1000x1000的矩阵），所以造成效率不佳的情形。
要解决这个问题，我们可以在cudaMalloc的时候稍微修改一下，让宽度变成适当的倍数就可以了。但是，适当的倍数是多少呢？幸运的是，我们并不需要知道这些细节。CUDA提供了一个cudaMallocPitch的函式，可以自动以最佳的倍数来配置内存。因此，我们可以把cudaMalloc的部份改成：
???size_tpitch_a,pitch_b,pitch_c;???cudaMallocPitch((void)&ac,&pitch_a,sizeof(float)n,n);???cudaMallocPitch((void)&bc,&pitch_b,sizeof(float)n,n);???cudaMallocPitch((void)&cc,&pitch_c,sizeof(float)n,n);
cudaMallocPitch函式会以适当的倍数配置内存，并把配置的宽度传回。因此，在把矩阵复制到显卡内存上时，要使用它传回的宽度：
???cudaMemcpy2D(ac,pitch_a,a,sizeof(float)lda,??????sizeof(float)n,n,cudaMemcpyHostToDevice);???cudaMemcpy2D(bc,pitch_b,b,sizeof(float)ldb,??????sizeof(float)n,n,cudaMemcpyHostToDevice);
呼叫kernel的部份也需要修改：
???matMultCUDA<<>>??????(ac,pitch_a/sizeof(float),bc,pitch_b/sizeof(float),??????cc,pitch_c/sizeof(float),n);
同样的，把计算结果复制回到主内存时，也要使用传回的宽度值：
???cudaMemcpy2D(c,sizeof(float)ldc,cc,pitch_c,??????sizeof(float)n,n,cudaMemcpyDeviceToHost);
这样就完成了。Kernel部份则不需要修改。
这样的修改有多大的效果呢？在GeForce8800GT上执行，结果如下：
???Maxerror:1.19209e-007?Averageerror:4.22751e-008???Timeused:0.1250??(16.00GFLOPS)
可以看到，执行速度又再大幅提高了三倍多！而这只是把内存的配置方式稍微修改一下而已。
虽然执行速度提高了很多，但是，和前面提到的理论值相比，其实还是有相当的差距。这是因为，前面也提到过，这样的做法需要n3+n2次的内存读取，和n2次的内存写入动作。由于n=1000，每个数字的大小是32bits，所以总共的内存存取数据量约为4GB。除以实际执行的时间0.125秒，得到的带宽数值是约32GB/s，这已经接近GeForce8800GT显卡内存的带宽了。由于我们计算时间的时候，把配置内存、以及数据的复制动作也计算进去，因此实际上花费在kernel的时间是更短的（约0.09秒）。因此，可以很明显的看出，这个程序的效率，是受限于内存带宽的。
进一步的改良 上一节的结论显示出，矩阵乘法的程序，效率是受限于内存带宽的。那有没有办法降低内存的存取次数呢？答案当然是有的，不然就不会有这一节了:)
要进一步降低内存带宽的使用，可以注意到，在上一节的方法中，虽然A矩阵的存取次数被减至最低，但是B矩阵的存取次数并没有减少。这是因为我们只将A矩阵的row加载到sharedmemory中，但是B矩阵的column也是有被重复使用的。理想上应该也可以避免重复加载才对。不过，由于B矩阵的column使用的时机，和A矩阵的row是不同的，所以并不能直接这样做。
解决方法是"blocking"。也就是把整个矩阵乘法的动作，切割成很多小矩阵的乘法。例如，要计算C矩阵的(0,0)~(15,15)的值，可以把它想成：
???A(0~15,0~15)B(0~15,0~15)+A(0~15,16~31)B(16~31,0~15)???+A(0~15,32~47)B(32~47,0~15)+...
这样一来，我们就可以把两个小矩阵加载到sharedmemory，则小矩阵本身的乘法就不需要再存取任何外部的内存了！这样一来，假设小矩阵的大小是k，则实际上需要的内存存取次数就会变成约2k2(n/k)3=2n3/k。
由于目前CUDA每个block的thread数目最多是512，因此k=16似乎是一个相当理想的数字（共256个threads）。因此，对于一个n=1000的矩阵来说，我们可以把内存存取的量减少到约500MB，也就是上一节的存取量的1/8。理论上，这样应该可以让效率提高八倍（假设没有遇到别的瓶颈）。
为了方便进行区块的计算，我们让每个block有16x16个threads，再建立(n/16)x(n/16)个blocks。把呼叫kernel的地方改成：
???intbx=(n+BLOCK_SIZE-1)/BLOCK_SIZE;???dim3blocks(bx,bx);???dim3threads(BLOCK_SIZE,BLOCK_SIZE);???matMultCUDA<<>>(ac,pitch_a/sizeof(float),??????bc,pitch_b/sizeof(float),cc,pitch_c/sizeof(float),n);
BLOCK_SIZE则是定义成16。dim3是CUDA的一种数据型态，表示一个3D的向量。在这里，我们透过dim3来建立16x16个threads的block，和(n/16)x(n/16)个blocks。
Kernel程序的部份，则改成：
???__global__staticvoidmatMultCUDA(constfloata,size_tlda,??????constfloatb,size_tldb,floatc,size_tldc,intn)???{??????__shared__floatmatA[BLOCK_SIZE][BLOCK_SIZE];??????__shared__floatmatB[BLOCK_SIZE][BLOCK_SIZE];??????constinttidc=threadIdx.x;??????constinttidr=threadIdx.y;??????constintbidc=blockIdx.xBLOCK_SIZE;??????constintbidr=blockIdx.yBLOCK_SIZE;??????inti,j;??????floatresults=0;??????floatcomp=0;??????for(j=0;j
注意到因为我们现在使用16x16的threads，因此threadIdx变量可以取得threadIdx.x和threadIdx.y，范围分别是0~15。blockIdx.x和blockIdx.y变量也是同样的情形，范围分别是0~n/16。
在程序中，因为矩阵的大小不一定会是16的倍数，因此需要使用if判断式检查是否超出矩阵范围。
这个版本在GeForce8800GT上的执行结果如下：
???Maxerror:1.19209e-007?Averageerror:4.22751e-008???Timeused:0.0780??(25.64GFLOPS)
速度虽然提高了，但是似乎并没有达到预期中的八倍。当然，前面提到过，我们在计算时间时，把一些复制内存、配置内存的动作也计算在内，这些动作的时间并不会缩短。实际上kernel的运行时间，大约是0.053秒左右（约略相当于38GFLOPS），比上一节的版本快了将近一倍。
如果这一版程序已经不再限于内存带宽，那为什么没有达到预期的效率呢？这是因为这一版程序已经是限于指令周期了。除了使用Kahan''sSummationFormula会需要更多的运算之外，程序中也有大量计算矩阵地址的乘法等等，这都会需要花费运算资源。另外，那些用来判断超出矩阵范围的if判断式，也会有一定的影响。
要把那些if判断式去掉，有一个方法是，在配置内存时，就配置成16的倍数，并在复制矩阵到显卡内存之前，先将它清为0。如下所示：
???intnewn=((n+BLOCK_SIZE-1)/BLOCK_SIZE)BLOCK_SIZE;???cudaMallocPitch((void)&ac,&pitch_a,??????sizeof(float)newn,newn);???cudaMallocPitch((void)&bc,&pitch_b,??????sizeof(float)newn,newn);???cudaMallocPitch((void)&cc,&pitch_c,??????sizeof(float)newn,newn);???cudaMemset(ac,0,pitch_anewn);???cudaMemset(bc,0,pitch_bnewn);
?这样一来，我们就可以把kernel中的if判断式都移除了：
???__global__staticvoidmatMultCUDA(constfloata,size_tlda,??????constfloatb,size_tldb,floatc,size_tldc,intn)???{??????__shared__floatmatA[BLOCK_SIZE][BLOCK_SIZE];??????__shared__floatmatB[BLOCK_SIZE][BLOCK_SIZE];??????constinttidc=threadIdx.x;??????constinttidr=threadIdx.y;??????constintbidc=blockIdx.xBLOCK_SIZE;??????constintbidr=blockIdx.yBLOCK_SIZE;??????inti,j;??????floatresults=0;??????floatcomp=0;??????for(j=0;j
这个版本的执行结果是：
???Maxerror:1.19209e-007?Averageerror:4.22751e-008???Timeused:0.0780??(25.64GFLOPS)
似乎没有改善。不过，实际上kernel的运行时间已经减少到0.042秒（约略相当于48GFLOPS）。
结论 有些读者可能会想，如果把block再变得更大（例如32x32）是否会有帮助呢？当然，由于最后的程序已经不再是受限于内存带宽（在0.042秒内存取500MB的数据约相当于12GB/s的带宽），所以把block再加大并不会有帮助了。而且，由于一个block内的thread数目最多只能到512个，将block变大也会造成很多额外负担。而且sharedmemory的大小也有限制（GeForce8800GT的sharedmemory大小限制是16384bytes），所以也不能任意增加block的大小。
最后一版程序的完整档案可以从这里下载。
GPU的硬件架构
这里我们会简单介绍，NVIDIA目前支持CUDA的GPU，其在执行CUDA程序的部份（基本上就是其shader单元）的架构。这里的数据是综合NVIDIA所公布的信息，以及NVIDIA在各个研讨会、学校课程等所提供的数据，因此有可能会有不正确的地方。主要的数据源包括NVIDIA的CUDAProgrammingGuide1.1、NVIDIA在Supercomputing''07介绍CUDA的session，以及UIUC的CUDA课程。
GPU的基本介绍 目前NVIDIA推出的显示芯片，支持CUDA的是G80系列的显示芯片。其中G80显示芯片支持CUDA1.0版，而G84、G86、G92、G94、G96则支援CUDA1.1版。基本上，除了最早的GeForce8800Ultra/GTX及320MB/640MB版本的GeForce8800GTS、Tesla等显卡是CUDA1.0版之外，其它GeForce8系列及9系列显卡都支持CUDA1.1。详细情形可以参考CUDAProgrammingGuide1.1的AppendixA。
所有目前支持CUDA的NVIDIA显示芯片，其shader部份都是由多个multiprocessors组成。每个multiprocessor里包含了八个streamprocessors，其组成是四个四个一组，也就是说实际上可以看成是有两组4D的SIMD处理器。此外，每个multiprocessor还具有8192个寄存器，16KB的sharememory，以及texturecache和constantcache。大致上如下图所示：
详细的multiprocessor信息，都可以透过CUDA的cudaGetDeviceProperties()函式或cuDeviceGetProperties()函式取得。不过，目前还没有办法直接取得一个显示芯片中有多少multiprocessor的信息。
在CUDA中，大部份基本的运算动作，都可以由streamprocessor进行。每个streamprocessor都包含一个FMA（fused-multiply-add）单元，可以进行一个乘法和一个加法。比较复杂的运算则会需要比较长的时间。
执行过程 在执行CUDA程序的时候，每个streamprocessor就是对应一个thread。每个multiprocessor则对应一个block。从之前的文章中，可以注意到一个block经常有很多个thread（例如256个），远超过一个multiprocessor所有的streamprocessor数目。这又是怎么回事呢？
实际上，虽然一个multiprocessor只有八个streamprocessor，但是由于streamprocessor进行各种运算都有latency，更不用提内存存取的latency，因此CUDA在执行程序的时候，是以warp为单位。目前的CUDA装置，一个warp里面有32个threads，分成两组16threads的half-warp。由于streamprocessor的运算至少有4cycles的latency，因此对一个4D的streamprocessors来说，一次至少执行16个threads（即half-warp）才能有效隐藏各种运算的latency。
由于multiprocessor中并没有太多别的内存，因此每个thread的状态都是直接保存在multiprocessor的寄存器中。所以，如果一个multiprocessor同时有愈多的thread要执行，就会需要愈多的寄存器空间。例如，假设一个block里面有256个threads，每个thread用到20个寄存器，那么总共就需要256x20=5,120个寄存器才能保存每个thread的状态。
目前CUDA装置中每个multiprocessor有8,192个寄存器，因此，如果每个thread使用到16个寄存器，那就表示一个multiprocessor同时最多只能维持512个thread的执行。如果同时进行的thread数目超过这个数字，那么就会需要把一部份的数据储存在显卡内存中，就会降低执行的效率了。
编者注：在NVIDIAGT200中的RegisterFile大小增加了一倍，在FP32下可用的registerfile为16K，FP64下是8K。
Sharedmemory 目前CUDA装置中，每个multiprocessor有16KB的sharedmemory。Sharedmemory分成16个bank。如果同时每个thread是存取不同的bank，就不会产生任何问题，存取sharedmemory的速度和存取寄存器相同。不过，如果同时有两个（或更多个）threads存取同一个bank的数据，就会发生bankconflict，这些threads就必须照顺序去存取，而无法同时存取sharedmemory了。
Sharedmemory是以4bytes为单位分成banks。因此，假设以下的数据：
???__shared__intdata[128];
那么，data[0]是bank0、data[1]是bank1、data[2]是bank2、…、data[15]是bank15，而data[16]又回到bank0。由于warp在执行时是以half-warp的方式执行，因此分属于不同的halfwarp的threads，不会造成bankconflict。
因此，如果程序在存取sharedmemory的时候，使用以下的方式：
???intnumber=data[base+tid];
那就不会有任何bankconflict，可以达到最高的效率。但是，如果是以下的方式：
???intnumber=data[base+4tid];
那么，thread0和thread4就会存取到同一个bank，thread1和thread5也是同样，这样就会造成bankconflict。在这个例子中，一个halfwarp的16个threads会有四个threads存取同一个bank，因此存取sharememory的速度会变成原来的1/4。
一个重要的例外是，当多个thread存取到同一个sharedmemory的地址时，sharedmemory可以将这个地址的32bits数据「广播」到所有读取的threads，因此不会造成bankconflict。例如：
???intnumber=data[3];
这样不会造成bankconflict，因为所有的thread都读取同一个地址的数据。
很多时候sharedmemory的bankconflict可以透过修改数据存放的方式来解决。例如，以下的程序：
???data[tid]=global_data[tid];???...???intnumber=data[16tid];
会造成严重的bankconflict，为了避免这个问题，可以把数据的排列方式稍加修改，把存取方式改成：
???introw=tid/16;???intcolumn=tid%16;???data[row17+column]=global_data[tid];???...???intnumber=data[17tid];
这样就不会造成bankconflict了。
编者注：sharememory在NVIDIA的文档中其实还有不同的叫法，例如PDC（ParallelDataCache）、PBSM（per-blocksharememory）。
Globalmemory 由于multiprocessor并没有对globalmemory做cache（如果每个multiprocessor都有自己的globalmemorycache，将会需要cachecoherenceprotocol，会大幅增加cache的复杂度），所以globalmemory存取的latency非常的长。除此之外，前面的文章中也提到过globalmemory的存取，要尽可能的连续。这是因为DRAM存取的特性所造成的结果。
更精确的说，globalmemory的存取，需要是"coalesced"。所谓的coalesced，是表示除了连续之外，而且它开始的地址，必须是每个thread所存取的大小的16倍。例如，如果每个thread都读取32bits的数据，那么第一个thread读取的地址，必须是164=64bytes的倍数。
如果有一部份的thread没有读取内存，并不会影响到其它的thread速行coalesced的存取。例如：
???if(tid!=3){??????intnumber=data[tid];???}
虽然thread3并没有读取数据，但是由于其它的thread仍符合coalesced的条件（假设data的地址是64bytes的倍数），这样的内存读取仍会符合coalesced的条件。
在目前的CUDA1.1装置中，每个thread一次读取的内存数据量，可以是32bits、64bits、或128bits。不过，32bits的效率是最好的。64bits的效率会稍差，而一次读取128bits的效率则比一次读取32bits要显著来得低（但仍比non-coalesced的存取要好）。
如果每个thread一次存取的数据并不是32bits、64bits、或128bits，那就无法符合coalesced的条件。例如，以下的程序：
???structvec3d{floatx,y,z;};???...???__global__voidfunc(structvec3ddata,floatoutput)???{??????output[tid]=data[tid].xdata[tid].x+?????????data[tid].ydata[tid].y+?????????data[tid].zdata[tid].z;???}
并不是coalesced的读取，因为vec3d的大小是12bytes，而非4bytes、8bytes、或16bytes。要解决这个问题，可以使用__align(n)__的指示，例如：
???struct__align__(16)vec3d{floatx,y,z;};
这会让compiler在vec3d后面加上一个空的4bytes，以补齐16bytes。另一个方法，是把数据结构转换成三个连续的数组，例如：
???__global__voidfunc(floatx,floaty,floatz,floatoutput)???{??????output[tid]=x[tid]x[tid]+y[tid]y[tid]+?????????z[tid]z[tid];???}
如果因为其它原因使数据结构无法这样调整，也可以考虑利用sharedmemory在GPU上做结构的调整。例如：
???__global__voidfunc(structvec3ddata,floatoutput)???{??????__shared__floattemp[THREAD_NUM3];??????constfloatfdata=(float)data;??????temp[tid]=fdata[tid];??????temp[tid+THREAD_NUM]=fdata[tid+THREAD_NUM];??????temp[tid+THREAD_NUM2]=fdata[tid+THREAD_NUM2];??????__syncthreads();??????output[tid]=temp[tid3]temp[tid3]+?????????temp[tid3+1]temp[tid3+1]+?????????temp[tid3+2]temp[tid3+2];???}
在上面的例子中，我们先用连续的方式，把数据从globalmemory读到sharedmemory。由于sharedmemory不需要担心存取顺序（但要注意bankconflict问题，参照前一节），所以可以避开non-coalesced读取的问题。
Texture CUDA支援texture。在CUDA的kernel程序中，可以利用显示芯片的texture单元，读取texture的数据。使用texture和globalmemory最大的差别在于texture只能读取，不能写入，而且显示芯片上有一定大小的texturecache。因此，读取texture的时候，不需要符合coalesced的规则，也可以达到不错的效率。此外，读取texture时，也可以利用显示芯片中的texturefiltering功能（例如bilinearfiltering），也可以快速转换数据型态，例如可以直接将32bitsRGBA的数据转换成四个32bits浮点数。
显示芯片上的texturecache是针对一般绘图应用所设计，因此它仍最适合有区块性质的存取动作，而非随机的存取。因此，同一个warp中的各个thread最好是读取地址相近的数据，才能达到最高的效率。
对于已经能符合coalesced规则的数据，使用globalmemory通常会比使用texture要来得快。
运算单元 Streamprocessor里的运算单元，基本上是一个浮点数的fusedmultiply-add单元，也就是说它可以进行一次乘法和一次加法，如下所示：
???a=bc+d;
compiler会自动把适当的加法和乘法运算，结合成一个fmad指令。
除了浮点数的加法及乘法之外，整数的加法、位运算、比较、取最小值、取最大值、及以型态的转换（浮点数转整数或整数转浮点数）都是可以全速进行的。整数的乘法则无法全速进行，但24bits的乘法则可以。在CUDA中可以利用内建的__mul24和__umul24函式来进行24bits的整数乘法。
浮点数的除法是利用先取倒数，再相乘的方式计算，因此精确度并不能达到IEEE754的规范（最大误差为2ulp）。内建的__fdividef(x,y)提供更快速的除法，和一般的除法有相同的精确度，但是在2216
此外CUDA还提供了一些精确度较低的内部函数，包括__expf、__logf、__sinf、__cosf、__powf等等。这些函式的速度较快，但精确度不如标准的函式。详细的数据可以参考CUDAProgrammingGuide1.1的AppendixB。
和主内存间的数据传输 在CUDA中，GPU不能直接存取主内存，只能存取显卡上的显示内存。因此，会需要将数据从主内存先复制到显卡内存中，进行运算后，再将结果从显卡内存中复制到主内存中。这些复制的动作会限于PCIExpress的速度。使用PCIExpressx16时，PCIExpress1.0可以提供双向各4GB/s的带宽，而PCIExpress2.0则可提供8GB/s的带宽。当然这都是理论值。
从一般的内存复制数据到显卡内存的时候，由于一般的内存可能随时会被操作系统搬动，因此CUDA会先将数据复制到一块内部的内存中，才能利用DMA将数据复制到显卡内存中。如果想要避免这个重复的复制动作，可以使用cudaMallocHost函式，在主内存中取得一块pagelocked的内存。不过，如果要求太大量的pagelocked的内存，将会影响到操作系统对内存的管理，可能会减低系统的效率。
|
|
