zoukankan      html  css  js  c++  java

  • 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]。

    人工智能芯片与自动驾驶
  • 相关阅读:
    Linux下截图工具
    Vue学习——router路由的实现原理
    Vue学习——vue的双向数据绑定原理
    JavaScript学习——面向对象（一）——创建对象（工厂模式和构造器模式）
    子组件给父组件的传值
    Vue组件
    JavaScript学习——事件对象Event
    JavaScript学习——事件处理程序
    JavaScript技巧——轮播图
    javascript——let关键字
  • 原文地址：https://www.cnblogs.com/wujianming-110117/p/14233691.html
Copyright © 2011-2022 走看看