并行编程实战—CUDA编程的精度操作

一、CUDA中的精度操作

所谓精度操作，从字面上来看，就是看数据处理的粒度，越小的粒度说明精度越高。当然，有经验的开发者迅速就会反应过来，那处理起来会越复杂，占用的空间相对也越多。CUDA中的精度操作与C++中的有些类似，它主要有以下几种：

1. 单精度类型
   FP32，类似于C++中的float，对应在CUDA中也是float，32位字长。同时，在CUDA中还支持矢量类型如float2,float4
2. 双精度类型
   FP64，类似于C++中的double,对应在CUDA中也是double，64位字长。当然也有矢量的double2等
3. 半精度类型
   FP16，C++中没有对应的数据类型，在CUDA中为__half,16位字长。

在CUDA编程中，提供了低/混合精度处理，它的目的是为了在低精度计算核心操作，再使用高精度操作生成结果的一种计算处理。举一个例子，在深度学习中，张量计算需要FP16，用来提高速度；而FP32用来保持权重数据，这样，通过不同的数据混合处理，则可以达到优化训练模型的目的。
在CUDA框架中，提供了很多内在的函数来实现相关的精度操作，比如half相关的__hfma()、__hmul()和__hadd()等。大家可以根据需要进行应用。

二、精度操作的应用

在CUDA的实际应用中，可以通过编译选项或相关的预处理指令来对精度进行控制。其主要的方式有以下几种情况：

1. 编译器选项处理

-ftz=true/false：非正规浮点数刷新为零,flush-to-zero
-prec-div=true/false：提高除法的精度
-prec-sqrt=true/false：提高平方根的精度
-fmad=true/false    乘加指令控制
-use_fast_math：使用快速但精度较低的数学函数

1. 编译指令处理
   它主要包括精度控制指令、收缩控制指令以及浮点访问控制几种，如下：：

//精度控制
#pragma float_control(precise, on/off):控制浮点精度模式
#pragma float_control(except, on/off):控制浮点异常
//收缩控制
#pragma fp_contract(on/off):控制乘加融合
//访问控制
#pragma STDC FENV_ACCESS ON/OFF:浮点环境访问

1. 标准或内在函数
   当然，在CUDA的框架中也内置上不少的相关函数，可以提供下着的精度计算处理，如“sinf、expf”用于单精度，“sin、exp”用于双精度等。这里需要说明一下，标准函数是指与C等一样的库函数如刚刚提及的sin等，它符合IEEE标准；而内在函数则是指带双下划线开头的函数，如__sin等，大家可以把它当成内部的标准函数，它可以直接映射到硬件指令。

如同在其它语言中一样，CUDA框架中也提供了精度之间转换，主要有：

1. 隐式（自动）转换
   即不需要开发者干预，自动将相关精度类型进行转换并忽略相关的精度损失，它的风险同其它语言类似，可能导致数据计算的精度下降甚至溢出导致错误，示例如下：

__half h = __float2half(1.0f);
float f = 2.2f;
double d = 6.0;

// 半精度加单精度时，半精度自动提升为单精度
float result1 = h + f;  // h自动提升为float

1. 显式转换
   利用CUDA框架中提供的接口进行转换，如：

float f1 = 3.0f;
__half hf = __float2half(f1);      // 单精度转半精度

精度的转换主要目的是为了提高效率，但需要在转换过程中小心控制风险，包括前面提到的精度损失和溢出等可能引起的异常等。所以开发者一定要善于利用内置的相关函数并了解相关的转换规则，这才是根本。

三、例程

精度操作不是一个专门的技术问题，而一个数据处理支持的问题，所以大家在应用是要视情况而定。相关的示例如下：

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


__host__ __device__ floatcaclError(floatstd, float cacl){
if (fabsf(std) < 1e-12) {
return fabsf(cacl - std);
    }
return fabsf((cacl - std) / std);
}

floatstdFunc(float x){
return sinf(x) * expf(x) / (1.0f + logf(1.0f + x*x));
}
//high
__global__ voidcaclHighPrec(constfloat* input, float* output, int n){
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx >= n) return;

float x = input[idx];
    output[idx] = sinf(x) * expf(x) / (1.0f + logf(1.0f + x*x));
}
//fast
__global__ voidcaclFastPrec(constfloat* input, float* output, int n){
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx >= n) return;

float x = input[idx];
    output[idx] = __sinf(x) * __expf(x) / (1.0f + __logf(1.0f + x*x));
}
//half
__global__ voidcaclHalfPrec(constfloat* input, __half* output, int n){
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx >= n) return;

float x = input[idx];
    __half hx = __float2half_rn(x);
    __half hOne = __float2half_rn(1.0f);

    __half hSin = hsin(hx);
    __half hExp = hexp(hx);
    __half hSqr = __hmul(hx, hx);
    __half hLog = hlog(__hadd(hOne, hSqr));
    __half hDenom = __hadd(hOne, hLog);
    __half hNum = __hmul(hSin, hExp);
    output[idx] = __hdiv(hNum, hDenom);
}

// mix
__global__ voidcaclPrecKernel(constfloat* input, float* output, int n){
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx >= n) return;

float sum = 0.0f;
float c = 0.0f;

for (int i = 0; i < 100; i++) {
float x = input[idx] * 0.01f * i;
float y = __sinf(x) - c;    
float t = sum + y;
        c = (t - sum) - y;       
        sum = t;
    }

    output[idx] = sum;
}

intmain(){
constint N = 1024 * 1024;
constsize_t size = N * sizeof(float);
constsize_t halfSize = N * sizeof(__half);


float* hInput = (float*)malloc(size);
float* hOutputH = (float*)malloc(size);
float* hOutputF = (float*)malloc(size);
float* hOutputM = (float*)malloc(size);
    __half* hOutputHalf = (__half*)malloc(halfSize);

for (int i = 0; i < N; i++) {
        hInput[i] = (float)i / 1000.0f;  // 0 到 ~1.024
    }


float* hStand = (float*)malloc(size);
for (int i = 0; i < N; i++) {
        hStand[i] = stdFunc(hInput[i]);
    }

float *dInput, *dOutputH, *dOutputF, *dOutputM;
    __half *dOutputHalf;

     cudaMalloc(&dInput, size);
     cudaMalloc(&dOutputH, size);
     cudaMalloc(&dOutputF, size);
     cudaMalloc(&dOutputM, size);
     cudaMalloc(&dOutputHalf, halfSize);

     cudaMemcpy(dInput, hInput, size, cudaMemcpyHostToDevice);

int threads = 256;
int blocks = (N + threads - 1) / threads;


// high
    caclHighPrec << <blocks, threads >> >(dInput, dOutputH, N);
//  fast
    caclFastPrec << <blocks, threads >> >(dInput, dOutputF, N);
// half
    caclHalfPrec << <blocks, threads >> >(dInput, dOutputHalf, N);
//mix
    caclPrecKernel << <blocks, threads >> >(dInput, dOutputM, N);

     cudaMemcpy(hOutputH, dOutputH, size, cudaMemcpyDeviceToHost);
     cudaMemcpy(hOutputF, dOutputF, size, cudaMemcpyDeviceToHost);
     cudaMemcpy(hOutputM, dOutputM, size, cudaMemcpyDeviceToHost);
     cudaMemcpy(hOutputHalf, dOutputHalf, halfSize, cudaMemcpyDeviceToHost);

double totalErrH = 0.0;
double totalErrF = 0.0;
double totalErrHalf = 0.0;
double mixErrH = 0.0;
double mixErrF = 0.0;
double mixErrHalf = 0.0;

for (int i = 0; i < N; i++) {
float err_high = caclError(hStand[i], hOutputH[i]);
        totalErrH += err_high;
if (err_high > mixErrH) mixErrH = err_high;

float err_fast = caclError(hStand[i], hOutputF[i]);
        totalErrF += err_fast;
if (err_fast > mixErrF) mixErrF = err_fast;

float val_half = __half2float(hOutputHalf[i]);
float err_half = caclError(hStand[i], val_half);
        totalErrHalf += err_half;
if (err_half > mixErrHalf) mixErrHalf = err_half;
    }

// clear
     cudaFree(dInput);
     cudaFree(dOutputH);
     cudaFree(dOutputF);
     cudaFree(dOutputM);
     cudaFree(dOutputHalf);

free(hInput);
free(hOutputH);
free(hOutputF);
free(hOutputHalf);
free(hOutputM);
free(hStand);

return0;
}

代码并不复杂，只要有一些其它语言的编程基础再加上CUDA相关的API说明看看，应该没有什么问题。

四、总结

其实在学习编程时，大家会不会有一种想法，在学习的知识越宽时，发现很多的技术或知识点是相通的。要么有联系，要么有改进，要么有扩展。总让人有一种似曾相识的感觉。其实这就对了，作为编程的底层，其基础和本质并没有脱离软硬件行业的底层建筑，所以就不大可能有完全革命性的颠覆。这也是和武侠小说中，一通百通大抵类似。