GPU编程和流式多处理器（三）

GPU编程和流式多处理器（三）

3. Floating-Point Support

快速的本机浮点硬件是GPU的存在理由，并且在许多方面，它们在浮点实现方面都等于或优于CPU。全速支持异常可以根据每条指令指定直接舍入，特殊功能单元可为六种流行的单精度先验函数，提供高性能的近似函数。相比之下，x86 CPU在微代码中实现异常，其运行速度可能比在规范化浮点算子上运行的速度慢100倍。舍入方向是由一个控制字指定的，该控制字需要数十个时钟周期来更改，并且SSE指令集中唯一的超越逼近函数是用于倒数和倒数平方根的，给出必须用牛顿-微分精确的12位近似值。使用Raphson迭代之前。

由于GPU的核心数量较多，但其较低的时钟频率会有所抵消，开发人员可以期望在公平竞争的环境下最多提高10倍（或大约10倍）的速度。如果一篇论文报告了从将优化的CPU实现移植到CUDA所实现的100倍或更高的速度，则可能是上述“指令集不匹配”之一。

3.1. Formats

图2描述了CUDA支持的三（3）个IEEE标准浮点格式：双精度（64位），单精度（32位）和半精度（16位）。值分为三个字段：正负号，指数和尾数。对于double，single和half，指数字段的大小分别为11位，8位和5位。相应的尾数字段为52、23和10位。

 

 

 Figure 2 Floating-point formats.

指数字段更改浮点值的解释。最常见的（“正常”）表示形式将一个隐式1位编码为尾数，然后将该值乘以2个e-bias，其中bias是在编码为浮点表示形式之前添加到实际指数的值。例如，单精度的偏差为127。

表7总结了如何对浮点值进行编码。对于大多数指数值（所谓的“正常”浮点值），尾数假定为隐式1，并将其乘以指数的偏置值。最大指数值保留用于无穷大和非数字值。除以零（或除以除法）会产生无穷大；执行无效运算（例如取平方根或负数的对数）会产生NaN。最小指数值保留给太小而无法用隐式前导1表示的值。随着所谓的正规数接近零，它们失去有效精度的位，这种现象称为逐渐下溢。表8给出了三种格式的某些极值的编码和值。

Table 7 Floating-Point Representations

 

 

 

 

 Table 8 Floating-Point Extreme Values

 

 

 

 

 

 

 Rounding舍入取整法

IEEE标准提供了四（4）个回合模式。

  • 四舍五入到最近（也称为“四舍五入”）
  • 向零舍入（也称为“截断”或“斩”）
  • 向下舍入（或“朝负无穷大舍入”）
  • 向上舍入（或“朝正无穷大舍入”）

舍入取整法是迄今为止最常用的舍入模式，在每次操作后将中间值舍入为最接近的可表示浮点值。向上舍入和向下舍入（“定向舍入模式”）用于间隔算术，其中一对浮点值用于括住计算的中间结果。为了正确地将结果括起来，该间隔的下限值和上限值必须分别四舍五入为负无穷大（“向下”）和正无穷大（“ up”）。

C语言不提供任何方式来按指令指定舍入模式，并且CUDA硬件不提供用于隐式指定舍入模式的控制字。因此，CUDA提供了一组内部函数来指定操作的舍入模式，如表9所示。

Table 9 Intrinsics for Rounding

 

 

 

 

 Conversion

通常，开发人员可以使用标准C结构，在不同的浮点表示形式和/或整数之间进行转换：隐式转换或显式类型转换。但是，如有必要，开发人员可以使用表10中列出的内在函数，执行C语言规范中不包含的转换，例如，使用定向舍入的转换。

 

 

 

 

 由于一半未在C编程语言中标准化，因此CUDA在__half2float（）和__float2half（）的接口中使用unsigned short。__float2half（）仅支持舍入到最近舍入模式。

float __half2float( unsigned short );

unsigned short __float2half( float );

3.2. Single Precision (32-Bit)

单精度浮点支持是GPU计算的主力军。GPU已经过优化，可以在此数据类型上原生提供高性能，11不仅适用于核心标准IEEE操作（例如加法和乘法），还适用于非标准操作（例如对先验的近似（例如sin（）和log（）））。32位值与整数保存在同一寄存器文件中，因此单精度浮点值和32位整数（使用__float_as_int（）和__int_as_float（））之间的强制转换是免费的。

加法，乘法和乘加

编译器自动将浮点值的+，–和*运算符转换为加，乘和乘加指令。__fadd_rn（）和__fmul_rn（）内部函数，可以被用于加入抑制融合和乘法操作进入乘加指令。

对等与除法

对于具有2.x或更高版本的计算能力的设备，使用--prec-div = true编译代码时，除法运算符符合IEEE规定。对于具有计算能力1.x的设备或对于具有计算能力2.x的设备，当使用--prec-div = false编译代码时，除法运算符和__fdividef（x，y）具有相同的精度，

 

 

 SM中的特殊功能单元（SFU）实现了六个常见超越功能的快速版本。

  • 正弦和余弦
  • 对数和指数
  • 倒数和倒数平方根

表11（摘自有关Tesla架构的论文）摘录了受支持的操作和相应的精度。SFU并不能完全实现精度，可以很好地近似于这些功能，并且速度很快。对于比优化的CPU等效端口（例如25倍或更多）快得多的CUDA端口，代码很可能依赖SFU。

Table 11 SFU Accuracy

 

 

 使用表12中给出的内在函数访问SFU。指定--fast-math编译器选项，将使编译器将常规C运行时runtime，调用替换为上面列出的相应SFU内部函数。

Table 12 SFU Intrinsics

 

 

 Miscellaneous

__saturate(x) returns 0 if x<0, 1 if x>1, and x otherwise.

3.3. Double Precision (64-Bit)

带有SM 1.3的CUDA中添加了双精度浮点支持（最早在GeForce GTX 280中实现），而SM 2.0则提供了大大改进的双精度支持（功能和性能）。CUDA对双精度的硬件支持具有全速异常功能，并且从SM 2.x开始，符合IEEE 754 c的本机融合乘法加法指令（FMAD）。2008年，仅执行一个舍入步骤。FMAD除了是本质上有用的操作之外，还可以使某些功能（与Newton-Raphson迭代收敛）完全准确。

与单精度运算一样，编译器会自动将标准C运算符转换为乘法，加法和乘加指令。__dadd_rn（）和__dmul_rn（）内部函数可以被用于加入抑制融合和乘法操作进入乘加指令。

3.4. Half Precision (16-Bit)

对于5位指数和10位有效数，半值具有足够的精度以用于HDR（高动态范围）图像，并可用于保存其它不需要浮点精度的值，例如角度。半精度值仅用于存储，而不用于计算，硬件仅提供指令以转换为32位或从32位转换。这些指令以__halftofloat（）和__floattohalf（）内部函数公开。

float __halftofloat( unsigned short );

unsigned short __floattohalf( float );

这些内部函数使用unsigned short，因为C语言尚未标准化半浮点类型。

3.5. Case Study: float→half Conversion

研究float → half转换操作，了解浮点编码和舍入细节的有用方法。这是一个简单的一元运算，可以专注于编码和舍入，而不会被浮点运算的细节和中间表示的精度所分散。

当从float转换为half时，对于任何太大而无法表示的float的正确输出是half infinity。任何太小而不能代表一半（甚至是反常的一半）的浮点数都必须固定为0.0。舍入为0.0的一半的最大浮点数为0x32FFFFFF或2.98，而舍入为无穷大的一半的最小浮点数为65520.0。此范围内的float值可以通过传播符号位转换为一半，重新偏置指数（因为float的8位指数偏差为127，一半的5位指数偏差为15），然后将float的尾数四舍五入到最接近的一半尾数。在所有情况下，舍入都是简单的，除非输入值恰好落在两个可能的输出值之间。在这种情况下，IEEE标准指定四舍五入到“最近的偶数”值。在十进制算术中，这意味着四舍五入为1.5到2.0，但也四舍五入为2.5到2.0，以及（例如）四舍五入到0.5到0.0。

清单3显示了一个C例程，该例程完全复制了CUDA硬件实现的浮点到一半转换操作。变量exp和mag包含输入指数和“幅值”，尾数和指数以及被屏蔽的符号位。可以对幅度执行许多操作，例如比较和舍入操作，而无需分离指数和尾数。

清单3中使用的宏LG_MAKE_MASK创建具有给定位数的掩码：#define LG_MAKE_MASK（bits）（（1 << bits）-1）。联合用于治疗相同的32位值作为浮子和无符号整型; 诸如*（（float *）（＆u））之类的习惯用法不是可移植的。该例程首先传播输入符号位并将其屏蔽掉。

提取幅度和指数后，该函数处理输入浮点为INF或NaN的特殊情况，并尽早退出。注意，INF是带符号的，但NaN具有规范的无符号值。将输入浮点值钳位到与可表示的半值相对应的最小值或最大值，然后重新计算钳位值的大小。不要被构造f32MinRInfin和f32MaxRf16_zero的复杂代码所迷惑；它们是分别具有值0x477ff000和0x32ffffff的常量。

例程的其余部分处理输出正常和异常的情况（输入异常在前面的代码中被钳位，因此mag对应于正常float）。与钳位代码一样，f32Minf16Normal是一个常量，其值为0x38ffffff。

要构建法线，必须计算新的指数（第92和93行），正确舍入的10位尾数将移入输出。要构造非正规数，必须将隐式1与输出尾数进行“或”运算，然后将所得尾数偏移与输入指数对应的量。对于正向和非正向，输出尾数的舍入分两步完成。舍入是通过添加一个1的掩码来实现的，该掩码的末尾刚好等于输出的LSB，如图3所示。

 

 

 Figure 3 Rounding mask (half).

如果输入的位12置1，则此操作将增加输出尾数；否则，将增加输出尾数。如果输入尾数全为1，则溢出会导致输出指数正确递增。如果向此调整的最高有效位再加1，将获得小学风格的四舍五入，其中平局决胜数更大。反之，即使设置了舍入到最接近值，如果设置了10位输出的LSB，则有条件地增加输出尾数（图8.4）。注意，这些步骤可以按任何顺序执行，也可以以许多不同的方式重新制定。

 

 

 Figure 4 Round-to-nearest-even (half).

Listing 3. ConvertToHalf().

/*

 * exponent shift and mantissa bit count are the same.

 *    When we are shifting, we use [f16|f32]ExpShift

 *    When referencing the number of bits in the mantissa,

 *        we use [f16|f32]MantissaBits

 */

const int f16ExpShift = 10;

const int f16MantissaBits = 10;

 

const int f16ExpBias = 15;

const int f16MinExp = -14;

const int f16MaxExp = 15;

const int f16SignMask = 0x8000;

 

const int f32ExpShift = 23;

const int f32MantissaBits = 23;

const int f32ExpBias = 127;

const int f32SignMask = 0x80000000;

 

unsigned short

ConvertFloatToHalf( float f )

{

    /*

     * Use a volatile union to portably coerce

     * 32-bit float into 32-bit integer

     */

    volatile union {

        float f;

        unsigned int u;

    } uf;

    uf.f = f;

 

    // return value: start by propagating the sign bit.

    unsigned short w = (uf.u >> 16) & f16SignMask;

 

    // Extract input magnitude and exponent

    unsigned int mag = uf.u & ~f32SignMask;

    int exp = (int) (mag >> f32ExpShift) - f32ExpBias;

 

    // Handle float32 Inf or NaN

    if ( exp == f32ExpBias+1 ) {    // INF or NaN

 

        if ( mag & LG_MAKE_MASK(f32MantissaBits) )

            return 0x7fff; // NaN

 

        // INF - propagate sign

        return w|0x7c00;

    }

 

    /*

     * clamp float32 values that are not representable by float16

     */

    {

        // min float32 magnitude that rounds to float16 infinity

 

        unsigned int f32MinRInfin = (f16MaxExp+f32ExpBias) <<

            f32ExpShift;

        f32MinRInfin |= LG_MAKE_MASK( f16MantissaBits+1 ) <<

            (f32MantissaBits-f16MantissaBits-1);

 

        if (mag > f32MinRInfin)

            mag = f32MinRInfin;

    }

 

    {

        // max float32 magnitude that rounds to float16 0.0

 

        unsigned int f32MaxRf16_zero = f16MinExp+f32ExpBias

            (f32MantissaBits-f16MantissaBits-1);

        f32MaxRf16_zero <<= f32ExpShift;

        f32MaxRf16_zero |= LG_MAKE_MASK( f32MantissaBits );

 

        if (mag < f32MaxRf16_zero)

            mag = f32MaxRf16_zero;

    }

 

    /*

     * compute exp again, in case mag was clamped above

     */

    exp = (mag >> f32ExpShift) - f32ExpBias;

 

    // min float32 magnitude that converts to float16 normal

    unsigned int f32Minf16Normal = ((f16MinExp+f32ExpBias)<<

        f32ExpShift);

    f32Minf16Normal |= LG_MAKE_MASK( f32MantissaBits );

    if ( mag >= f32Minf16Normal ) {

        //

        // Case 1: float16 normal

        //

 

        // Modify exponent to be biased for float16, not float32

        mag += (unsigned int) ((f16ExpBias-f32ExpBias)<<

            f32ExpShift);

 

        int RelativeShift = f32ExpShift-f16ExpShift;

 

        // add rounding bias

        mag += LG_MAKE_MASK(RelativeShift-1);

 

        // round-to-nearest even

        mag += (mag >> RelativeShift) & 1;

 

        w |= mag >> RelativeShift;

    }

    else {

        /*

         * Case 2: float16 denormal

         */

 

        // mask off exponent bits - now fraction only

        mag &= LG_MAKE_MASK(f32MantissaBits);

 

        // make implicit 1 explicit

        mag |= (1<<f32ExpShift);

 

        int RelativeShift = f32ExpShift-f16ExpShift+f16MinExp-exp;

 

        // add rounding bias

        mag += LG_MAKE_MASK(RelativeShift-1);

 

        // round-to-nearest even

        mag += (mag >> RelativeShift) & 1;

 

        w |= mag >> RelativeShift;

    }

    return w;

}

实际上，开发人员应使用__floattohalf（）内在函数将float转换为一半，编译器会将其转换为单个F2F机器指令。提供此示例例程的目的纯粹，为了帮助理解浮点布局和舍入。同样，检查所有INF / NAN和非标准值的特殊情况代码，有助于说明IEEE规范的这些功能，自诞生以来就一直引起争议：使硬件变慢，成本更高，或者由于增加的硅片面积和工程设计而使硬件变慢验证工作。

清单3中的ConvertFloatToHalf（）例程被合并到名为float_to_float16.cu的程序中，该程序针对每个32位浮点值测试其输出。

3.6. Math Library

CUDA包括一个以C运行时库为模型的内置数学库，但有一些小区别：CUDA硬件不包括舍入模式寄存器（取而代之的是，舍入模式是按指令编码的），作为引用当前舍入模式的rint（），始终舍入为最接近值。此外，硬件不会引发浮点异常。异常运算的结果（例如，取负数的平方根）将编码为NaN。

表13列出了数学库函数以及每个函数的最大误差（单位为ulps）。在float上操作的大多数函数的函数名称后均带有“ f”，例如，计算正弦函数的函数如下。

double sin( double angle );

float sinf( float angle );

Table 13 Math Library

     

ULP ERROR

FUNCTION

OPERATION

EXPRESSION

32

64

x+y

Addition

x+y

01

0

x/y

Multiplication

x*y

01

0

x/y

Division

x/y

22

0

1/x

Reciprocal

1/x

12

0

acos[f](x)

Inverse cosine

cos–1 x

3

2

acosh[f](x)

Inverse hyperbolic cosine

  

4

2

asin[f](x)

Inverse sine

sin–1x

4

2

asinh[f](x)

Inverse hyperbolic sine

  

3

2

atan[f](x)

Inverse tangent

tan–1x

2

2

atan2[f](y,x)

Inverse tangent of y/x

  

3

2

atanh[f](x)

Inverse hyperbolic tangent

tanh–1

3

2

cbrt[f](x)

Cube root

  

1

1

ceil[f](x)

“Ceiling,” nearest integer greater than or equal to x

⌊x⌋

0

  

copysign[f](x,y)

Sign of y, magnitude of x

  

n/a

  

cos[f](x)

Cosine

cos x

2

1

cosh[f](x)

Hyperbolic cosine

  

2

  

cospi[f](x)

Cosine, scaled byΠ

cosΠx

2

  

erf[f](x)

Error function

  

3

2

erfc[f](x)

Complementary error function

  

6

4

erfcinv[f](y)

Inverse complementary error function

Return x for which y=1-erff(x)

7

8

erfcx[f](x)

Scaled error function

ex2[erff(x))

6

3

erfinv[f](y)

Inverse error function

Return x for which y=erff(x)

3

5

exp[f](x)

Natural exponent

ex

2

1

exp10[f](x)

Exponent (base 10)

10x

2

1

exp2[f](x)

Exponent (base 2)

2x

2

1

expm1[f](x)

Natural exponent, minus one

ex–1

1

1

fabs[f](x)

Absolute value

|x|

0

0

fdim[f](x,y)

Positive difference

x–1,x>y+0,x≤NAN,X or y NAN

0

0

floor[f](x)

“Floor,” nearest integer less than or equal to x

⌊x⌋

0

0

fma[f](x,y,z)

Multiply-add

xy + z

0

0

fmax[f](x,y)

Maximum

x, x y or isNaN(y)y, otherwise

0

0

fmin[f](x,y)

Minimum

x, x y or isNaN(y)y, otherwise

0

0

fmod[f](x,y)

Floating-point remainder

  

0

0

frexp[f](x,exp)

Fractional component

  

0

0

hypot[f](x,y)

Length of hypotenuse

  

3

2

ilogb[f](x)

Get exponent

  

0

0

isfinite(x)

Nonzero if x is not ±INF

    

n/a

isinf(x)

Nonzero if x is ±INF

    

n/a

isnan(x)

Nonzero if x is a NaN

    

n/a

j0[f](x)

Bessel function of the first kind (n=0)

J0(x)

93

73

j1[f](x)

Bessel function of the first kind (n=1)

J1(x)

93

73

jn[f](n,x)

Bessel function of the first kind

Jn(x)

  

*

ldexp[f](x,exp)

Scale by power of 2

x2exp

0

0

lgamma[f](x)

Logarithm of gamma function

ln(Γ(x))

64

44

llrint[f](x)

Round to long long

  

0

0

llround[f](x)

Round to long long

  

0

0

lrint[f](x)

Round to long

  

0

*

lround[f](x)

Round to long

  

0

0

log[f](x)

Natural logarithm

ln(x)

1

1

log10[f](x)

Logarithm (base 10)

log10 x

3

1

log1p[f](x)

Natural logarithm of x+1

ln(x + 1)

2

1

log2[f](x)

Logarithm (base 2)

log2 x

3

1

logb[f](x)

Get exponent

  

0

0

modff(x,iptr)

Split fractional and integer parts

  

0

0

nan[f](cptr)

Returns NaN

NaN

  

n/a

nearbyint[f](x)

Round to integer

  

0

0

nextafter[f](x,y)

Returns the FP value closest to x in the direction of y

    

n/a

normcdf[f](x)

Normal cumulative distribution

  

6

5

normcdinv[f](x)

Inverse normal cumulative distribution

  

5

8

pow[f](x,y)

Power function

xy

8

2

rcbrt[f](x)

Inverse cube root

  

2

1

remainder[f](x,y)

Remainder

  

0

0

remquo[f](x,y,iptr)

Remainder (also returns quotient)

  

0

0

rsqrt[f](x)

rsqrt[f](x) Reciprocal

  

2

1

rint[f](x)

Round to nearest int

  

0

0

round[f](x)

Round to nearest int

  

0

0

scalbln[f](x,n)

Scale x by 2n (n is long int)

x2n

0

0

scalbn[f](x,n)

Scale x by 2n (n is int)

x2n

0

0

signbit(x)

Nonzero if x is negative

  

n/a

0

sin[f](x)

Sine

sin x

2

1

sincos[f](x,s,c)

Sine and cosine

*s=sin(x);*c=cos(x);

2

1

sincospi[f](x,s,c)

Sine and cosine

*s=sin(Πx);*c=cos(Πx);

2

1

sinh[f](x)

Hyperbolic sine

  

3

1

Sine, scaled by Π

sin Πx

  

2

1

sqrt[f](x)

Square root

  

36

0

tan[f](x)

Tangent

tan x

4

2

tanh[f](x)

Hyperbolic tangent

  

2

1

tgamma[f](x)

True gamma function

γ(x)

11

8

trunc[f](x)

Truncate (round to integer toward zero)

  

0

0

y0[f](x)

Bessel function of the second kind (n=0)

Y0(x)

93

73

y1[f](x)

Bessel function of the second kind (n=1)

Y1(x)

93

73

yn[f](n,x)

Bessel function of the second kind

Yn(x)

  

..

 

 

 这些在表13中表示为例如sin [f]。

 

posted @ 2021-01-05 08:52  吴建明wujianming  阅读(537)  评论(0)  编辑  收藏  举报