【CUDA】（七）调整指令级原语

（七）调整指令级原语

本篇笔记参考如下：

https://blog.csdn.net/gao_zhennan/article/details/120717424?ops_request_misc=elastic_search_misc&request_id=8a12ee6315cc5ba5c95b4dc10ccb3baf&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-120717424-null-null.142

本章重点介绍计算密集型应用。处理器的计算吞吐量可以用它在一段时间内执行操作的数量来衡量。因为GPU有很多SIMT指令和计算核心，所以其峰值计算吞吐量通常比其他的处理器要高。但并不是所有的指令都是平等的。如果结果不能正确收敛或没有获得预期的结果，那么程序运行速度再快也没有用。对应用程序的吞吐量和正确性进行优化时，理解不同低级原语的性能、数值精确度和线程安全性方面的优缺点是很重要的。

对于以下代码段

double value = in1 * in2 + in3;

这种乘法后紧跟加法的算术模式被称为乘法加，或者MAD，其应用非常广泛。任何利用向量和矩阵的应用程序都可能包含许多MAD指

令来执行线性代数中的点积运算、矩阵乘法运算，以及其他运算。一个简单的编译器会把一个MAD指令转化成两个算术指令：先进行乘法运算紧接着进行加法运算。因为这种模式很常见，所以现代运算结构（包括NVIDIA GPU）都支持MAD指令。因此，执行一个MAD的结果是循环次数减少了一半。

7.1 CUDA指令概述

指令是处理器中的一个逻辑单元。对两个功能等效指令的选择可以影响很多应用程序的特性，包括性能、精确度和正确性。当通过严格的数字验证请求，把遗留应用程序传输到CUDA时，就要特别留意这些问题了。

本节涵盖了显著影响CUDA内核生成指令的3大因素：浮点运算、内置和标准函数、原子操作。浮点运算是针对于非整数值的计算，并且会影响CUDA程序的精确度和性能。

7.1.1 浮点指令

自从浮点运算采用IEEE——754标准后，所有的主流处理器厂商都使用这一标准，包括NVIDIA。这个标准规定将二进制浮点数据编码成3段：符号段（sign），一个比特位；指数段（exponent），多个比特位；以及尾数或分数段（fraction），多个比特位。

IEEE754提供了四种精度规范, 其中最常用的是 单精度浮点型 和 双精度浮点型 , 但IEEE754并没有规定32位浮点数类型需要叫做 float, 或64位浮点数需要叫做 double. 它只是提供了一些关于如何存储不同精度浮点数的规范和标准. 不过一般情况下, 如果我们提到 float, 其实指的就是IEEE754标准中的32位单精度浮点数. 如果我们提到 double, 其实指的就是IEEE754标准中的64位双精度浮点数

32位浮点数内存占用示意图如下：

sign: 符号位, 即图中蓝色的方块

biased exponent: 偏移后的指数位, 即图中绿色的方块

fraction: 尾数位, 即图中红色的方块

1. 符号位: sign

符号位: 占据最高位(第31位)这一位, 用于表示这个浮点数是正数还是负数, 为0表示正数, 为1表示负数.

2. 偏移后的指数位: biased exponent

指数位占据第30位到第23位这8位. 也就是上图的绿色部分.

用于表示以2位底的指数. 至于这个指数的作用, 后文会详细讲解, 这里只需要知道: 8位二进制可以表示256种状态, IEEE754规定, 指数位用于表示[-127, 128]范围内的指数.

为了表示起来更方便, 浮点型的指数位都有一个固定的偏移量(bias), 用于使 指数 + 这个偏移量 = 一个非负整数.

规定: 在32位单精度类型中, 这个偏移量是127. 在64位双精度类型中, 偏移量是1023. 所以, 这里的偏移量是127

3. 尾数位:fraction

尾数位: 占据剩余的22位到0位这23位. 用于存储尾数.在以二进制格式存储十进制浮点数时, 首先需要把十进制浮点数表示为二进制格式。

举例：

拿十进制数20.5

十进制浮点数20.5 = 二进制10100.1

我们需要把这个二进制数转换为以2为底的指数形式:二进制10100.1 = 1.01001 * 2^4

转换时, 对于乘号左边, 加粗的那个二进制数1.01001, 需要把小数点放在左起第一位和第二位之间. 且第一位需要是个非0数. 这样表示好之后, 其中的1.01001就是尾数.

用 二进制数 表示 十进制浮点数 时, 表示为 尾数*指数的形式, 并把尾数的小数点放在第一位和第二位之间, 然后保证第一位数非0, 这个处理过程叫做 规范化(normalized)

规范化之后的这个数: 1.01001 * 2^4

其中1.01001是尾数, 而4就是偏移前的指数(unbiased exponent), 上文讲过, 32位单精度浮点数的偏移量(bias)为127, 所以这里加上偏移量之后, 得到的偏移后指数(biased exponent)就是 4 + 127 = 131, 131转换为二进制就是1000 0011

现在还需要对尾数做一些特殊处理：

1）隐藏高位1

尾数部分的最高位始终为1. 比如这里的 **1.**01001, 这是因为前面说过, 规范化之后, 尾数中的小数点会位于左起第一位和第二位之间. 且第一位是个非0数. 而二进制中, 每一位可取值只有0或1, 如果第一位非0, 则第一位只能为1. 所以在存储尾数时, 可以省略前面的 1和小数点. 只记录尾数中小数点之后的部分, 这样就节约了一位内存. 所以这里只需记录剩余的尾数部分: 01001

2）低位补0

有时候尾数会不够填满尾数位(即图中的红色格子). 比如这里的, 尾数01001不够23位，此时, 需要在低位补零, 补齐23位.

因此我们可以得到：

十进制浮点数 20.5 的:

符号位是: 0

偏移后指数位是: 1000 0011

补零后尾数位是: 0100 1000 0000 0000 000

现在, 把这三部分按顺序放在32位浮点数容器中, 就是 0 1000 0011 0100 1000 0000 0000 000

编写代码test.c如下

#include <stdio.h>
#include <stdlib.h>

int main()
{
    float a = 3.1415927f;
    float b = 3.1415928f;
    if(a == b)
    {
        printf("a and b are equal\n");
    }
    else
    {
        printf("a and b are not equal\n");
    }   
}

得到输出如下

a and b are equal

由上我们可以看出对于机器来说数值的精确度是有限的，并且用浮点类型存储的数据是离散且有限的。

浮点编程中需要考虑的另一个方面是浮点数的粒度问题。像上面所讨论的，浮点数的粒度比整数来说要好。然而浮点数只能在离散的区间间隔内存储数据。随着浮点数值离零越来越远（在正负两个方向上），表示数值的区间也会随之增大。

值得注意的是，随着x值的增大，精度会大幅降低。浮点数之间的区间间隔意味着在任何可能产生极端数值的应用中，对数值进行四舍五入会对输出有很大的影响。

CUDA和其他遵循IEEE——754标准的编程模式支持两种浮点精度：32位和64位。这些不同的格式也分别被称为单精度和双精度。因为双精度浮点数的位数是单精度浮点数的两倍，所以双精度可以表示更多的数值。这意味着双精度浮点数既有更好的细粒度又有比单精度值更大的数值范围。

修改代码如下：

#include <stdio.h>
#include <stdlib.h>

int main()
{
    double a = 3.1415927;
    double b = 3.1415928;
    if(a == b)
    {
        printf("a and b are equal\n");
    }
    else
    {
        printf("a and b are not equal\n");
    }   
}

可以得到以下输出

a and b are not equal

当用双精度变量存储时，a和b最近似的表示值是不同的。

7.1.2 内部函数和标准函数

除了单精度和双精度操作的区别，CUDA还将所有算数函数分成内部函数和标准函数。标准函数用于支持可对主机和设备进行访问并标准化主机和设备的操作。标准函数包含来自于C标准数学库的数学运算，如sqrt、exp和sin。

CUDA内置函数只能对设备代码进行访问。在编程中，如果一个函数是内部函数或是内置函数，那么在编译时对它的行为会有特殊响应，从而产生更积极的优化和更专业化的指令生成。事实上，许多三角函数是直接在GPU硬件上实现的，因为它们中的大部分是用图形应用计算的（变换、旋转和其他在3D可视化应用上的操作）。

7.1.3 原子操作指令

一条原子操作指令用来执行一个数学运算，此操作是一个独立不间断的操作，且没有其他线程的干扰。通常一个“加法”操作在硬件层面分为三步：读取旧值、计算新值、写回新值。在并发环境下，如果两个线程同时读取了同一个旧值，它们的写回操作会互相覆盖，导致结果错误。 原子指令 会将这三步封装成一个不可分割的整体，确保在当前线程完成“读-改-写”之前，其他线程无法干扰该内存地址。

CUDA 提供了丰富的原子函数，主要分为以下几类：

A. 算术运算 (Arithmetic Functions)

• atomicAdd(): 原子加法（支持整型和浮点型）。
• atomicSub(): 原子减法（主要支持整型）。
• atomicExch(): 原子交换，将新值写入并将旧值返回。
• atomicMin() / atomicMax(): 原子最小值/最大值。
• atomicInc() / atomicDec(): 原子自增/自减（通常带有一个边界值，循环计数时常用）。

B. 位运算 (Bitwise Functions)

主要用于操作整型数据的位标志：

• atomicAnd(): 按位与。
• atomicOr(): 按位或。
• atomicXor(): 按位异或。

C. 比较与交换 (Compare-and-Swap)

• atomicCAS(): 这是最基础、最强大的原子指令。
  • 语法: old = atomicCAS(address, compare, val);
  • 逻辑: 如果 address 处的值等于 compare，则将 val 写入；无论是否写入，都返回旧值。
  • 用途: 可以用来实现自定义的复杂原子操作（如双精度浮点数的原子乘法）。

支持的内存范围

原子操作可以在以下两种内存空间执行：

• 全局内存 (Global Memory)：所有线程块都可以参与，常用于全局结果汇总。
• 共享内存 (Shared Memory)：仅限于同一个线程块内的线程，速度比全局内存快得多，常用于线程块内的局部规约（Reduction）或直方图统计。

性能影响：序列化 (Serialization)

虽然原子操作保证了正确性，但它是有代价的：

• 冲突与排队：当大量线程竞争同一个地址时，操作会被序列化（即一个接一个执行）。这会极大地降低并行效率。
• 优化建议：尽量减少对同一个地址的竞争。例如，可以先在共享内存中使用原子操作进行“局部汇总”，最后再由一个线程将结果原子加到全局内存中。

7.2 程序优化指令

用于优化程序的指令，有很多的选择：单精度或双精度浮点值、标准或内部函数、原子函数或不安全访问。

7.2.1 单精度与双精度的比较

双精度变量相较于单精度变量来说，可以在一个更精细的粒度和更广泛的范围上表示不同的数值。以代码floating-point-accuracy.cu为例

/*部分代码*/
__global__ void kernel(float *F, double *D)
{
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid == 0)
    {
        *F = 12.1;
        *D = 12.1;
    }
}

    float *deviceF;
    float h_deviceF;
    double *deviceD;
    double h_deviceD;

    float hostF = 12.1;
    double hostD = 12.1;

    CHECK(cudaMalloc((void **)&deviceF, sizeof(float)));
    CHECK(cudaMalloc((void **)&deviceD, sizeof(double)));
    kernel<<<1, 32>>>(deviceF, deviceD);
    CHECK(cudaMemcpy(&h_deviceF, deviceF, sizeof(float),
                     cudaMemcpyDeviceToHost));
    CHECK(cudaMemcpy(&h_deviceD, deviceD, sizeof(double),
                     cudaMemcpyDeviceToHost));

}

可以得到输出如下

Host single-precision representation of 12.1   = 12.10000038146972656250
Host double-precision representation of 12.1   = 12.09999999999999964473
Device single-precision representation of 12.1 = 12.10000038146972656250
Device double-precision representation of 12.1 = 12.09999999999999964473
Device and host single-precision representation equal? yes
Device and host double-precision representation equal? yes

在这个特殊的例子中，双精度数值比单精度数值稍微更接近于真实数值。

双精度数值的精确性是以空间和性能消耗为代价的。以代码floating-point-perf.cu为例

static void run_float_test(size_t N, int niters, int blocksPerGrid,
                           int threadsPerBlock, double *toDeviceTime,
                           double *kernelTime, double *fromDeviceTime,
                           float *sample, int sampleLength)
{
    int i;
    float *h_floatInputs, *h_floatOutputs;
    float *d_floatInputs, *d_floatOutputs;

    h_floatInputs = (float *)malloc(sizeof(float) * N);
    h_floatOutputs = (float *)malloc(sizeof(float) * N);
    CHECK(cudaMalloc((void **)&d_floatInputs, sizeof(float) * N));
    CHECK(cudaMalloc((void **)&d_floatOutputs, sizeof(float) * N));

    for (i = 0; i < N; i++)
    {
        h_floatInputs[i] = (float)i;
    }

    double toDeviceStart = seconds();
    CHECK(cudaMemcpy(d_floatInputs, h_floatInputs, sizeof(float) * N,
                     cudaMemcpyHostToDevice));
    *toDeviceTime = seconds() - toDeviceStart;

    double kernelStart = seconds();
    lots_of_float_compute<<<blocksPerGrid, threadsPerBlock>>>(d_floatInputs,
            N, niters, d_floatOutputs);
    CHECK(cudaDeviceSynchronize());
    *kernelTime = seconds() - kernelStart;

    double fromDeviceStart = seconds();
    CHECK(cudaMemcpy(h_floatOutputs, d_floatOutputs, sizeof(float) * N,
                     cudaMemcpyDeviceToHost));
    *fromDeviceTime = seconds() - fromDeviceStart;

    for (i = 0; i < sampleLength; i++)
    {
        sample[i] = h_floatOutputs[i];
    }

    CHECK(cudaFree(d_floatInputs));
    CHECK(cudaFree(d_floatOutputs));
    free(h_floatInputs);
    free(h_floatOutputs);
}

/**
 * Runs a full test of double-precision floating-point, including transferring
 * inputs to the device, running the single-precision kernel, and copying
 * outputs back.
 **/
static void run_double_test(size_t N, int niters, int blocksPerGrid,
                            int threadsPerBlock, double *toDeviceTime,
                            double *kernelTime, double *fromDeviceTime,
                            double *sample, int sampleLength)
{
    int i;
    double *h_doubleInputs, *h_doubleOutputs;
    double *d_doubleInputs, *d_doubleOutputs;

    h_doubleInputs = (double *)malloc(sizeof(double) * N);
    h_doubleOutputs = (double *)malloc(sizeof(double) * N);
    CHECK(cudaMalloc((void **)&d_doubleInputs, sizeof(double) * N));
    CHECK(cudaMalloc((void **)&d_doubleOutputs, sizeof(double) * N));

    for (i = 0; i < N; i++)
    {
        h_doubleInputs[i] = (double)i;
    }

    double toDeviceStart = seconds();
    CHECK(cudaMemcpy(d_doubleInputs, h_doubleInputs, sizeof(double) * N,
                     cudaMemcpyHostToDevice));
    *toDeviceTime = seconds() - toDeviceStart;

    double kernelStart = seconds();
    lots_of_double_compute<<<blocksPerGrid, threadsPerBlock>>>(d_doubleInputs,
            N, niters, d_doubleOutputs);
    CHECK(cudaDeviceSynchronize());
    *kernelTime = seconds() - kernelStart;

    double fromDeviceStart = seconds();
    CHECK(cudaMemcpy(h_doubleOutputs, d_doubleOutputs, sizeof(double) * N,
                     cudaMemcpyDeviceToHost));
    *fromDeviceTime = seconds() - fromDeviceStart;

    for (i = 0; i < sampleLength; i++)
    {
        sample[i] = h_doubleOutputs[i];
    }

    CHECK(cudaFree(d_doubleInputs));
    CHECK(cudaFree(d_doubleOutputs));
    free(h_doubleInputs);
    free(h_doubleOutputs);
}

编译代码并运行

nvcc -arch=sm_120 floating-point-perf.cu -o perf
./perf

得到结果如下

Running 3472128 blocks with 256 threads/block over 888864768 elements
Input   Diff Between Single- and Double-Precision
------  ------
0       1.16110611328622326255e-01
1       1.42341757498797960579e-01
2       1.45135404032771475613e-01
3       1.47929050144739449024e-01
4       1.03847696445882320404e-01
5       1.84766342732473276556e-01
6       1.48497488888096995652e-01
7       1.20041135203791782260e-01
8       1.38459781592246145010e-01
9       1.49065927878837101161e-01

For single-precision floating point, mean times for:
  Copy to device:   0.463959 s
  Kernel execution: 0.135360 s
  Copy from device: 1.732521 s
For double-precision floating point, mean times for:
  Copy to device:   0.848217 s (1.83x slower than single-precision)
  Kernel execution: 9.921548 s (73.30x slower than single-precision)
  Copy from device: 5.593063 s (3.23x slower than single-precision)

数据显示了单精度与双精度结果之间的差值（Diff），数量级大约在 $10^{-1}$。

• 误差积累：对于 $8.8$ 亿个元素的大规模运算，单精度的累积误差变得非常明显。
• 离散区间问题：浮点数只能在离散区间内存储数据，且数值越大，表示数值的区间（粒度）也随之增大。在处理近 9 亿个元素时，单精度由于尾数位不足（仅 23 位），在频繁加法或复杂运算中会产生严重的舍入误差。

双精度内核执行比单精度慢了 73.30 倍（9.92s vs 0.135s）。

• 硬件吞吐量限制：大多数消费级 GPU（如我的 RTX 5060 Ti）的 FP32（单精度）单元数量远多于 FP64（双精度）单元。通常 FP64 的吞吐量仅为 FP32 的 $1/32$ 甚至更低。
• 寄存器压力：双精度变量占用 64 位，是单精度的两倍。这会导致核函数运行过程中寄存器使用量翻倍，从而降低占用率（Occupancy），限制了网格级并发的能力。
• 指令带宽：双精度指令本身的发射和执行周期在硬件层面就比单精度长。

数据搬运的耗时也呈现出明显的倍数关系：

• H2D (Copy to device): 双精度慢了 1.83x。这非常符合逻辑，因为双精度数据量恰好是单精度的 2 倍，传输耗时理论上应翻倍。
• D2H (Copy from device): 双精度慢了 3.23x。这里的耗时增加超过了 2 倍，可能的原因是双精度在写回 CPU 时，PCIe 总线的带宽利用率或主机端的内存写入开销更高。

一般情况下，如果应用程序对精确度要求很高的话，那么必须使用双精度数值。否则，使用单精度数值可以获得性能提升。

7.2.2 标准函数与内部函数的比较

标准函数和内部函数在数值精确度和性能上的表现是不同的。标准函数支持大部分的数学运算。但是，许多等效的内部函数能够使用较少的指令、改进的性能和更低的数值精确度，实现相同的功能。

在foo.cu上编写以下代码

#include "../common/common.h"
#include <stdio.h>
#include <stdlib.h>

__global__ void intrinsic(float *ptr)
{
    *ptr = __powf(*ptr, 2.0f);
}

__global__ void standard(float *ptr)
{
    *ptr = powf(*ptr, 2.0f);
}

输入以下命令生成ptx文件：

nvcc --ptx -o foo.ptx foo.cu

可以看到

对于内部函数__powf只需要十几行代码

.visible .entry _Z9intrinsicPf(
	.param .u64 _Z9intrinsicPf_param_0
)
{
	.reg .f32 	%f<5>;
	.reg .b64 	%rd<3>;


	ld.param.u64 	%rd1, [_Z9intrinsicPf_param_0];
	cvta.to.global.u64 	%rd2, %rd1;
	ld.global.f32 	%f1, [%rd2];
	lg2.approx.f32 	%f2, %f1;
	add.f32 	%f3, %f2, %f2;
	ex2.approx.f32 	%f4, %f3;
	st.global.f32 	[%rd2], %f4;
	ret;

}

而powf函数则需要一百多行：

.visible .entry _Z8standardPf(
	.param .u64 _Z8standardPf_param_0
)
{
	.reg .pred 	%p<13>;
	.reg .f32 	%f<84>;
	.reg .b32 	%r<13>;
	.reg .b64 	%rd<3>;


	ld.param.u64 	%rd2, [_Z8standardPf_param_0];
	cvta.to.global.u64 	%rd1, %rd2;
	mov.f32 	%f83, 0f3F800000;
	cvt.rzi.f32.f32 	%f10, %f83;
	add.f32 	%f11, %f10, %f10;
	mov.f32 	%f12, 0f40000000;
	sub.f32 	%f13, %f12, %f11;
	abs.f32 	%f1, %f13;
	ld.global.f32 	%f2, [%rd1];
	abs.f32 	%f3, %f2;
	setp.lt.f32 	%p1, %f3, 0f00800000;
	mul.f32 	%f14, %f3, 0f4B800000;
	selp.f32 	%f15, %f14, %f3, %p1;
	selp.f32 	%f16, 0fC1C00000, 0f00000000, %p1;
	mov.b32 	%r1, %f15;
	add.s32 	%r2, %r1, -1060439283;
	and.b32  	%r3, %r2, -8388608;
	sub.s32 	%r4, %r1, %r3;
	mov.b32 	%f17, %r4;
	cvt.rn.f32.s32 	%f18, %r3;
	mov.f32 	%f19, 0f34000000;
	fma.rn.f32 	%f20, %f18, %f19, %f16;
	add.f32 	%f21, %f17, 0fBF800000;
	add.f32 	%f22, %f17, 0f3F800000;
	rcp.approx.ftz.f32 	%f23, %f22;
	add.f32 	%f24, %f21, %f21;
	mul.f32 	%f25, %f24, %f23;
	mul.f32 	%f26, %f25, %f25;
	sub.f32 	%f27, %f21, %f25;
	add.f32 	%f28, %f27, %f27;
	neg.f32 	%f29, %f25;
	fma.rn.f32 	%f30, %f29, %f21, %f28;
	mul.rn.f32 	%f31, %f23, %f30;
	mov.f32 	%f32, 0f3B52E7DB;
	mov.f32 	%f33, 0f3A2C32E4;
	fma.rn.f32 	%f34, %f33, %f26, %f32;
	mov.f32 	%f35, 0f3C93BB73;
	fma.rn.f32 	%f36, %f34, %f26, %f35;
	mov.f32 	%f37, 0f3DF6384F;
	fma.rn.f32 	%f38, %f36, %f26, %f37;
	mul.rn.f32 	%f39, %f38, %f26;
	mov.f32 	%f40, 0f3FB8AA3B;
	fma.rn.f32 	%f41, %f25, %f40, %f20;
	sub.f32 	%f42, %f20, %f41;
	fma.rn.f32 	%f43, %f25, %f40, %f42;
	fma.rn.f32 	%f44, %f31, %f40, %f43;
	mov.f32 	%f45, 0f32A55E34;
	fma.rn.f32 	%f46, %f25, %f45, %f44;
	mul.f32 	%f47, %f39, 0f40400000;
	fma.rn.f32 	%f48, %f47, %f31, %f46;
	fma.rn.f32 	%f49, %f39, %f25, %f48;
	add.rn.f32 	%f50, %f41, %f49;
	neg.f32 	%f51, %f41;
	add.rn.f32 	%f52, %f50, %f51;
	neg.f32 	%f53, %f52;
	add.rn.f32 	%f54, %f49, %f53;
	mul.rn.f32 	%f55, %f50, %f12;
	neg.f32 	%f56, %f55;
	fma.rn.f32 	%f57, %f50, %f12, %f56;
	fma.rn.f32 	%f58, %f54, %f12, %f57;
	cvt.rni.f32.f32 	%f59, %f55;
	sub.f32 	%f60, %f55, %f59;
	add.f32 	%f61, %f58, %f60;
	mov.f32 	%f62, 0f3AAF85ED;
	mov.f32 	%f63, 0f391FCB8E;
	fma.rn.f32 	%f64, %f63, %f61, %f62;
	mov.f32 	%f65, 0f3C1D9856;
	fma.rn.f32 	%f66, %f64, %f61, %f65;
	mov.f32 	%f67, 0f3D6357BB;
	fma.rn.f32 	%f68, %f66, %f61, %f67;
	mov.f32 	%f69, 0f3E75FDEC;
	fma.rn.f32 	%f70, %f68, %f61, %f69;
	mov.f32 	%f71, 0f3F317218;
	fma.rn.f32 	%f72, %f70, %f61, %f71;
	fma.rn.f32 	%f73, %f72, %f61, %f83;
	cvt.rzi.s32.f32 	%r5, %f59;
	setp.gt.f32 	%p2, %f59, 0f00000000;
	selp.b32 	%r6, 0, -2097152000, %p2;
	add.s32 	%r7, %r6, 2130706432;
	mov.b32 	%f74, %r7;
	mul.f32 	%f75, %f73, %f74;
	shl.b32 	%r8, %r5, 23;
	sub.s32 	%r9, %r8, %r6;
	mov.b32 	%f76, %r9;
	mul.f32 	%f77, %f75, %f76;
	abs.f32 	%f78, %f55;
	setp.gt.f32 	%p3, %f78, 0f43180000;
	setp.lt.f32 	%p4, %f55, 0f00000000;
	selp.f32 	%f79, 0f00000000, 0f7F800000, %p4;
	selp.f32 	%f4, %f79, %f77, %p3;
	setp.eq.f32 	%p5, %f2, 0f3F800000;
	@%p5 bra 	$L__BB1_7;

	setp.gtu.f32 	%p6, %f3, 0f7F800000;
	@%p6 bra 	$L__BB1_6;
	bra.uni 	$L__BB1_2;

$L__BB1_6:
	mov.f32 	%f82, 0f40000000;
	add.rn.f32 	%f83, %f2, %f82;
	bra.uni 	$L__BB1_7;

$L__BB1_2:
	setp.eq.f32 	%p7, %f2, 0f00000000;
	setp.eq.f32 	%p8, %f3, 0f7F800000;
	or.pred  	%p9, %p7, %p8;
	@%p9 bra 	$L__BB1_5;
	bra.uni 	$L__BB1_3;

$L__BB1_5:
	setp.eq.f32 	%p12, %f1, 0f3F800000;
	add.f32 	%f81, %f2, %f2;
	mov.b32 	%r10, %f81;
	and.b32  	%r11, %r10, 2147483647;
	selp.b32 	%r12, %r10, %r11, %p12;
	mov.b32 	%f83, %r12;
	bra.uni 	$L__BB1_7;

$L__BB1_3:
	setp.geu.f32 	%p10, %f2, 0f00000000;
	mov.f32 	%f83, %f4;
	@%p10 bra 	$L__BB1_7;

	setp.eq.f32 	%p11, %f1, 0f3F800000;
	neg.f32 	%f80, %f4;
	selp.f32 	%f83, %f80, %f4, %p11;

$L__BB1_7:
	st.global.f32 	[%rd1], %f83;
	ret;

}

区分标准函数和内部函数的不仅有性能，它们的计算精度也是不同的。为了测试性能和精确度的不同。我们采用代码intrinsic-standard-comp.cu

在该程序中的核函数中，先使用标准函数powf，再使用内部函数__powf,利用它们反复计算输入值的平方根。这个例子也使用主机上的C标准数学库来执行相同的计算，并使用主机上的结果作为基准值。

编译并运行

nvcc -arch=sm_120 intrinsic-standard-comp.cu -o intrinsic-standard-comp
./intrinsic-standard-comp

得到以下结果

Host calculated                 66932852.000000
Standard Device calculated      66932852.000000
Intrinsic Device calculated     66932804.000000
Host equals Standard?           Yes diff=0.000000e+00
Host equals Intrinsic?          No diff=4.800000e+01
Standard equals Intrinsic?      No diff=4.800000e+01

Mean execution time for standard function powf:    0.001979 s
Mean execution time for intrinsic function __powf: 0.000953 s

从输出结果看：

• Host (CPU): 66932852.0
• Standard (powf): 66932852.0（与 Host 完全一致）
• Intrinsic (__powf): 66932804.0（误差为 48.0）

Standard Functions (powf)：遵循 IEEE 754 标准，力求最高的数值精度。在 GPU 上，这些函数通常包含更多的指令来处理边界情况和舍入，以确保结果与 CPU 端的计算高度一致。

Intrinsic Functions (__powf)：这是 CUDA 提供的“硬件加速”版本。它们直接映射到 GPU 的 SFU（Special Function Units，特殊函数单元）。

• 为了极致的吞吐量，内建函数简化了计算过程。
• 根据 NVIDIA 文档，__powf(x, y) 实际上是通过 $2^{y \cdot \log_2(x)}$ 这种方式快速计算的。

在大多数情况下，将程序员编写的内核代码转换为GPU指令集这一过程是由CUDA编译器完成的。程序员很少会有检查或手动修改指令的想法。CUDA编译器中有两种方法可以控制指令级优化类型：编译器标志、内部或标准函数调用。

例如：内部函数__fdividef与运算符“/”相比，在执行浮点数除法时速度更快但数值精确度相对较低。

nvcc的–fmad选项可全局性地启用或禁用FMAD整个编译单元的优化。默认情况下，nvcc使用“–fmad=true”以启用FMAD指令来优化性能。“–fmad=false”的意思是阻止编译器混合任何乘法和加法

核心标志：-use_fast_math：将所有标准函数（如 sinf(), cosf(), powf(),expf()）自动替换为对应的内建函数（如 __sinf(), __powf()）。

或者可以根据需要单独使用以下子标志。-use_fast_math 实际上就是以下标志的组合：

标志	说明	影响
-ftz=true	Flush-to-Zero	将极小的非规范化浮点数视为 0，提升处理速度，但会丢失极小值的精度。
-prec-div=false	Faster Division	使用精度较低但速度更快的除法算法。
-prec-sqrt=false	Faster Sqrt	使用精度较低但速度更快的平方根算法。
-fmad=true	Fused Multiply-Add	默认开启。将 $x \cdot y + z$ 合并为一条指令，既提升速度通常还能提升精度。

7.2.3 了解原子指令

7.2.3.1 从头开始

通过使用一个原子函数，每个由CUDA提供的原子函数可以重复被执行：原子级比较并交换（CAS）运算符。

CAS将3个内容作为输入：内存地址、存储在此地址中的期望值，

以及实际想要存储在此位置的新值，然后执行以下几步：

1.读取目标地址并将该处地址的存储值与预期值进行比较。

a.如果存储值与预期值相等，那么新值将存入目标位置。

b.如果存储值与预期值不等，那么目标位置不会发生变化。

2.不论发生什么情况，一个CAS操作总是返回目标地址中的值。注意，使用返回值可以用来检查一个数值是否被替换成功。

场景模拟：假设内存地址 0x100 里的值是 A。

1. 线程 1 读到了 A，想把它改成 B。
2. 在线程 1 计算出 B 的过程中，线程 2 突然插一脚，把地址里的值改成了 C。
3. 如果线程 1 直接覆盖写入，那么线程 2 的修改就丢失了。

CAS 的逻辑是： 线程 1 说：“我要把 0x100 改成 B，但我有个前提：这个位置现在必须还是 A。如果被别人偷偷改成了 C，那我就不改了，你把现在的真实值告诉我，我重新考虑。”

步骤 1：比较与写入（原子性保证）

这是 CAS 的核心。在硬件层面，这一步是不可中断的。

• 命中 (Success)：存储值 == 预期值。硬件认为：“很好，自你上次读取以来，没人动过这个变量。” 于是放心地把新值存进去。
• 错失 (Fail)：存储值 != 预期值。硬件认为：“有人捷足先登了！” 为了保护数据一致性，它拒绝执行写入。

步骤 2：总是返回“旧值”

无论更新成功还是失败，CAS 都会把当前内存里的实际值吐出来给你。

• 成功时：返回的值等于你的“预期值”。
• 失败时：返回的值不等于你的“预期值”。这时你拿到的就是“别人修改后的新值”。

CUDA 官方只提供了 atomicAdd（加法）等基础操作。如果你想做一个原子乘法（CUDA 原生不支持），就必须用 CAS。

7.2.3.2 CUDA原子函数

1. 原子算术运算 (Atomic Arithmetic)

这类函数最常用于统计、累加和直方图计算。

函数名	支持的类型	描述
atomicAdd()	int, unsigned int, unsigned long long, float, double	将目标地址的值增加 val
atomicSub()	int, unsigned int	将目标地址的值减少 val
atomicExch()	int, unsigned int, unsigned long long, float	写入新值，并返回旧值
atomicMin()	int, unsigned int, long long, unsigned long long	比较并保留较小的值
atomicMax()	int, unsigned int, long long, unsigned long long	比较并保留较大的值
atomicInc()	unsigned int	原子自增，超过阈值则重置为 0 (用于循环计数)
atomicDec()	unsigned int	原子自减，低于 0 或超过阈值则设为阈值

2. 原子位运算 (Atomic Bitwise)

主要用于标志位设置、掩码操作或位图计算。这些操作仅支持整数类型。

函数名	支持的类型	逻辑操作
atomicAnd()	int, unsigned int, unsigned long long	*address &= val
atomicOr()	int, unsigned int, unsigned long long	`*address
atomicXor()	int, unsigned int, unsigned long long	*address ^= val

3. 比较并交换 (Atomic CAS)

这是最底层的原子操作，是实现所有自定义原子逻辑的基础。

• atomicCAS(int\* address, int compare, int val)
• 支持类型：int, unsigned int, unsigned long long, unsigned short int。

核心逻辑回扣： 只有当 *address == compare 时，才会执行 *address = val。无论是否成功，都返回 *address 的旧值。

7.2.3.3 原子操作的成本

1. 硬件路径成本（访问延迟）

在现代 NVIDIA GPU 架构中，原子操作并不是在多处理器（SM）内部的寄存器中完成的，而是在 L2 缓存（L2 Cache） 层面处理的。

• 普通操作：线程读取数据到寄存器，计算，然后写回。
• 原子操作：SM 发送一个原子指令请求给 L2 缓存控制器。L2 控制器负责锁定该缓存行、执行加法/比较逻辑，然后再更新。
• 代价：由于数据必须跨越 SM 到达 L2 缓存，其**往返延迟（Round-trip Latency）**远高于寄存器操作。

2. 串行化成本（争用瓶颈）

这是原子操作最昂贵的开销。当成千上万个线程尝试同时修改同一个内存地址时，硬件无法并行处理。

• 冲突（Contention）：如果一个 Warp（32个线程）中的所有线程都对同一个地址执行 atomicAdd，这些操作将被串行化。这意味着原本 1 个时钟周期能完成的事，现在至少需要 32 个周期。
• 吞吐量下降：随着竞争同一个地址的线程数增加，GPU 的并行优势会线性下降，最终演变成极其低效的单线程模式。

3. 指令流水线成本

不同的原子操作对指令流水线的影响不同：

• 专用硬件支持：像 atomicAdd (int) 这种常用指令，硬件有专门的电路支持，效率较高。
• CAS 循环成本：如果使用 atomicCAS 来实现自定义原子操作（如前面的原子乘法），它涉及到一个 do-while 循环。如果冲突频繁，线程会不断重复“读取-比较-失败-重试”的过程，这会极大地消耗指令吞吐量。

操作类型	相对延迟	并行度	备注
寄存器计算	1x (极低)	极高	GPU 满载运行
无冲突原子操作	10x - 50x	高	访问不同地址，L2 负载均衡
高冲突原子操作	100x - 1000x+	极低	所有线程排队等待同一个地址

对atomic-ordering.cu，创建并运行相关程序。这个小应用程序包含了两个核函数，分别叫作atomics和unsafe。在一个共享变量上，atomics核函数在每个线程上执行原子加法，保存目标地址中原来存储的数值。

对atomics核函数来说，执行原子加法

__global__ void atomics(int *shared_var, int *values_read, int N, int iters)
{
    int i;
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid >= N) return;

    values_read[tid] = atomicAdd(shared_var, 1);

    for (i = 0; i < iters; i++)
    {
        atomicAdd(shared_var, 1);
    }
}

对unsafe核函数来说，执行相同的加法，但并不使用原子函数。

__global__ void unsafe(int *shared_var, int *values_read, int N, int iters)
{
    int i;
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid >= N) return;

    int old = *shared_var;
    *shared_var = old + 1;
    values_read[tid] = old;

    for (i = 0; i < iters; i++)
    {
        int old = *shared_var;
        *shared_var = old + 1;
    }
}

编译并执行atomic-ordering.cu

nvcc -arch=sm_120 atomic-ordering.cu -o atomic-ordering
./atomic-ordering

得到以下输出

(base) root@63320d55395d:~/cuda_learn/cuda_program/CodeSamples/chapter07# ./atomic-ordering 
In total, 30 runs using atomic operations took 0.018436 s
  Using atomic operations also produced an output of 6400064
In total, 30 runs using unsafe operations took 0.002244 s
  Using unsafe operations also produced an output of 100001
Threads performing atomic operations read values 0 1 2 3 4 5 6 7 8 9
Threads performing unsafe operations read values 0 0 0 0 0 0 0 0 0 0

期望结果：总共有 $N=64$ 个线程，每个线程执行 iters + 1 次加法（循环内 iters 次，循环外 1 次）。所以期望值是 $64 \times (100000 + 1) = \mathbf{6,400,064}$。

原子操作 (atomics)：输出 6,400,064。结果完全正确。

不安全操作 (unsafe)：输出 100,001。结果错得离谱（只有期望值的约 1.5%）。

Atomic 时间：0.018436 s

Unsafe 时间：0.002244 s

结论：原子操作慢了约 8.2 倍。

Unsafe：GPU 线程完全不顾及彼此，直接向 L2 缓存发起写请求，硬件流水线可以满载运行，代价是数据混乱。

Atomic：正如我们之前讨论的，atomicAdd 强制执行了 L2 缓存层面的串行化。64 个线程在同一个内存地址上“排队”。由于这种排队竞争，执行速度大幅下降。

7.2.3.4 限制原子操作的性能成本

当必须执行原子操作时，使用有些方法可以减少性能损失。你可以使用局部操作来增强全局原子操作，这些局部操作能从同一线程块的线程中产生一个中间结果。这需要使用本地较低延迟的资源，如shuffle指令或共享内存，在使用原子操作把局部结果结合到最终全局结果之前，需要先从每个线程块产生局部结果。

1. 局部规约：使用共享内存

这是最常用且最高效的优化手段。

• 原理：与其让数万个线程直接竞争全局内存（L2 缓存）的一个地址，不如让每个线程块（Block）先在自己内部的共享内存里进行累加。
• 优势：共享内存在 SM 内部，延迟极低，且冲突仅限制在块内（最多 1024 个线程）。最后，每个块只需向全局内存提交一次结果。

2. 利用 Warp 级原语

在线程访问全局内存之前，先利用 Warp Shuffle (__shfl_sync) 指令在 Warp 内部（32 个线程）进行求和。

• 做法：一个 Warp 的 32 个线程先互相交换数据，计算出总和。
• 结果：由该 Warp 的“带头线程”（Leader thread）执行一次 atomicAdd。这直接将全局内存的冲突压力降低了 32 倍，且不消耗任何内存带宽，完全在寄存器层面完成。

3. 分散热点：填充与分桶

如果你的算法是在做直方图统计，且某些“桶”极其繁忙（例如图像中背景色对应的像素极多）：

• 策略：为热点地址建立多个副本。例如，创建 4 个计数器副本，线程根据 tid % 4 随机选择一个进行原子加法。
• 代价：最后需要额外的一步来汇总这些副本的值，但这通常比在高冲突下等待要快得多。

4. 优化内存作用域

在现代 NVIDIA 架构（Volta 及以后）中，使用更精准的作用域标志：

• atomicAdd_block：如果你只需要在当前线程块内保证原子性，使用块级作用域可以告诉硬件不需要同步其他 SM 的缓存，从而减少一致性协议的开销。

5. 算法层面的权衡：减少原子操作频率

• 合并操作：检查代码逻辑，看是否可以将多次小幅度的原子加法合并为一次大的原子加法。
• 混合模式：对于大多数数据，先尝试非原子的局部累加；仅在处理跨块边界或最后合并阶段才调用原子指令。

7.2.3.5 原子级浮点支持

早于 Pascal 架构 (Compute Capability < 6.0)：

• 仅原生支持单精度浮点数 float 的 atomicAdd()。
• 双精度 (double)：完全没有硬件指令支持。开发者必须使用 atomicCAS() 指令手动编写循环来模拟双精度原子加法。

Pascal 架构 (CC 6.0, 如 P100/GTX 10系列)：

• 里程碑：首次引入了原生硬件支持的 双精度 (double) atomicAdd()。
• 这极大地提升了高性能计算（HPC）应用的速度，因为不再需要慢速的 CAS 循环。

Ampere 架构 (CC 8.0, 如 A100)：

• 引入了对 半精度 (__half 和 __half2) 的原生原子加法支持。
• 这对深度学习（尤其是 FP16 训练）的梯度累加非常关键。

浮点原子操作在硬件实现上面临两个挑战：

1. 非结合律 (Non-associativity)：

   在数学上，$(a + b) + c = a + (b + c)$。但在浮点运算中，由于舍入误差的存在，运算顺序的不同会导致结果微小的差异。

   • 原子操作的问题：原子操作不保证线程执行顺序。这意味着在并发竞争下，即便输入数据相同，两次运行得到的浮点总和可能会有极细微的差别（最后几个比特位）。

2. 电路复杂性：

   整数加法器很简单，但浮点加法需要对阶、尾数求和、规范化等步骤，在 L2 缓存层面集成浮点逻辑单元（ALU）需要消耗更多的硬件晶体管。

7.3 总结

本章探讨了 CUDA 指令级原语对程序性能与准确性的影响，重点分析了单双精度浮点数在表示能力与硬件吞吐量上的巨大差异，揭示了标准函数与硬件加速内置函数（Intrinsic）在 IEEE 754 合规性与 SFU 执行速度之间的权衡，并详细阐述了原子指令（尤其是 CAS）在保证多线程数据安全时的硬件实现成本及通过局部规约（如共享内存和 Warp Shuffle）进行优化的策略。