访存优化_4、多维数组，插桩与循环分块、莫顿码、多核缓存伪共享

多维数组

C语言中的静态数组：

• float a[n]; 可以在栈上分配有 n 个元素的 一维 数组。

• 通过 a[ i ] 访问第 i 个元素。

• float a[n][m]; 可以在栈上 分配 n 行 m 列的二维数组。

• 通过 a[ i ][j] 访问第 i 行，第 j 列的元素。

 

众所周知，内存是一维的，因此任何二维数组，都必须被扁平化，才能储存在内存中。

对于 float a[3][4] 编译器实际上会把他变成一维数组 float a[3*4] ，然后把 a[ i ][j] 翻译为 a[ i * 4 + j] 。

C语言中的动态数组：

• float *a = malloc(n * sizeof (float)); 可以在堆上分配有 n 个元素的一维数组。

• 通过 a[ i ] 访问第 i 个元素。

• float *a = malloc(n * m * sizeof (float)); 可以在堆上 分配 n 行 m 列的二维数组。

• 通过 a[ i * m + j] 访问第 i 行，第 j 列的元素。

• 释放时，统一用 free(a)

 而C++由于RAII思想，不需要free，在作用域结束后自动释放内存。

• vector<float> a(n); 可以在堆上分配有 n 个元素的一维数组。

• 通过 a[ i ] 访问第 i 个元素。

• vector<float> a(n * m); 可以在堆上 分配 n 行 m 列的二维数组。

• 通过 a[ i * m + j] 访问第 i 行，第 j 列的元素。

二级指针 ！= 二维数组

float **a = malloc(n * sizeof(float *));
for (int i = 0; i < m; i++) a[i] = malloc(m * sizeof(float));
// 通过 a[i][j] 访问第 i 行，第 j 列的元素。
for (int i = 0; i < m; i++) free(a[i]);
free(a);

这样做了的话内存都是分散的，每次都要解开两层指针，非常低效，而且要用m+1次malloc。

分配 n*m 二维数组，正确的方式永远是：float *a = malloc(n * m * sizeof(float));

也不要用 vector<vector<float>>，这个是一样的，甚至更烂——vector 本身是24字节。

但是 float [n][m] 和 array<array<float, n>, m> 就没问题，因为他们是静态大小，编译器可以检测到并自动扁平化，转换成乘法和加法来计算地址。

 为什么说二级指针是低效的储存和索引方式：

        如上图所示，二级指针指进去是一个指针的指针，这个指针的列表再去访问a[3]再去访问a[3][0]。实际上是两次内存的读取。而一个a[i*m+j]只涉及到一次内存的读取。而且乘法加法这种运算是要比访存快的多的。

        

从内存的分配上也可以看出自动扁平化的方法快得多  ，这是因为std::vector<std::vector<std::vector<float>>> 实际上是三维数组，它在内存中不是连续存储的，而是需要额外的指针来维护各个维度之间的地址关系。这意味着每次访问三维数组的元素时需要进行多次指针解引用操作，而这些操作会消耗较多的时间。

二维数组的行主序和列主序：

定义：row-major中，在同一行的元素在内存中是相邻的；column-major中，同一列的元素在内存中是相邻的。

二维数组因为涉及到扁平化的问题，所以分为了两种方法。

以下符号约定：i 行号，j 列号；n 行数，m 列数

 C/C++ 编译器把静态数组 a[i][j] 翻译为 a[i * m + j]，所以是行主序。

Fortran 等会把矩阵 A(i, j)  翻译为 a[i + j * n]，所以是列主序。

手动扁平化时，如果用 a[i * m + j] 就是行主序。

手动扁平化时，如果用 a[i + j * n] 就是列主序。

简单来说：哪个索引连续，就是什么主序。

        实际上，如果你认为 C++ 中用 a[j][i] 来访问矩阵元素的话，那 C++ 就是列主序了，索引出现的先后顺序并不是重点，索引的意义才重点。

 /*   我总搞混的行主序和按列读取：

 在CUDA当中也是按照行主序，行优先的模式存储

 所以在cuda当中图a说明的是按行读取的模式，即在第一次循环中warp读取0-31行中列索引为0的元素，在第二次循环中读取列索引为1的元素，而这种读取方法在C语言中是隔着地址读取的.图b说明的是按列读取的模式，即在第一次循环中warp读取0-31列中行索引为1的元素。见上图可知这种读取是连续的。更加高效。

*/

constexpr size_t nx = 1 << 13;
constexpr size_t ny = 1 << 13;

std::vector<float> a(nx* ny);

void BM_yx_loop_yx_array(benchmark::State& bm) {
    for (auto _ : bm) {
#pragma omp parallel for collapse(2)
        for (int y = 0; y < ny; y++) {
            for (int x = 0; x < nx; x++) {
                a[y * nx + x] = 1;
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_yx_loop_yx_array);

void BM_xy_loop_yx_array(benchmark::State& bm) {
    for (auto _ : bm) {
#pragma omp parallel for collapse(2)
        for (int x = 0; x < nx; x++) {
            for (int y = 0; y < ny; y++) {
                a[y * nx + x] = 1;
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_xy_loop_yx_array);

void BM_yx_loop_xy_array(benchmark::State& bm) {
    for (auto _ : bm) {
#pragma omp parallel for collapse(2)
        for (int y = 0; y < ny; y++) {
            for (int x = 0; x < nx; x++) {
                a[y + x * ny] = 1;
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_yx_loop_xy_array);

void BM_xy_loop_xy_array(benchmark::State& bm) {
    for (auto _ : bm) {
#pragma omp parallel for collapse(2)
        for (int x = 0; x < nx; x++) {
            for (int y = 0; y < ny; y++) {
                a[y + x * ny] = 1;
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_xy_loop_xy_array);

这是因为 YX 序的数组，X 方向在内存空间中的排布是连续的。

而 YX 序的循环，其 X 是内层循环体，因此在先后执行的时间上是连续的。

如果是 XY 序的循环，其 X 是外层循环体，在先后执行的时间上是不连续的。

从而在硬件看来，以 YX 序遍历，就和顺序访问一维数组没什么两样，从而缓存预取能正常运作，甚至编译器可以优化成一个 nx*ny 的一层循环。

 可以看到连续（合并）访问真的很重要。

 ndarray类：

原作者在这里提供了ndarray类：

用 ndarray<2, float> 声明一个二维浮点数组，ndarray<3, int> 声明一个三维整型数组。

这里的 ndarray 通过 a(x, y) 来索引，看起来像 Fortran，但是实际上还是 YX 序，和静态数组的 a[y][x] 一样的，指标出现的顺序并无关紧要。

插桩与循环分块： 

插桩这个操作，就是说在结构网格中，从一个点往周围固定范围读取值，并根据一定权重累加，然后用于修改自身的值。

比如求一个场的梯度、散度、旋度、拉普拉斯。如果场是用结构网格（structured grid ）表示，那就是一个插桩操作。

插桩的内核（ kernel ）指的就是这个“周围范围”的形状（如右图三个例子）和每个地方读取到值对修改自身值的权重等信息。

有点类似于神经网络中卷积的概念。

 这是一个X方向模糊的程序，nblu*2+1的顺序读取操作的长度相对较短，执行速度很快。当该顺序读取操作执行完毕后，CPU的预测机制可能会将注意力转移到其他可能的指令上，而不再关注该顺序读取的下一步。

因此，当下一步执行的指令需要依赖于顺序读取操作的结果时，CPU就无法预测下一步要读取的位置。这是因为顺序读取操作的完成相对较快，导致CPU无法提前加载所需的数据到缓存中，从而影响了预测能力。

在这种情况下，CPU可能需要等待数据从内存中加载到缓存中，这可能会导致额外的延迟和性能下降。

所以作者用了prefetch进行预取，提前告知程序下一步要取什么。但是也可以见到这样做之后速度变慢了许多。这是因为没进新货一次x++，系统就进行一次预取，但是预取一下子取了32个元素，就造成了性能的浪费。

于是作者将原来X的循环进行了分块。

void BM_x_blur_prefetched(benchmark::State &bm) {
    for (auto _: bm) {
#pragma omp parallel for collapse(2)
        for (int y = 0; y < ny; y++) {
            for (int x = 0; x < nx; x++) {
                _mm_prefetch(&a(x + 16, y), _MM_HINT_T0);
                float res = 0;
                for (int t = -nblur; t <= nblur; t++) {
                    res += a(x + t, y);
                }
                b(x, y) = res;
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_x_blur_prefetched);

void BM_x_blur_tiled_prefetched(benchmark::State &bm) {
    for (auto _: bm) {
#pragma omp parallel for collapse(2)
        for (int y = 0; y < ny; y++) {
            for (int xBase = 0; xBase < nx; xBase += 16) {
                _mm_prefetch(&a(xBase + 16, y), _MM_HINT_T0);
                for (int x = xBase; x < xBase + 16; x++) {
                    float res = 0;
                    for (int t = -nblur; t <= nblur; t++) {
                        res += a(x + t, y);
                    }
                    b(x, y) = res;
                }
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_x_blur_tiled_prefetched);

分块后就是没前进16个元素再进行一次预取。而且每一个x的小区间还是[xbase , xbase + 16)

为什么不用if(x%16)是因为要避免分支。

也可以采取stream直写进一步优化。

但通过对比可以发现，Y方向的插桩比X方向插桩慢好多。

因为 X 方向的插桩所读取的数据，在内存中是连续的。

而 Y 方向的插桩所读取的数据，在内存看来表现为跳跃 nx 大小来访问，是不连续的。

 上图为X方向的插桩，展开为一维数组是这样的

 

 而Y方向的插桩是不连续的，是跳跃着访问的。

void BM_y_blur(benchmark::State &bm) {
    for (auto _: bm) {
#pragma omp parallel for collapse(2)
        for (int y = 0; y < ny; y++) {
            for (int x = 0; x < nx; x++) {
                float res = 0;
                for (int t = -nblur; t <= nblur; t++) {
                    res += a(x, y + t);
                }
                b(x, y) = res;
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_y_blur);

void BM_y_blur_tiled(benchmark::State &bm) {
    for (auto _: bm) {
        constexpr int blockSize = 32;
#pragma omp parallel for collapse(2)
        for (int yBase = 0; yBase < ny; yBase += blockSize) {
            for (int xBase = 0; xBase < nx; xBase += blockSize) {
                for (int y = yBase; y < yBase + blockSize; y++) {
                    for (int x = xBase; x < xBase + blockSize; x++) {
                        float res = 0;
                        for (int t = -nblur; t <= nblur; t++) {
                            res += a(x, y + t);
                        }
                        b(x, y) = res;
                    }
                }
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_y_blur_tiled);

对于这种不连续的访问情况，可以用循环分块来解决：

这样的话就可以把小方块装进缓存里了

但是这种方法存在问题：比如 插桩的读取范围跨越了 block 的边界，反而会变得跨度更大，导致缓存彻底失效。 

 应对这种，就是只对X进行分块：

插桩的读取范围跨越了 block 的边界，反而会变得跨度更大，导致缓存彻底失效。

 

void BM_y_blur_tiled_only_x(benchmark::State &bm) {
    for (auto _: bm) {
        constexpr int blockSize = 32;
#pragma omp parallel for collapse(2)
        for (int xBase = 0; xBase < nx; xBase += blockSize) {
            for (int y = 0; y < ny; y++) {
                for (int x = xBase; x < xBase + blockSize; x++) {
                    float res = 0;
                    for (int t = -nblur; t <= nblur; t++) {
                        res += a(x, y + t);
                    }
                    b(x, y) = res;
                }
            }
        }
        benchmark::DoNotOptimize(a);
    }
}
BENCHMARK(BM_y_blur_tiled_only_x);

同理，对于矩阵的转置也是这样的

 

具体可以看一下上面这张图，最开始A是顺序的读取，所以A是高效的，但B每次都要重新读取，是低效的（每次都要读取一整个缓存行）。在经过分块之后，A还是高效的，B也变成高效率的了。而且只要二者的 分块大小的和小于Cache的大小，那么就可以保证高效。 

但是这样也会有一个问题，就是如果编程序的人的块分的比较小，那么Cache就不能全部利用上。

所以开发了莫顿码：

        行主序是指按照行来存储多维数组或矩阵的元素。在行主序中，相邻元素在内存中的地址是相邻的，但对于二维及更高维的数组来说，具有相邻位置的元素在逻辑上可能相距较远，因此具有较差的局部性。这是因为在行主序中，相邻的元素可能在内存中的不同页或不同缓存行中，导致缓存未命中（cache miss）增加，降低了内存访问效率。

        莫顿码（Morton Order），也称为Z序编码，是一种将多维空间中的点映射到一维序列的方法。莫顿码通过将多维空间划分为一系列二进制位，然后将每个维度的坐标值编码为对应的二进制位，从而将多维坐标转化为一个整数。在莫顿码中，具有相邻位置的元素在逻辑上也是相邻的，因此具有较好的局部性。在内存中，莫顿码相邻的元素通常也会在物理上靠近，因此内存访问更加连续，减少了缓存未命中的可能性，提高了内存访问效率。

 

• 如 x=x1x2x3, y=y1y2y3

• 则他们的莫顿码： m( x,y )=y1x1y2x2y3x3

• 二维莫顿编码可以把两个长度为 n 的二进制数，交错打包成一个长度 2*n 的二进制数。而莫顿编码的逆运算，就是莫顿解码：

• mdec (m1m2m3m4)=(m2m4, m1m3)

   

莫顿码的几何意义在于，以 ( x,y )= mdec (t) 为参数方程，可以生成一条分形（自相似）的 Z 字型曲线。 该曲线的有一个 性质，上面两个 ( x,y ) 坐标相近的点，他们的 t 也相近（大概率）。

t可以近似看为时间，可以视为空间上相邻的点，在时间上也相近。

意义：可以用一维的 t 遍历二维的网格，然后用 mdec (t) 求出要访问元素的 ( x,y ) 坐标，这样可以保证的数据在时间 t 上是接近的，同时二维空间上 ( x,y ) 也是接近的，有利于访存局域性

 

 所以在循环中使用莫顿码来进行循环，最大的好处就是他是自适应的。

 如上图，可以将1视为一个块，或者1234视为一个块，或者上图的全部都视为一个块。

         按照Z曲线遍历，首先只需要一个一维循环，其循环变量就是莫顿码，区间范围则是 [0, nx*ny/blockSize^2)。

        然后，通过莫顿解码，获取 X，Y 分量。

 最后 TBB自带莫顿遍序。

示例：矩阵乘

• a( i , j) 始终在一个地址不动（一般）。

• b( i , t) 每次跳跃 n 间隔的访问（坏）。

• c(t, j) 连续 的 顺序访问（好）。

• 因为存在不连续的 b 和一直不动的 a ，导致矢量化失败，一次只能处理一个标量， CPU 也无法启动指令级并行（ ILP ）。

• a( i , j) 连续 32 次顺序访问（好）。

• b( i , t) 连续 32 次顺序访问 （ 好 ）。

• c(t, j) 32 次 在一个地址不动 （ 一般 ）。

• 这样就消除不连续的 访问了 ，从而内部的 i 循环可以顺利矢量化，且多个循环体之间没有依赖关系， CPU 得以启动指令级并行，缓存预取也能正常工作，快好多！

还有小内核卷积：

 最后unroll一下会更好。

多核下的缓存：

 

• 需要注意，如果多个核心同时访问的地址非常接近，这时候会变得很慢！

• 这是因为 CPU 之间通信的最小单位也是 缓存行（ 64 字节），如果两个核心访问到了的同一缓存行，假设一个核心修改了该缓存行的前 32 字节，另一个修改了后 32 字节，同时写回的话，结果要么是只有前 32 字节，要么是只有后 32 字节，而不能两个都正确写入。

• 所以 CPU 为了安全起见，同时只能允许一个核心写入同一地址的缓存行。从而导致读写这个变量的速度受限于三级缓存的速度，而不是一级缓存的速度。

 

• 要想消除 错误共享很简单，只需要把每个核心写入的地址尽可能分散开了就行了。比如这里，我们把每个核心访问的地方跨越了 16KB ，这样 CPU 就知道每个核心之间不会发生冲突，从而可以放心地放在一级缓存里，不用担心会不会和其他核心共用了一个缓存行了。

不过错误共享只会发生在写入的情况，如果多个核心同时读取两个很靠近的变量，是不会产生冲突的，也没有性能损失