0_Simple__fp16ScalarProduct

▶ 使用cuda内置无符号整数结构（__half2）及其汇编函数，计算两个向量的内积。

▶ 源代码

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "cuda_fp16.h"
#include "helper_cuda.h"

// 将数组 v 进行二分规约加法，使用 __forceinline__ 强制内联
__forceinline__ __device__ void reduceInShared(half2 * const v)
{
    if (threadIdx.x < 64)
        v[threadIdx.x] = __hadd2(v[threadIdx.x], v[threadIdx.x + 64]);
    __syncthreads();
    for (int i = 32; i > 0; i /= 2)
    {
        if (threadIdx.x < 32)
            v[threadIdx.x] = __hadd2(v[threadIdx.x], v[threadIdx.x + i]);
        __syncthreads();
    }
}

// 将数组 a 与 b 相加后进行规约加法，输入还包括指向结果的指针 h_result 及数组大小
__global__ void scalarProductKernel(half2 const * const a, half2 const * const b, float * const h_result, size_t const size)
{    
    __shared__ half2 shArray[128];
    const int stride = gridDim.x * blockDim.x;

    shArray[threadIdx.x] = __float2half2_rn(0.f);                               // 浮点数转无符号整数，这里相当于初始化为 0
    
    half2 value = __float2half2_rn(0.f);                                        
    for (int i = threadIdx.x + blockDim.x + blockIdx.x; i < size; i += stride)  // 半精度混合乘加，value = a[i] * b[i] + value
        value = __hfma2(a[i], b[i], value);                                     
    shArray[threadIdx.x] = value;
    __syncthreads();

    reduceInShared(shArray);                                                    // 规约得 a 和 b 的内积，因为使用了内联，共享内存指针可以传入

    if (threadIdx.x == 0)                                                       // 0 号线程负责写入结果
    {
        half2 result = shArray[0];
        h_result[blockIdx.x] = (float)(__low2float(result) + __high2float(result));
    }
}

void generateInput(half2 * a, size_t size)                                      // 生成随机数组
{
    for (size_t i = 0; i < size; ++i)
    {
        unsigned temp = rand();
        temp &= 0x83FF83FF;                                                     // 2214560767(10), 10000011111111111000001111111111(2)
        temp |= 0x3C003C00;                                                     // 1006648320(10),   111100000000000011110000000000(2)
        a[i] = *(half2*)&temp;
    }
}

int main(int argc, char *argv[])
{
    srand(time(NULL));
    const int blocks = 128, threads = 128;
    size_t size = blocks * threads * 16;

    int devID = 0;
    cudaDeviceProp devProp;
    cudaGetDeviceProperties(&devProp, devID);
    if (devProp.major < 5 || (devProp.major == 5 && devProp.minor < 3))
    {
        printf("required GPU with compute SM 5.3 or higher.
");
        return EXIT_WAIVED;
    }

    half2 *h_vec[2], *d_vec[2];
    float *h_result, *d_result;
    for (int i = 0; i < 2; ++i)
    {
        cudaMallocHost((void**)&h_vec[i], size * sizeof*h_vec[i]);
        cudaMalloc((void**)&d_vec[i], size * sizeof*d_vec[i]);
    }
    cudaMallocHost((void**)&h_result, blocks * sizeof*h_result);
    cudaMalloc((void**)&d_result, blocks * sizeof*d_result);
    for (int i = 0; i < 2; ++i)
    {
        generateInput(h_vec[i], size);
        cudaMemcpy(d_vec[i], h_vec[i], size * sizeof*h_vec[i], cudaMemcpyHostToDevice);
    }
    scalarProductKernel << <blocks, threads >> >(d_vec[0], d_vec[1], d_result, size);
    cudaMemcpy(h_result, d_result, blocks * sizeof * h_result, cudaMemcpyDeviceToHost);
    cudaDeviceSynchronize();

    float result = 0;
    for (int i = 0; i < blocks; ++i)
        result += h_result[i];
    printf("Result: %f 
", result);

    for (int i = 0; i < 2; ++i)
    {
        cudaFree(d_vec[i]);
        cudaFreeHost(h_vec[i]);
    }
    cudaFree(d_result);
    cudaFreeHost(h_result);
    getchar();
    return EXIT_SUCCESS;
}

● 输出结果

GPU Device 0: "GeForce GTX 1070" with compute capability 6.1

Result: 853856.000000

▶ 涨姿势

● CUDA 无符号半精度整数，就是用 unsigned short 对齐到 2 Byte 来封装的

typedef struct __align__(2) { unsigned short x; } __half;

typedef struct __align__(4) { unsigned int x; } __half2; 

#ifndef CUDA_NO_HALF
typedef __half half;
typedef __half2 half2;
#endif

● 关于 __inline__ 和 __forceinline__

参考stackoverflow。https://stackoverflow.com/questions/19897803/forceinline-effect-at-cuda-c-device-functions

与C中__forceinline__类似，忽略编译器的建议，强制实现内联函数。如果函数只调用累次那么优化没有效果，但是如果调用了多次（如内联函数出现在循环中），则会产生明显的提升。另外，在递归中一般不用。

● 关于 __CUDACC__ 和 __CUDA_ARCH__

■ 参考 stackoverflow【https://stackoverflow.com/questions/8796369/cuda-and-nvcc-using-the-preprocessor-to-choose-between-float-or-double】

■ __CUDACC__ 使用 nvcc 进行编译时有定义。

■ __CUDA_ARCH__ 编译主机代码时无定义（无论是否使用 nvcc）；编译设备代码时有定义，且值等于编译命令指定的计算能力号。

■ 范例代码：（为了方便查看，使用了缩进）

#ifdef __CUDACC__
    #warning using nvcc

    template <typename T>　　　　　　　　　　　　　　　　　　// 一般的核函数
    __global__ void add(T *x, T *y, T *z)
    {
        int idx = threadIdx.x + blockDim.x * blockIdx.x;
        z[idx] = x[idx] + y[idx];
    }

    #ifdef __CUDA_ARCH__
        #warning device code trajectory
        #if __CUDA_ARCH__ > 120
            #warning compiling with datatype double
            template void add<double>(double *, double *, double *);
        #else
            #warning compiling with datatype float
            template void add<float>(float *, float *, float *);
        #endif
    #else
        #warning nvcc host code trajectory
    #endif
#else
    #warning non - nvcc code trajectory
#endif

■ 编译及输出结果

$ ln -s cudaarch.cu cudaarch.cc
$ gcc -c cudaarch.cc -o cudaarch.o
cudaarch.cc:26:2: warning: #warning non-nvcc code trajectory

$ nvcc -arch=sm_11 -Xptxas="-v" -c cudaarch.cu -o cudaarch.cu.o
cudaarch.cu:3:2: warning: #warning using nvcc
cudaarch.cu:14:2: warning: #warning device code trajectory
cudaarch.cu:19:2: warning: #warning compiling with datatype float
cudaarch.cu:3:2: warning: #warning using nvcc
cudaarch.cu:23:2: warning: #warning nvcc host code trajectory
ptxas info    : Compiling entry function '_Z3addIfEvPT_S1_S1_' for 'sm_11'
ptxas info    : Used 4 registers, 12+16 bytes smem

$ nvcc -arch=sm_20 -Xptxas="-v" -c cudaarch.cu -o cudaarch.cu.o
cudaarch.cu:3:2: warning: #warning using nvcc
cudaarch.cu:14:2: warning: #warning device code trajectory
cudaarch.cu:16:2: warning: #warning compiling with datatype double
cudaarch.cu:3:2: warning: #warning using nvcc
cudaarch.cu:23:2: warning: #warning nvcc host code trajectory
ptxas info    : Compiling entry function '_Z3addIdEvPT_S1_S1_' for 'sm_20'
ptxas info    : Used 8 registers, 44 bytes cmem[0]

● 用到的汇编函数

// 表明主机和设备共有代码
#define __CUDA_FP16_DECL__ __host__ __device__

// 浮点数转无符号整数
__CUDA_FP16_DECL__ __half2 __float2half2_rn(const float f)
{
    __half2 val;
    asm("{.reg .f16 low;
"
        "  cvt.rn.f16.f32 low, %1;
"
        "  mov.b32 %0, {low,low};}
" : "=r"(val.x) : "f"(f));
    return val;
}

// 计算无符号整数 a + b
#define BINARY_OP_HALF2_MACRO(name)                                             
    do                                                                          
    {                                                                           
        __half2 val;                                                            
        asm("{"#name".f16x2 %0,%1,%2;
}" :"=r"(val.x) : "r"(a.x), "r"(b.x));   
        return val;                                                             
    }                                                                           
    while(0);                                                                   

__CUDA_FP16_DECL__ __half2 __hadd2(const __half2 a, const __half2 b)
{
    BINARY_OP_HALF2_MACRO(add);
}

// 计算无符号整数 a * b + c
#define TERNARY_OP_HALF2_MACRO(name)                                                        
    do                                                                                      
    {                                                                                       
        __half2 val;                                                                        
        asm("{"#name".f16x2 %0,%1,%2,%3;
}" : "=r"(val.x) : "r"(a.x), "r"(b.x), "r"(c.x)); 
        return val;                                                                         
    }                                                                                       
    while(0);                                                                               

__CUDA_FP16_DECL__ __half2 __hfma2(const __half2 a, const __half2 b, const __half2 c)
{
    TERNARY_OP_HALF2_MACRO(fma.rn);
}

// 将无符号整数的低 2 字节转化为浮点数
__CUDA_FP16_DECL__ float __low2float(const __half2 l)
{
    float val;
    asm("{.reg .f16 low,high;
"
        "  mov.b32 {low,high},%1;
"
        "  cvt.f32.f16 %0, low;}
" : "=f"(val) : "r"(l.x));
    return val;
}

// 将无符号整数的高 2 字节转化为浮点数
__CUDA_FP16_DECL__ float __high2float(const __half2 l)
{
    float val;
    asm("{.reg .f16 low,high;
"
        "  mov.b32 {low,high},%1;
"
        "  cvt.f32.f16 %0, high;}
" : "=f"(val) : "r"(l.x));
    return val;
}