CUDA 半浮點數運算
cudaSamples裡面0_Simple裡面有個關於fp16的例子，做fp16向量的點積的。自己簡單實現一個，做個對自己的測試。
1、關於fp16定義
CUDA 7.5 新特性介紹 -- FP16（即fp16是cuda7.5引入的，需要計算能力達到5.3或以上）
16-bit Floating Point (半精度浮點)
從 Tegra X1 開始，NVIDIA 的 GPU 將支援原生的 FP16 計算指令，理論上可以獲得兩倍於 FP32 （單精度浮點）的效能，適用於大規模的神經網路或者計算機視覺相關的應用。而從 CUDA 7.5 開始，開發者可以很容易的將原有的 FP32 的程式碼移植為 FP16：直接使用內建的half
CUDA 7.5 主要提供以下三種 FP16 相關的功能:
-
新的 cuda_fp16.h 標頭檔案定義了 half 和 half2 型別，併為 FP32 和 FP16 之間的型別轉換提供了half2float() 與float2half() 兩個函式。
-
新的 ”cublasSgemmEx()“ 介面實現了混合精度的矩陣乘法（在輸入 FP16 的情況下以 32 位的精度進行計算）以此在保證精度的前提下處理兩倍於原有規模的矩陣運算。
-
對於現有的 Tegra X1 裝置以及未來的 GPU 型號（如下一代 Pascal 架構），CUDA 7.5 中的 cuda_fp16.h 標頭檔案提供了一系列的 intrinsics 來幫助開發者實現 高效的 FP16 計算（FP16x2 SIMD 指令）。另外 cuBLAS 也新加入了一個高度優化的 cublasHgemm() 實現，以在這類裝置上提供高效能的半精度浮點的矩陣乘法。
-
符號位: 1 bit
-
指數位: 5 bits
-
有效數位: 11 bits (10 位顯式)
半精度浮點範圍：CUDA 7.5 中定義的 half2 結構在一個32位的字中儲存了兩個半精度浮點數，如下圖所示。在 GPU 中這些 half2 的型別將以寬度為 2 的 SIMD 方式進行計算，這也是解釋了為什麼 FP16 的效能可以兩倍於 FP32。
2、關於IEEEp16標準的詳細定義
如上節最後一個連線的wiki頁面
一些half二進位制表示所對應的值
0 01111 0000000000 = 1 0 01111 0000000001 = 1 + 2−10 = 1.0009765625 (next smallest float after 1) 1 10000 0000000000 = −2 0 11110 1111111111 = 65504 (max half precision) 0 00001 0000000000 = 2−14 ≈ 6.10352 × 10−5 (minimum positive normal) 0 00000 1111111111 = 2−14 - 2−24 ≈ 6.09756 × 10−5 (maximum subnormal) 0 00000 0000000001 = 2−24 ≈ 5.96046 × 10−8 (minimum positive subnormal) 0 00000 0000000000 = 0 1 00000 0000000000 = −0 0 11111 0000000000 = infinity 1 11111 0000000000 = −infinity 0 01101 0101010101 = 0.333251953125 ≈ 1/3
3、自己寫的fp16加法
加法函式
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "cuda_fp16.h"
#include <stdio.h>
#include <iostream>
using namespace std;
#define CHECK(call) \
{ \
const cudaError_t error = call; \
if (error != cudaSuccess) \
{ \
printf("Error: %s: %d, ", __FILE__, __LINE__); \
printf("code: %d, reason: %s\n", error, cudaGetErrorString(error)); \
system("pause"); \
} \
}
__global__ void myHalf2Add(half2 *a, half2 *b, half2 *c, int size)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
c[i] = __hadd2(a[i], b[i]);
}
__global__ void float22Half2Vec(float2 * src, half2 *des, int size)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
des[i] = __float22half2_rn(src[i]);
}
__global__ void half22Float2Vec(half2 *src, float2 *des, int size)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
des[i] = __half22float2(src[i]);
}
int main()
{
const int blocks = 128;
const int threads = 128;
size_t size = blocks*threads * 2;
float *vec1 = new float[size];
float *vec2 = new float[size];
float *res = new float[size];
for (size_t i = 0; i < size; i++)
{
vec2[i] = vec1[i] = i;
}
float * vecDev1, *vecDev2, *resDev;
CHECK(cudaMalloc((void **)&vecDev1, size * sizeof(float)));
CHECK(cudaMalloc((void **)&vecDev2, size * sizeof(float)));
CHECK(cudaMalloc((void **)&resDev, size * sizeof(float)));
CHECK(cudaMemcpy(vecDev1, vec1, size * sizeof(float), cudaMemcpyHostToDevice));
CHECK(cudaMemcpy(vecDev2, vec2, size * sizeof(float), cudaMemcpyHostToDevice));
half2 *vecHalf2Dev1, *vecHalf2Dev2, *resHalf2Dev;
CHECK(cudaMalloc((void **)&vecHalf2Dev1, size * sizeof(float) / 2));
CHECK(cudaMalloc((void **)&vecHalf2Dev2, size * sizeof(float) / 2));
CHECK(cudaMalloc((void **)&resHalf2Dev, size * sizeof(float) / 2));
float22Half2Vec << <128, 128 >> > ((float2*)vecDev1, vecHalf2Dev1, size);
float22Half2Vec << <128, 128 >> > ((float2*)vecDev2, vecHalf2Dev2, size);
myHalf2Add << <128, 128 >> > (vecHalf2Dev1, vecHalf2Dev2, resHalf2Dev, size);
half22Float2Vec << <128, 128 >> >(resHalf2Dev, (float2*)resDev, size);
//half22Float2Vec << <128, 128 >> >(vecHalf2Dev1, (float2*)resDev, size);
//CHECK(cudaMemcpy(res, resDev, size * sizeof(float), cudaMemcpyDeviceToHost));
CHECK(cudaMemcpy(res, resDev, size * sizeof(float), cudaMemcpyDeviceToHost));
for (int i = 0; i < 128; i++)//打印出前64個結果，並與CPU結果對比
{
cout << vec1[i] << " + " << vec2[i] << " = " << vec1[i] + vec2[i] << " ? " << res[i] << endl;
}
for (int i = 128*128; i < 128*128+128; i++)//打印出前64個結果，並與CPU結果對比
{
cout << vec1[i] << " + " << vec2[i] << " = " << vec1[i] + vec2[i] << " ? " << res[i] << endl;
}
delete[] vec1;
delete[] vec2;
delete[] res;
CHECK(cudaFree(vecDev1));
CHECK(cudaFree(vecDev2));
CHECK(cudaFree(resDev));
CHECK(cudaFree(vecHalf2Dev1));
CHECK(cudaFree(vecHalf2Dev2));
CHECK(cudaFree(resHalf2Dev));
system("pause");
return 0;
}程式碼是非常簡單的。
half2定義和運算所需要的標頭檔案為cuda_fp16.h
以下為計算結果的一部分，由於使用位元組少，計算精度低，所以會有很多資料存在誤差。
wiki裡面有關於在各個資料範圍內，資料的誤差。其中資料位0-2048範圍內的整數時，資料是完全準確的。基本上所有資料的誤差都在千分之一量級。此處注意half型資料精度低，表示的資料範圍也比較低，最大能表示65520。
4、進一步使用所需要的文件資料
這裡是cuda的半精度內建函式和定義的相關內容。包含：
1.half和half2的算術運算
2.half和half2的比較函式
3.half和half2精度轉換和資料傳輸（包括float2在內的各種資料型別與half和half2的相互轉換）
4.half和half2的數學函式
half和half2的定義
5、遇到的問題
使用過程中一定要注意將compute_61,sm_61設定正確(需要將所有低於要求版本的選項都刪掉)，否則nvcc將預設使用該cuda版本支援的最低架構，cuda8.0將設定為20.低於所需要的5.3，因而編譯失敗。
錯誤表現為
1>F:/cuda/vsCuda/learn/fp16ScalarProductLearn/fp16ScalarProductLearn.cu(29): error : identifier "__hadd2" is undefined
6、其他應用
cublas提供了對fp16的支援，並且官方稱速度是單精度的2倍。
以下是有人在TX1上做的關於fp16+Faster R-CNN測試
也可以參考我以前的部落格實現cublasHgemm運算
很多應用都是受限於記憶體頻寬，且許多應用都會受益於低精度資料儲存時進行高精度計算。NVIDIA CEO黃建勳宣佈以後的GPU都會支援混合精度計算。
使用半精度可以在記憶體中儲存2倍大的模型，受限於記憶體頻寬的應用能獲取2倍加速（TX1也是）。
特性：
1.cuda_fp16.h中定義了half和half2相關內容
3.cublasHgemm() 提供了fp16乘法。Drive PX withTegra X1 GPUs都支援。
4.Windows Remote Desktop使用cuda7.5 windows端可以在沒有NVIDIA gpu的情況下使用遠端Windows伺服器來進行cuda開發（但沒說怎麼做）