# CUDA 的半精度浮點數運算

半精度浮點數，又稱為 float16，由於半精度浮點數的理論計算量比單精度浮點數少，所以理論上可以減少運算時間，但半精度浮點數需要硬體支援，所以一般的硬體無法使用。而較新的 Nvidia 的 GPU 則支援半精度浮點。

自 CUDA 7.5 後開始支援半精度浮點計算，它要求 GPU 的計算力 (compute capability) 至少大於等於 5.3，可自行到 [wiki](https://en.wikipedia.org/wiki/CUDA) 查詢自己 GPU 對應的計算力，或是直接在 CUDA 中通過 cudaDeviceProp 查詢。

    int dev_id = 0;
    cudaDeviceProp dev_prop;
    cudaGetDeviceProperties(&dev_prop, dev_id);
    printf("The compute capability is %d.%d\n", dev_prop.major, dev_prop.minor);

</br>

一般而言在 CUDA 設備上使用半精度浮點數有三個步驟

1. 在 host 端將 float32 轉成 float16 並將資料推入 GPU
2. 對 float16 做運算，有兩種狀況
* 存在可以使用的 float16 硬體加速指令，使用指令做運算
* 不存在可以使用或不想使用 float16 硬體加速指令，先轉成 float32 做運算，最後轉成 float16 輸出
3. 將 float16 資料推入 host，並在 host 端將 float16 轉成 float32

</br>

## 轉換浮點數型態

CUDA 的半精度浮點數使用 [IEEE 754 表準](https://en.wikipedia.org/wiki/Half-precision_floating-point_format)，格式如下。

<div id="fp16-format" align="center">
<br/>
<img src="https://i.imgur.com/EupG9Zc.png" width="384"/>
</div>

</br>

CUDA 內相互轉換 float16 和 float32 有以下三種方法

1. 在 kernel 內轉換
在 CUDA kernel 內轉換，好處是簡單且沒有版本問題，缺點是要多宣告更多 CUDA 的記憶體，且轉換只能在 GPU 上執行。

    __global__ void float_to_half(half *out, float *in) {
        int idx = threadIdx.x;
        half v0 = (half)(in[idx]);       // 直接轉換
        half v1 = __float2half(in[idx]); // 用指令轉換
        out[idx] = v0;
    }

2. 在 host 端轉換
可以直接在你的 C/C++ 程式碼中使用 ```__float2half``` 指令，但要求 CUDA 版本至少大於等於 9.2。
3. 使用第三方庫在 host 端轉換
這是最保守的方法，不需要考慮版本問題，也不用浪費記憶體。你可以使用我撰寫的庫 [simple half](https://github.com/CGLemon/simple_half) 完成轉換工作，它可以滿足基本需求。

</br>

## float16 精度與溢位的問題

由於 float16 表示範圍有限，它的確有較嚴重的精度和溢位問題，導致算出的數值有可感知的誤差，但大部分出來的數值差距小於 1% ，如果你的任務對於精度有要求，則應當放棄 float16 ，反之你的任務對於誤差較不敏感，例如神經網路的計算，則極力推薦 float16，因為在算力足夠的 GPU 上加速效果明顯。

## 設計 kernel

理論而言，我們可以用模板設計一個 float16 和 float32 共用的 kernel，如此不需要新增任何工作量，就能直接使用 float16 運算，如下

template <typename T>
__global__ void add(T *C, T *A, T *B) {
    int idx = threadIdx.x;
    T a = A[idx];
    T b = B[idx];
    C[idx] = a + b;
}

但實際上，由於 float16 屬於硬體加速指令，需要硬體、程式庫和編譯器的配合（指 ```gencode``` 選項），在當前大環境中還是需要支援較就的硬體，編譯器無法同時對所有的硬體優化，因此 kernel 本身不能設計成純粹的 float16 計算 （模板還是需要同時編譯 float16 和 float32，因此會有成純粹的 float16 kernel），需要 float16 和 float32 混合使用，例如

template <typename T>
__global__ void add(T *C, T *A, T *B) {
    int idx = threadIdx.x;
    float a = (float)(A[idx]); // 永遠都轉換成 float32
    float b = (float)(B[idx]); // 永遠都轉換成 float32
    float c = a + b; // 永遠用 float32 作用算
    C[idx] = (T)c; // 轉換成對應的型態 
}

但如果你只針對自己需要的平台，例如只設計給 40 系列的，那當然可以使用使用版本一的 kernel

</br>

而 blas 本身不受編譯器影響，可以根據當前需要切換

if (fp16) {
    cublasHgemm(...); // float16 的矩陣乘法
} else {
    cublasSgemm(...); // float32 的矩陣乘法
}

</br>

## 實例

假設有 A 、B 和 C 三個矩陣，執行的任務為 A = A + B 且 O = A x C ，假設加法為簡單運算，大部分時間都消耗在乘法上，那加法就不太需要對 float16 實做專門優化的版本，因此 ```add_kernel``` 是 float16 和 float32 共用的，而乘法則須使用專門的 float16 指令 ```cublasHgemm```。

// 檔案名稱 main.cu
// 編譯 nvcc main.cu -o demo -lcublas
// 執行 ./demo

#include <stdio.h>

#include "half.h" // https://github.com/CGLemon/simple_half

#include <cublas_v2.h>
#include <cuda_fp16.h>

int div_up(int a, int b) {
    return (a + b - 1) / b;
}

// float16 和 float32 通用的 kernel
template <typename T>
__global__ void add_kernel(T *C, T *A, T *B, int total_elements) {
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    if (idx < total_elements) {
        float a = (float)(A[idx]);
        float b = (float)(B[idx]);
        float c = a + b;
        C[idx] = (T)c;
    }
}

// float16 和 float32 通用的函數
template <typename T>
void add(T *C, T *A, T *B, int total_elements) {
    int block_size = 32;
    int blocks = div_up(total_elements, block_size);
    add_kernel<<<blocks, block_size>>>(C, A, B, total_elements);
}

// 宣告矩陣的大小
#define N (10)
#define SIZE (N*N)

void test_fp16(cublasHandle_t handle) {
    // host 端的矩陣
    half_float_t A[SIZE];
    half_float_t B[SIZE];
    half_float_t C[SIZE];
    half_float_t O[SIZE];

    for (int i = 0; i < SIZE; ++i) {
        A[i] = GetFp16(1);
        B[i] = GetFp16(2);
        C[i] = GetFp16(3);
    }

    // device (GPU) 端的矩陣
    half *A_gpu;
    half *B_gpu;
    half *C_gpu;
    half *O_gpu;

    size_t bufsize = SIZE * sizeof(half);
    cudaMalloc(&A_gpu, bufsize);
    cudaMalloc(&B_gpu, bufsize);
    cudaMalloc(&C_gpu, bufsize);
    cudaMalloc(&O_gpu, bufsize);

    // host -> device
    cudaMemcpy(A_gpu, A, bufsize, cudaMemcpyHostToDevice);
    cudaMemcpy(B_gpu, B, bufsize, cudaMemcpyHostToDevice);
    cudaMemcpy(C_gpu, C, bufsize, cudaMemcpyHostToDevice);

    // 執行 A = A + B 之矩陣加法
    add<half>(A_gpu, A_gpu, B_gpu, SIZE);

    // 執行 O = A x C 之矩陣乘法
    half_float_t alpha = GetFp16(1);
    half_float_t beta = GetFp16(0);
    cublasHgemm(handle,
                CUBLAS_OP_N,
                CUBLAS_OP_T,
                N, N, N,
                (half*)&alpha,
                C_gpu, N,
                A_gpu, N,
                (half*)&beta,
                O_gpu, N);

    // device -> host
    cudaMemcpy(O, O_gpu, bufsize, cudaMemcpyDeviceToHost);

    // 輸出答案
    for (int i = 0; i < SIZE; ++i) {
        printf("%.2f ", GetFp32(O[i]));
    }
    printf("\n");

    cudaFree(A_gpu);
    cudaFree(B_gpu);
    cudaFree(C_gpu);
    cudaFree(O_gpu);
}

int main(void) {
    cublasHandle_t handle;
    cublasCreate(&handle);
    test_fp16(handle);
    cublasDestroy(handle);
    return 0;
}