CUDA2.4-原理之性能优化及浮点运算

本部分来自于《大规模并行处理器编程实战》第六章、第七章。打算不再看这本书了，准备看《programming massively parallel processors 2nd》，即它的第二版，第一版是09年的，第二版是13年的，虽说第二版可是里面涉及的是cuda4.0 和5.0，然而现在2015年7月，cuda都7.0了，正所谓赶速度，完全赶不上啊。虽然说本书好，不过一个不小心，你费老大劲做的优化，发现其实新版本的cuda或者硬件完全不需要，果然有关cuda的最好的资料其实还是官方文档，因为这些完全赶不上速度啊。

一、性能优化

        CUDA kernel函数的执行速度很大程度上取决于每个设备的资源约束。而且不同的应用程序中，不同的约束可能决定并成为限制因素。在特定的CUDA设备上，可以通过一种资源代替另一种资源来提高英程序的性能。合理的策略有可能提高性能，也可能不起作用，所以需要测试，本章说的这些可以用来培养程序员对算法的直觉，如何来提高整体的性能。

1.1 线程执行问题

        相比较前面几个博文，没有过多的讨论每个块中线程的执行时间问题。在本书发售的时候，Nvidia公司是通过对块中的线程进行捆绑执行的，即执行的不是单个线程而是一个warp（包含32个线程），这样做降低硬件成本，而且一定程度上优化了存储器访问的服务。对于划分warp来说，如果是一维的，那么对于最后一个warp的划分如果不满32个，会将其他块中的线程拉过来补全成32个再执行。如果对于多维线程的块来说，划分warp前会把维度映射成一个线性顺序，y 和z 的坐标小的放前面，大的放后面。假如一个块有二维的线程，那么将所有threadIdx.y是1的线程放在threadIdx.y是0的线程后面（注意这里是y 不是x，即按照 Fortran语言的顺序，前面的变化快，后面的变化慢），以此类推（这里要注意这种访问形式，在后面的1.2很有用，相当于二维的矩阵是转置的，不是正常的那种形式，注意，不过如果将y想象成矩阵的行，x想像成矩阵的列，那就没问题了）。


图1

上图就是将二维的线程块映射成一维的顺序形式。（注意顺序）。而对于三维的块来说，首先将threadIdx.z为0的所有线程按照线性顺序排列，以此类推。将warp中所有线程执行玩一个指令，再去执行另一个指令，这种执行方式叫做SIMT（single Instruction multiple thread，单指令多线程）。如果在一个warp中所有的线程执行的是相同的指令，那么工作效率最高，如果说在if-then-else结构中，决策条件基于threadIdx的值，比如if（threadIdx>2）｛｝这种的，那么会导致线程分支，即按照两个控制流分支路径，这样线程0、1、2和3、4、5就不同；如果循环条件基于threadIdx的值，则循环也会引起线程分支。这种用法自然会导致很多重要的并行算法的产生。

        这里举个归约数组和的算法：即对于一个有着N个元素的数组来说，想要求的所有元素的和，那么归约采取的是两两结合，即假如10个元素（(1,2),(3,4),(5,6),(7,8),(9,10)），然后第二次（（1,2,3,4）,(5,6,7,8),(9,10)）。这种类型的，两两归约，因为前面计算一次之后就只有10/2=5个数值，所以第二次只需要计算5/2=2次（每次执行完一次，再重新分配，比如前面的第二次（9,10）因为缺少2个，就无法被4整除，所以第二次只计算了2次，然后（1,2,3,4,5,6,7,8）,(9,10)）接着最后得到最终答案；代码如下：

1.__shared__float partialSum[]

2.unsigned int t = threadIdx.x;
3.for(unsigned int stride = 1;
4.    stride < blockDim.x; stride *= 2)
5.{
6.__syncthreads();
7.if(t %( 2 * stride) == 0)
8.   partialSUm[t] += paritialSum[t + stride];
9.}
执行过程如下：

图2

注意：上述的归约其实导致块中有一般的线程从来不执行，是很浪费的，所以需要修改kernel函数，留待作业。

在上面的代码中，就是当threadIdx.x的值为偶数的时候才执行加法，所以会导致那些不执行第8行代码的线程需要通过另外一条路径。下面是修改的kernel函数，即不再采用两两相加的方法，这样的是具有较少线程分支的kernel函数：

1.__shared__float partialSum[]

2.unsigned int t = threadIdx.x;
3.for(unsigned int stride = blockDim.x>>1;
4.    stride > 0; stride >>= 1)
5.{
6.__syncthreads();
7.if(t < stride)
8.   partialSUm[t] += paritialSum[t + stride];
9.}
修改之后的代码的性能有所不同，上述代码中通过移位来代替除法操作，降低开销。执行过程如下：

图3

因为第一次迭代中，线程0-255都执行加法，而线程256-511不执行加法。由于warp中包含的32个线程对应的threadIdx.x值是连续的，因此都0-7个warp所有线程都执行加法，第8-15个warp则跳过加法，由于warp中所有线程都通过相同的路径，所以没有线程分支。可是还是因为有if的存在，kernel函数的分支并未完全消除。在执行第5次迭代的时候，第8行代码的线程个数低于32.也就是，最后5次迭代中分别只有16、8、4、2、1个线程执行加法运算，所以仍然存在分支。

1.2全局存储器的带宽
        制约CUDA kernel函数的一个重要因素就是全局存储器的访问数据。之前有讨论过如何减少访问的流量来达到加速的目的，这里接着讨论存储器合并技术。使得更加有效的将数据从全局存储器中移动到共享存储器和寄存器上。因为cuda系统采用的是DRAM的全局存储器，这种DRAM单元为了加快数据访问的速度，采用并行进程的方式，即当DRAM芯片中的传感器接收到请求的单元的索引的时候，会顺带把其附近的单元的电位一起传送过来，如果应用程序在访问单元改变前能够充分利用这种来自多个连续单元的数据，则会比真正的随机顺序的单元访问要块得多。所以需要kernel函数安排数据的访问顺序。（本书发布时，现在不知道cuda的设备还是不是采用DRAM）在G80/GT200中，考虑这样一个事实：同一个waro中的线程在任何给定的时间内都执行同一条指令，也就是说当同一个wari中所有线程执行同一条指令访问全局存储器中的连续单元时，这种访问模式是最好的。如线程0访问单元N（N应该是与16个字的边界对齐，也就是说，N的低16位应该为0），线程1访问单元N+1,以此类推。那么所有的这些访问请求将会被合并和结合成对所有连续单元的单个请求，使得DRAM在传输数据的速度接近全局存储器带宽的峰值。


图4

上图为合并前和后的访问图。假设合并前，同一个warp中所有线程第0次迭代时访问0-31行中第0列的元素，第1次迭代访问0-31行中第1列的元素，这样没法合并；可以第0次迭代的时候第0个warp中的线程读取第0-31列的第0行元素。为何这样更有利，这需要了解矩阵元素在全局存储器中是如何存放的。

       在全局存储器中其实和主机内存一样，地址是线性的，也就是如果全局存储器包含1024个单元，那么访问地址其实是0-1023.对于GT200来说，其地址范围为0-2^32-1。在C和CUDA中矩阵元素都是以行为主的约定存放的，即如下图形式：


图5

M(0,0)和M(0,1)之间隔了4个元素


图6

上图中，在第1次迭代中，同一个warp中线程由循环发起访问矩阵中所有列的第0个元素。因为所有列的第0个元素存放在全局存储器中连续的单元中。这样硬件检测到，就会合并成一个综合的访问，这样DRAM就支持高速的数据访问；下面就是合并前的访问形式：


图7

下面是个矩阵乘法的技术，如果说算法本身要求kernel函数以行的方式遍历数据，那么可以使用共享存储器来实现存储空间合并，先通过线程协作方式，即以合并的方式加载数据到共享存储器，然后在共享存储器上不管是以行还是列都不会引起性能的太大差异，因为共享存储器本身能实现告诉的片上存储器，无需通过合并来提高数据的访问速度。


图8

__global__ void MatrixMulKernel(float* Md,float* Nd,float* Pd, int Width)
{
1.__shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
2.__shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

3.int bx = blockIdx.x;  int by = blockIdx.y;
4.int tx = threadIdx.x; int ty = threadIdx.y;

//标识要处理的Pd元素对应的行和列
5.int Row = by * TILE_WIDTH + ty;
6.int Col = bx * TILE_WIDTH + tx;

7.float Pvalue = 0;
//循环Md和Nd中要计算Pd元素的块
8.for(int m = 0; m < Width/TILE_WITDH; ++m){

//相互协作把Md和Nd的块加载到共享存储器中
9.    Mds[ty][tx] = Md[Row * Width + (m * TILE_WITDH + tx)];
10.   Nds[ty][tx] = Nd[(m * TILE_WIDTH + ty) * Width + Col];
11.   __syncthreads();

12.for(int k = 0;k <TILE_WITDH; ++k)
13.   Pvalue += Mds[ty][k] * Nds[k][tx];

14.Pd[Row][Col] = Pvalue;
  }
}
（个人：注意ty，tx的顺序，因为tx变化快，而且刚好矩阵在存储器中是行存储的，刚好对上）上述代码是使用共享存储器的分块矩阵乘法的kernel。在第10行Nd的索引的所有项对于warp中所有线程都有相同的值（最先变化的是tx，在本博文最上面说了），除了Col中的tx项，即同一个warp中所有线程对应的tx值是相邻的；所以一个warp中相邻的线程访问一行中相邻的元素，会将它们合并（个人：方法就是先固定ty，假设只有tx是变量，然后观察访问的全局存储器的地址是否相邻即可）（个人：这里的Mds和Nds都是按照块读取的，即一次是读取一个矩阵，不是单独的Md行和Nd列，参考上面那个图的形式）。

1.3 SM资源的动态划分
    SM中的执行资源有寄存器、线程块槽和线程槽。比如在GT200中每个SM支持1024个线程槽，每个槽容纳一个线程，运行时划分这些线程槽给线程块，如果每个线程块有256个线程，就将1024个线程槽划分给4个块，如果每个线程块有128个线程，就划分1024个线程槽给8个线程块。这种是动态划分，不是固定划分。所以有好处也有坏处，坏处是会导致资源限制之间有轻微的相互影响，导致资源利用率降低，因为SM中块槽支持8个，所以最多容纳8个块，算下来每个块至少要有128个线程才能充分利用。

        寄存器文件是另一种动态划分的资源。在每个CUDA设备中，寄存器的数量不是在编程模型中指定的，而是在实现中可变的。比如G80的每个SM支持8192个寄存器文件。比如kernel中每个线程使用10个寄存器，那么块大小为16*16，那么一个块用到16*16*10=2560个。那么每个SM中按照寄存器的限制只能放入3个块，所以会导致总的线程为768个，虽然这时候块槽和线程槽都不是限制因素。假如每个线程使用11个寄存器，那么总的三个快需要8448个寄存器，超过了寄存器限制，按照前面说的，SM通过减少块的形式来处理（下图b），所以使得寄存器数目减少16*16*11=2816，SM中只有两个块，线程也变成了512个进驻，也就是说，多使用一个自动变量（存放入寄存器），在G80上执行的程序中warp的并行性减少1/3。参考资源《cuda occupany calculatar，nvidia 2009》，这是一个可下载的excel表格，在给定kernel函数使用资源的情况下，对某个特殊的设备实现可以计算出在每个SM上运行的线程的实际数量。


图9.资源限制的相互影响

         不过在有的情况下可以通过增加自动变量提高单个线程的执行速度，因为将其放入寄存器中，就无需去访问耗时的存储器了，比如在kernel函数中，全局存储器的加载和使用有4个独立的指令，在G80中，处理每条指令需要4个时钟周期，所以这4条指令需要16个时钟周期，如果全局存储器有200个时钟周期的延时，为了充分利用执行单元，至少需要200/(4*4)=14个warp来保证零开销调度。（即增加寄存器是为了加快执行时间）。假设将独立指令数目从4增加到8，那么间歇性访问存储器需要32个时钟周期，如果全局存储器中还是200个时钟周期的延时，那么为了充分利用执行单元需要200/(4*8)=7个warp保证零开销调度，也就是上面的即使块数目从3减到2，此时warp的数据也从24减少到16，可是每个SM中仍有足够的warp在执行，反而性能其实还是提高的。这种参数分配需要不断的调试，去测试，为了减少这种人力的趺粗。Ryoo等人提出了一种用于实现实验过程自动化的方法（参考文献见最下面）。

1.4数据预取

         通过使用之前介绍的分块数据使用共享存储器访问解决全局存储器带宽的限制问题，使得在一些warp在等待的时候，或者说只需要先计算一半，等剩下数据到了在计算另一半结果。而对于那些所有线程都等待其存储器访问结果的时候，这种方法也有不足了，即当所有线程在存储器访问指令和已访问的数据使用指令之间的独立指令很少，就会这样。

       一种方法解决上述问题是：当使用当前数据元素的时候预取下一个数据元素，通过预取技术和分块技术的结合来解决带宽限制和长延时的问题。如下图：


图10

         上图A对应图8的分块矩阵乘法的kernel函数，对应的图8下面的代码的第9，10行就是将当前块加载到共享存储器，不过第9行其实是执行了两部分：将Md的元素从全局存储器的单元加载到寄存器；将寄存器的内容存储到共享存储器单元中。这两部分之间没有独立的指令，所以那些加载当前块的warp可能需要等待很长时间，之后才能计算当前块，所以，可能需要很多的warp才能让浮点单元保持忙碌。图10中B就是修改的预取版本。在进入循环之前，将第一个块加载到寄存器，进入循环后，将已加载的数据移动到共享存储器中，对应图8下面代码的第12和13行的点积循环，（个人，类似于cpp常写的索引自加操作）当循环迭代时，这个当前迭代中“下一个块”变成下一次迭代的“当前块”。在执行点积循环在两者之间增加很多独立的指令，减少了线程必须等待全局存储器访问数据的时间。虽然这样用到额外的寄存器会减少SM上的块，不过如果显著减少了每个线程等待其全局存储器加载数据的时间，这技术还是有优势的。

1.5指令混合

         在当前（本书的时间）的CUDA GPU 中，每个处理器的内核都只有有限的指令处理带宽。不管是浮点计算还是加载指令还是分支指令，都需要占用指令处理带宽。

for(int k = 0;k < BLOCK_SIZE; ++k)
    Pvalue += Ms[ty][k] * Ns[k][tx];
   (a)循环引入指令开销

Pvalue += Ms[ty][0] * Ns[0][tx]+ ...
            Ms[ty][15] * Ns[15][tx];
   (b)循环展开可以提高指令的混合水平
上面的代码中，对于（a）来说，循环引入的更新循环计数变量k 和每次迭代后的分支指令，而且在矩阵Ms和Ns的索引中用到循环技术变量k,会引入地址运算指令。相比较有2条浮点运算指令、循环分支指令、2条地址运算指令，1条循环增量计数指令。也就是说只有1/3的指令是执行浮点运算的 。所以这样的性能限制在带宽峰值的1/3以内。相比较下面的展开可以消除分支指令、循环计数和地址运算指令（常量偏移即可）。在未来（基于本书的时间），编译器将会支持自动展开，不过在工具未成熟之前，还是需要在源码中展开的。

1.6 线程粒度
       在性能优化中一个重要的算法决定因素是线程的粒度。应该尽可能的把更多工作放在每个线程中并采取更少的线程是有优势的额，特别是线程之间有重复的工作时。


图11 用矩形块增大线程的粒度

上图就是之前那个矩阵乘法的例子。线程粒度优化的可能性在于多个快重复的加载Md的每个块。在之前的算法中两个块生成Pd的两个块需要重复加载Md的同一行。如果将两个线程块合并成一个，就可以消除这种冗余。在新的线程块中，每个线程现在负责计算Pd的两个元素。修改kernel使得函数一次可以计算两个点积。两个点积都使用相同的行Mds和不同的列Nds进行运算。这样将全局存储器的访问次数减少1/4。不过这样做的负面效果就是新的kernel函数现在使用更多的寄存器和共享存储器。因此，在SM上能运行的块的数目也可能相应减少，线程块的总数会减少一半，这对于维度较小的矩阵，会导致并行性不足。对于G80/G200中一个2048*2048的矩阵乘法，合并水平方向上4个相邻的块来计算水平方向上相邻的块，性能最好。

1.7可度量的性能


图12 不同技术的性能度量

上述用到的技术有分块、循环展开、数据预取和线程粒度，左边的是8*8的块，右边是16*16的块，每个块的大小，比较了基础代码、粒度调整中合并两个块后的代码和合并4个块后的代码。结论是： 

       1）块的大小在性能中起到主要作用。在块大小未达到16*16时，循环展开和数据预取都不起作用。因为对于8*8的块来说，全局存储器的带宽处于饱和状态。并且1*2的矩形分块将全局存储器的访问次数减少了1/4，而1*4的减少了3/8，不过在1*8的矩形块只能减少7/16，因为使用更大的块能提高的性能呈递减趋势，这使得加大块的大小失去了吸引力。

      2）当块的大小足够大时（比如这里的16*16）,那么减轻全局存储器带宽的饱和状态，循环展开和数据预取就变得重要了。

     3）虽然数据预取对于1*1的块有用，但是在1*2的矩形分块中基本不起作用。对于1*4的矩形分块，数据预期中一个16*16的块需要使用的寄存器已经超过SM中寄存器的总数，这使得代码在G80上没法执行。不过通过多种技术调整可以看出在G80上执行速度可以从18 GFLOPS增加到120 GFLOPS。


二、浮点运算

2.1浮点格式

        浮点标准有两个IEEE-754（1985年）和IEEE 754-2008【IEEE，2008】，几乎绝大部分计算机厂商都接受的。理解这些概念对于应用程序的开发还是很重要的。最开始，浮点数系统采用位模式来表示数值。在IEEE浮点标准中，一个数值的表示由3部分组成：符号位（S）、阶码（E）和尾数（M）。每个（S，E，M）模式根据下面的格式可以标识一个唯一的数值：


S代表的意义：当S=0时表示是一个正数，当S=1时表示一个负数。

这里为了举例容易点，假设每个浮点数包含1位符号位，3位阶码和2位尾数。例如，0.5D表示5*10^(-1)，0.1B表示1*2*(-1)。

2.1.1M的规范化表示

在上面的公式中，要求M的范围为。这使得对每个浮点数来说，它的尾数是唯一的。0.5D的唯一尾数就是M=1.0：

0.5D=1.0B*2(-1)

另一种表示为0.1B*2^0，不过根据表示规则，这种表示法的尾数值太小，10.0B*2^(-2)的尾数值太大。在这种限定下，每个浮点数只有一个符合规格的尾数值。满足这种限制条件的数称之为规格化数。因为所有的尾数值都满足1.xx这种形式的限制条件，所以可以省略前面的"1."，比如0.5对应的尾数在两位尾数表示法中为00，这是1.00中省略1.得到的。也就是在IEEE中，n位的尾数其实可以有效的表示（n+1）位尾数的值。

2.1.3 E的余码表示

       E的位数决定了数的表示范围。IEEE标准采用E的编码约定。如果阶码E采用n位表示，那么对于阶码在二进制补码表示法的基础上再加上（2^(n-1)-1）构成它的余码表示法。二进制补码采取这样的方式来表示，即负数的补码是通过对值的每位取反，在加1得到的。例如对于三位阶码来说，在使用补码之后加上的是2^(3-1)-1 = 011.即-1的余码是先去符号，为001，然后取反110，加1 111，然后加上011，结果为010.

      余码的表示法的好处是在无符号比较器上可以用来比较有符号数。如图


图12.按照二进制补码的顺序排序的余3码（二进制补码、十进制值、余码表示）

如上图，当吧3位余码看做无符号数时，它可以表示的数值在-3到3之间单调递增。-3的余码是000,3的余码是110，所以如果用无符号比较器比较余3码001和100的大小，那么001要比100小，这对于硬件实现是一个期望的属性，因为无符号比较器比有符号比较器小，而且速度快。


图13 按照余码的顺序排序余3码

现在用6位的格式来表示0.5D：

0.5D=001000，其中S=0，E=010，M=1.00（或者00，省略1.）

也就是说0.5D的6位表示形式是001000.

一般来说，如果采取的是规格化的尾数和余码表示的阶码，那么带n位阶码的数对应的值为：


2.2能表示的数

         一种数值格式所能表示的数是指在这种格式中能精确表示的数值。例如，3位无符号整形表示格式所能的数如下图


图14 3位无符号格式所能表示的数；数轴表示法

在这种格式中，-1和9都不能表示。下图为目前为止能表示的所有数值及两种变换形式。为了有效控制表的大小，采用的格式是5位。这种格式中包含位符号位S，2位阶码（采用1位余码编码），2位尾数M（省略了1.）。非零列给出了目前为止讨论过的格式所能表示的数。这种格式不能表示0这个数。


图15 非0、下溢出和非规范化格式所有表示的数


图16非0表示法的表示的数

在图16中就是图15中的能够表示的正数，负数是相对0对称的，所以省略。从这个数轴中，的5条结论：

       1）能表示的数之间的主间隔取决于阶码位。在图16中，在0的每一边都有3个主间隔，这是因为阶码有两位。主间隔基本上位于2的幂之间，由于阶码有2位，因此这两位阶码可以形成3个不停的幂（2^-1=0.5D,2^0=1.0D,2^1=2.0D），每个都是能表示的数的间隔的开始。

       2）每个间隔中能表示的数的个数取决于尾数位。如果尾数位2位，则在每个间隔内能表示4个数。一般而言，如果尾数为N位，则在每个间隔内能表示2^N个数。如果一个要表示的值落在其中的一个间隔内，舍入后它就成为间隔内一个可表示的数。所以，每个间隔内能表示的数的个数越多，在这个区间中要表示的数也越精确。所以，尾数的位数决定了表示的精度。

      3）这种格式不能表示0.。

      4）越靠近0的地方，可表示的数离得越近。向0方向移动，每个间隔是前一个间隔的一半。能表示的数之间离得越近，表示数时就越准确。

      5）在0附近，能表示的数出现了空白。这是因为规格化的尾数对应的范围已经把0排除掉。当要表示0-0.5之间的数时，引入的误差远远大于（4倍）0.5-1.0这些更大的数之间引入的误差。一般而言，在这种表示法中，如果尾数的位数为N，那么接近0的间隔内引入的误差要比下一个间隔内引入的误差大2^N倍。对于基于非常小的数据值且依赖于收敛条件的准确检测的数值算法，这种缺陷会导致执行时不稳定且结果不准确。比如在除法中，表示这些很小的数时引入的误差在除法过程中，会被放大很多倍。

      在规格化浮点数系统中，一种可以容纳0的方法是下溢出（abrupt underflow）约定，如图15第2列，当阶码位E=0是，对应的数是0。在5位的格式中，这种方法在0附近提出了8个能表示的数（4个正数 4个负数），把它们都设置成0.虽然这种方法能表示0，但在0附近能表示的数之间出现了更大的空白，如下图：


图17 下溢出格式能表示的数

显然，与图16相比，这些数值算法的准确性依赖于0附近较小数的准确表示。


图18非规格化格式能表示的数

       实际上，IEEE标准采用的是非规格化的方法。这种方法在0附近放宽了对规格化数据的要求。如图15，当E=0时，不再假定尾数是1.xx这种形式，假定它是0.xx形式。阶码所表示的值与之前的间隔仍相同。例如，在图15中，非规格化数00001的阶码位00，尾数位01,。假定尾数是0.01，而阶码值与之前间隔的阶码值相同：是0而不是-1.也就是说，00001所表示的值是0.01*2^0=2^(-2)。图18展示了非规格化格式所能表示的数。现在，在0附近所能表示的数之间是均匀的间隔。从直观上看，非规格化约定吧在非0表示法中最后一个间隔内能表示的4个数分散开，用于弥补0附近的空白区域。这样就可以消除前面两种方法中不期望出现的空白。注意，在最后两个间隔中能表示的数之间的距离相同。一般而言，如果n位的阶码是0，那么这个值是


总之，浮点数表示法的精度取决于引入的最大误差，这种误差是要表示一个浮点数时，将它转化成能够表示的数而引入的误差。误差越小，精度越高。尾数的位数增加，可以提供浮点数表示法的精度。将图18中能表示的数增加一位尾数，可以将最大误差减少一半，从而提高精度。因此，一个数值系统用的尾数的位数越多，它的精度也越高。

2.3 特殊的位模式与精度

      在IEEE中，如果阶码的所有位都是1，而尾数为0，那么该数表示一个无穷大的值。如果尾数不为0，表示NaN（not a number）。


其他的所有数都是规格化的浮点数。单精度的浮点数有1位符号位（S）、8位阶码（E）和23位尾数（M）。双精度的浮点数则有1位符号位（S）、11位阶码（E）和52位尾数（M）。与单精度的浮点数格式相比，由于双精度浮点数的尾数多29位，因此可表示的数的最大误差减少到1/(2^29)！由于多了3位阶码，因此双精度浮点数格式也扩展了该区间能够表示的数的个数。这也扩展了能表示的数值范围，能表示非常大和非常小的数。所有能表示的数都落在-（负无穷大）到+（正无穷大）之间。溢出时可能会出现（如除以0）。任何能表示的数除以+或者-，结果都是0.

       其输入值没有意义的操作产生一个NaN--如0/0、0*、/、-。它们也用于表示一些数据。这些数据在程序中没有正确的初始化。在IEEE标准中NaN有两种类型：signaling和quiet。signaling NaN（SNaN）通过最重要的尾数位清零来表示，而quiet NaN则通过设置最重要的尾数位来表示。如果用一个signaling NaN作为算术运算操作的输入，则会引起异常。例如（1.0+SNaN）。只有说程序员希望在浮点计算过程中确保，只要使用任何NaN值时，程序就会中断执行，这时候signaling NaN才能使用。这时候意味着程序在执行过程中出错了。在执行关键任务的应用程序中，在没有通过单独的方式进行执行过程的有效性验证前，这种执行过程没法继续。例如，软件工程师常用signaling NaN来标记所有没有初始化的数据。这种实践确保在执行过程中检测所使用的哪些数据没有初始化。目前的GPU硬件（2010年的时候）还不支持signaling NaN，这是因为在大规模并行执行过程中，支持准确的signaling NaN是有难度的。

       quiet NaN在作为算术运算操作的输入时得到的额结果是quiet NaN。例如，操作（1.0+quiet NaN）得到的结果是一个quiet NaN。quiet NaN 常用的如下的应用程序中，：用户可以预先知道程序的运行结果，并在通过不同的输入时可以得到更为有效的结果时，决定是否要重新运行程序。当输出结果时，quiet NaN是以“NaN”的形式输出，因此用户就能在输出文件中很容易的找到它们。这也是在CUDA应用程序中如何检测数据损坏的方法。

2.4算术运算的准确度和舍入

         当两个输入操作上要进行浮点加法或减法运算时，而两者的阶码不同时，需要将阶码较小的操作数的尾数右移。比如（1.00*2^1）+(1.00*2^(-2))。因为阶码不同，所以先右移结果为（1.00*2^1）+(0.001*2^1)。理论上结果为1.001，可是在5位表示法中没法表示，所以按照硬件的位数，最好的就是产生一个非常接近的表示的数1.01*2^1.这样引入了（0.001*2^1）的误差。这个误差是最低有效位的一般。称这个误差为0.5 ULP（unit in the last place）.如果用于算术和舍入操作的硬件很完美，那么引入最大的误差也不会超过0.5 ULP。不过在某些更复杂的硬件算术运算单元中，都是通过迭代逼近算法实现的，比如分段函数和超越函数。不过在更新的GPU中算术运算显然要比之前的准确的多。

2.5算法的优化

     在计算数学中有个知识点就是在比如求和的时候，需要先对不同的数进行排序，或者说对不同的数的加法有顺序，这就是为了防止大数吃小数，先让最小的数不断的相加，就像滚雪球一样，这样才能不使得最后的结果丢失了小数，这也是为什么在大规模的并行数值算法中频繁的使用排序的重要原因。



ps：1、SIMT和SIMD是两种不同的实现技术，在CPU中的SSE（streaming SIMD extensions，单指令多数据扩展）。在SSE实现中，一条指令处理多个数据元素，先收集数据再将数据打包到单个的寄存器中，如果有严格对齐需求，则程序员需要用到额外的指令。在SIMT中，所有线程都在自己的寄存器中处理数据，减轻了数据收集和打包方面的负担。在SIMT中，选择控制流分支相对于SSE也要简单很多，不需要程序员用额外的指令来处理控制流。

        2、（书本当前时间）未来的cuda社保会为全局存储器中的数据提供更大容量的片上缓存，这种缓存会自动合并kernel函数中的访问模式，程序员会很少需要手动重新排列它们的访问模式，不过在目前的1亿次和仍在使用的上一代cuda设备中，这种合并技术仍然还是很有用的。

       3、Ryoo, S., Rodrigues, C. I., Baghsorkhi, S. S., & Stone, S. S. (2008). Optimization principles and application performance evaluation of a multithreaded GPU using CUDA. In Proceedings of the 13th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (pp. 73–82). Salt Lake City, UT. 

             Ryoo, S., Rodrigues, C. I., Stone, S. S., Baghsorkhi, S. S., Ueng, S. Z.,Stratton, J. A., et al. (2008). Program pruning for a multithreaded GPU. In Code generation and optimization: Proceedings of the Sixth Annual IEEE/ACM International Symposium on code generation and optimization (pp. 195–204). Boston, MA.








posted @ 2015-07-09 19:59  仙守  阅读(1615)  评论(0)  编辑  收藏  举报