请教：数组索引连续化的性能问题

[chenxinfeng4]

我又来打扰了。我记得小彭老师说过，连续的内存访问速度比较快。可是我的高维数组进行抽提时，低维也有很长片段的连续，为啥还是很慢。需要用什么 stream 或者 prefetch 来优化吗？

# array0 是一个体像素(voxel)，对其进行3D上的crop。crop 的尺寸为 64x64x64
array2 = np.ascontiguousarray(array0[:9, x:x+64, y:y+64, z:z+64])  #280fps, 4.9 GByte/s 
array3 = C++版本的memcopy(array, x, y, z, ....)   #同样 280fps

这个 4.9 GByte/s 的内存访问速度就很离谱，我还没有做任何计算。

https://github.com/chenxinfeng4/hpc_issue_array_continguous

[archibate]

我看一下。

很长片段的连续

只有 64 个元素的连续，怎么好意思叫“很长片段”的？一般预取都可以连续一整页（4KB）。

:9

这是什么？4D 体素切片？咱们不是说好了用 HWC 嘛？

[archibate]

上一次和你说SOA好，但是对于这种切片操作，最实惠的反而又是AOS了。

优化前（CXYZ）

100%|████████████| 1000/1000 [00:02<00:00, 369.85it/s]
100%|████████████| 1000/1000 [00:02<00:00, 387.13it/s]

优化后（XYZC）

100%|████████████| 1000/1000 [00:01<00:00, 503.44it/s]
100%|████████████| 1000/1000 [00:01<00:00, 508.49it/s]

原理：本来一次只能连续拷贝 64 个元素的（Z），把 C 排到 Z 后面，就可以一次连续拷贝 64*9 个元素了（ZC）。
于是，从 6.8 GB/s 提升到了 8.9 GB/s（我的内存理论带宽上限是 41.6 GB/s，因为拷贝是双向的所以拷贝的上限是 20.8 GB/s，不用 stream 的话则是三倍带宽，理论上限只有 13.8 GB/s，我再尝试 stream + prefetch 优化下）。

[archibate]

无法使用 stream 或 prefetch，因为 numpy 不提供让数组对齐到 64 字节（缓存行）的方式。

[archibate]

可以了

#include <cstdio>
#include <cstdint>
#include <cstring>
#include <type_traits>
#include <immintrin.h>
#include <stdexcept>

template <int Offset>
static void memcpy_streamed(char *dst, char const *src, char const *nextSrc, size_t size) {
    const size_t width = 64;
    /* auto chDst = (char *)dst; */
    /* auto chSrc = (char const *)src; */
    /* auto chDstTemp = (char *)(uintptr_t)chDst */
    /* auto chDstEnd = (char *)dst + size; */
    /* while (size-- % width) { */
    /*     *chDst++ = *chSrc++; */
    /* } */
    // 如果 Numpy 的数组首地址离对齐到 64 字节差 16 字节
    // 为了对齐到缓存行以便使用 stream 和 prefetch，需要先把前面 3 个 16 字节当作边角料处理掉
    if (Offset) {
        for (int i = 0; i < 4-Offset; i++) {
            _mm_stream_si128((__m128i *)dst, _mm_loadu_si128((__m128i *)src));
            dst += 16;
            src += 16;
        }
        size -= Offset * 16;
    }
    /* if ((uintptr_t)dst % 64) [[unlikely]] { */
    /*     printf("%d, dst: %zd\n", Offset, (uintptr_t)dst % 64); */
    /*     throw std::runtime_error("not aligned dst"); */
    /* } */
    int batch = size / width;
    auto vecDst = (__m256i *)dst;
    auto vecSrc = (__m256i const *)src;
    auto vecDstEnd = vecDst + batch * width / sizeof(__m256i);
    while (vecDst != vecDstEnd - 2) {
        auto v1 = _mm256_loadu_si256(vecSrc);
        auto v2 = _mm256_loadu_si256(vecSrc + 1);
        _mm256_stream_si256(vecDst, v1);
        _mm256_stream_si256(vecDst + 1, v2);
        vecSrc += 2;
        vecDst += 2;
    }
        auto v2 = _mm256_loadu_si256(vecSrc + 1); // 故意反过来破坏当前预取序列
        auto v1 = _mm256_loadu_si256(vecSrc);
        _mm_prefetch(nextSrc, _MM_HINT_T0); // 告诉他接下来该处理下一行了
        _mm256_stream_si256(vecDst + 1, v2);
        _mm256_stream_si256(vecDst, v1);
        vecSrc += 2;
        vecDst += 2;
    // 还需要也把后面 1 个 16 字节当作边角料处理掉
    dst = (char *)vecDst;
    src = (char const *)vecSrc;
    if (Offset) {
        for (int i = 0; i < Offset; i++) {
            _mm_stream_si128((__m128i *)dst, _mm_loadu_si128((__m128i *)src));
            dst += 16;
            src += 16;
        }
    }
    _mm_prefetch(nextSrc + width, _MM_HINT_T0);
}

template<typename T>
void concat(T *source, T *destination, const size_t DIM_SRC_N,
            const size_t DIM_SRC_X, const size_t DIM_SRC_Y, const size_t DIM_SRC_Z,
            const size_t DIM_DST_X, const size_t DIM_DST_Y, const size_t DIM_DST_Z,
            const size_t START_X,   const size_t START_Y,   const size_t START_Z
            )
{
    // 开启 openmp 多线程加速, 好像2个线程就饱和了 // 小彭老师：我靠，memcpy是典型的内存瓶颈（membound）操作，当然没法用并行加速了
    // #pragma omp parallel for num_threads(2)
    // T (*vectorSRC)[DIM_SRC_Y][DIM_SRC_Z] = reinterpret_cast<T(*)[DIM_SRC_Y][DIM_SRC_Z]>(source);
    // T (*vectorDST)[DIM_DST_Y][DIM_DST_Z] = reinterpret_cast<T(*)[DIM_DST_Y][DIM_DST_Z]>(destination);

    size_t size_cp = DIM_SRC_N * DIM_DST_Z * sizeof(T);
    if (size_cp % 64) throw std::runtime_error("cannot handle bianjiaoliao for now");
    // 实测发现 Numpy 的数组首地址总是对齐到 16 字节
    if ((uintptr_t)destination % 16) { printf("%zd\n", (uintptr_t)destination % 64); throw std::runtime_error("dst must be 16B-aligned"); }

    auto threadconcat = [&] (auto offset) {
        #pragma omp parallel for num_threads(4)
        for (size_t ix = 0; ix < DIM_DST_X; ix++){
            for (size_t iy = 0; iy < DIM_DST_Y; iy++){
                // T * ptrSRC = &vectorSRC[ix][iy][0];
                // T * ptrDST = &vectorDST[ix+START_X][iy+START_Y][START_Z];
                T * ptrSRC = &source[DIM_SRC_N * ((ix + START_X) * DIM_SRC_Y * DIM_SRC_Z + (iy + START_Y) * DIM_SRC_Z + START_Z)];
                T * ptrNextSRC = &source[DIM_SRC_N * (((ix + (iy + 1) / DIM_DST_Y) + START_X) * DIM_SRC_Y * DIM_SRC_Z + ((iy + 1) % DIM_DST_Y + START_Y) * DIM_SRC_Z + START_Z)];
                T * ptrDST = &destination[DIM_SRC_N * (ix * DIM_DST_Y * DIM_DST_Z + iy * DIM_DST_Z)];
    /* printf("!!!pd%zd\n", (uintptr_t)ptrDST%64); */
                memcpy_streamed<offset.value>((char *)ptrDST, (char *)ptrSRC, (char *)ptrNextSRC, size_cp);
                // for (int iz = 0; iz < DIM_DST_Z; iz++){
                //     ptrDST[iz] = ptrSRC[iz];
                // }
            }
        }
    };

    int offset = (uintptr_t)destination % 64 / 16;
    /* printf("!!!of%d\n", offset); */
    switch (offset) {
    case 0: threadconcat(std::integral_constant<int, 0>()); break;
    case 1: threadconcat(std::integral_constant<int, 1>()); break;
    case 2: threadconcat(std::integral_constant<int, 2>()); break;
    case 3: threadconcat(std::integral_constant<int, 3>()); break;
    };
}

[archibate]

不并行：

100%|████████████| 1000/1000 [00:02<00:00, 498.93it
100%|████████████| 1000/1000 [00:01<00:00, 775.83it

4线程并行：

100%|██████████| 1000/1000 [00:01<00:00, 520.58it/s]
100%|██████████| 1000/1000 [00:01<00:00, 977.46it/s]

达到 17.2 GB/s 接近 20.8 GB/s 理论上限。

[archibate]

哦还有你这个是不是knn之类的，START_X是任意数的，如果你能保证START_X都是64的整数倍，可以把src数组block化，这样就完全没有跳跃访问了。

[chenxinfeng4]

感谢解答。
HWC 的图像格式，我现在退化为HW，只考虑灰度信息。因为我的视频数据原始是灰色，在后处理各方面都快一些。

我说一下我用这个是做什么。我提取切片用来做深度学习，切片尺寸为KXYZ (K_64_64_64).

下图是示意图。可以看到下面的这个框的坐标点就是 (K_64_64_64) 【第一步】。接着从这个2D的 KHW 图像中，按坐标索引出 (K_64_64_64) 的体像素【第二步】。 总之，这个 (K_64_64_64)是将2D图像升维到3D，提升了神经网络的识别效果（paper: VoxelPose）。瓶颈就是这个 (K_64_64_64) 的体像素制备速度非常慢，我还在一点点优化。

这个issue 还在【第一步】，还有【第二步】今后我还需要慢慢优化。以后有需求可能再提升到彩色， KCXYZ (K_3_64_64_64)，那时再回顾上一个issue的技巧。

再次感谢解答。

[chenxinfeng4]

但是对于这种切片操作，最实惠的反而又是AOS了

我学到了。切片时，让连续的片段尽可能长，使用（XYZC）会更快。

[chenxinfeng4]

我再尝试 stream + prefetch 优化下

震惊！！原来 std::memcpy 也是可以被优化的吗？因为它没有考虑内存对齐问题？ Numpy 数组也有对齐，可怕。

不并行：

100%|████████████| 1000/1000 [00:02<00:00, 498.93it
100%|████████████| 1000/1000 [00:01<00:00, 775.83it

从一开始 numpy 的CXYZ 的 387 it, 通过 XYZC + stream + prefetch 提升到 775it。提升到 2x 速率，太生猛了。

[chenxinfeng4]

哦还有你这个是不是knn之类的，START_X是任意数的，如果你能保证START_X都是64的整数倍，可以把src数组block化，这样就完全没有跳跃访问了。

START_X的确是任意数。

我的所有操作都是内存 I/O 密集的任务。没有计算。我记得小彭老师介绍cuda编程的时候说过N卡的显存速度快与主存。我输出的数据最后也要搬到显存上做深度学习，那我是不是可以直接在cuda上做切片。。。小彭老师能预估一下这个性能提升大概有多大？

[archibate]

原来 std::memcpy 也是可以被优化的吗？

memcpy 只有在数据量大于 4096 时才会尝试用 stream 指令，否则依然会使用 store 指令。而你这里都是一小段一小段分别 memcpy，他是不知道你是一个大数据量拷贝的。

[archibate]

我的所有操作都是内存 I/O 密集的任务。没有计算。我记得小彭老师介绍cuda编程的时候说过N卡的显存速度快与主存。我输出的数据最后也要搬到显存上做深度学习，那我是不是可以直接在cuda上做切片。。。小彭老师能预估一下这个性能提升大概有多大？

显存确实比内存快，但是从内存到显存，却比两者都慢！因为需要经过 PCIe 总线。

如果说：CPU=湖泊，PCIe=河流，GPU=大海。

你的问题就变成：我有很多货物，但是有很大一部分是要扔掉的。是提前在湖泊里挑选出需要的货物留下，在运河里开精简的货物，还是先把所有的货物一次性运到大海里再进行挑选。显然河流是瓶颈，在湖泊里提前挑选好可以尽可能降低河流的负担。

设显存到显存的时间成本为1，内存到内存的时间成本为2，内存到显存的时间成本为3。

所以如果你是一次性切片的话，那么两种方案的速度比较如下：

• 在 CPU 上切片好后拷给 GPU 计算：264 + 364 + 1*64
• CPU 直接全拷给 GPU 切片并计算：2164 + 3164 + 1*64

显然在 CPU 上提前切片比较划算，毕竟切片不是计算密集型任务。

总结：GPU to GPU > CPU to CPU > CPU to GPU > GPU to CPU

其中 CPU to GPU 比 GPU to CPU 还要稍微快一点，因为后者需要强制同步：https://stackoverflow.com/questions/50302112/why-is-it-faster-to-transfer-data-from-cpu-to-gpu-rather-than-gpu-to-cpu

但是有一个例外，就是如果你需要对一个 164 体素，多次切片的话，那么两种方案的速度比较如下：

• 在 CPU 上切片好后拷给 GPU 计算：(2*64 + 3*64) * n + 1*64
• CPU 直接全拷给 GPU 切片并计算：2*164 + 3*164 + n*1*64

也就是每一个集装箱都有用，但是要派发给不同的国家，这时就是要先全部挪到大海里，再统一重新分配一下高效，不提前全部通过河流，每挑选一次就需要通过河流一次，反而不划算了。

[archibate]

哦对了，你这个果真是uint64的数据？如果不需要，可以改用uint32的，还有你这个图片，看起来像是可以变成稀疏化体素的，不知道是自己写还是因为别的第三方库就要求这种XYZC的稠密格式？有没有办法改用专业的openvdb？

[archibate]

还有一件事，你知不知道你这个切片其实还是有一定计算的，那就是整数索引的乘法计算，因为你的数组大小都不保证是2的幂次方，不仅没有对齐保障也需要一直计算整数乘法（虽然最后还会是内存瓶颈）。使用专业的openvdb稀疏存储，那么每个儿子块都是8x8x8，可以用高效的>>3和&7来实现三维索引线性化。openvdb最重要的一点是他可以把全部为0的8x8x8区块用一个空指针表示，对于你这种只有一小块有东西，大部分是0的体素能够节省不少空间，可能你是要伺候第三方库，只提供了稠密的固定CXYZ的接口？

[chenxinfeng4]

这是因为我们的体像素本身就是要稠密的，后续的深度学习也需要用稠密数组。 虽然你看到图像有很多黑像素，但是在 那个立方形感兴趣框里面，大部分是有效像素。说体像素，一般是说ct 那种透视扫描的，是真实的。但我们是用2d图像投影到3d空间，是类似投影仪的光线在空气中，是一种伪体像素。但依然是稠密的。

…

---- 回复的原邮件 ---- | 发件人 | ***@***.***> | | 日期 | 2023年08月05日 18:22 | | 收件人 | ***@***.***> | | 抄送至 | ***@***.***>***@***.***> | | 主题 | Re: [parallel101/simdtutor] 请教：数组索引连续化的性能问题 (Issue #3) | 还有一件事，你知不知道你这个切片其实还是有一定计算的，那就是整数索引的乘法计算，因为你的数组大小都不保证是2的幂次方，不仅没有对齐保障也需要一直计算整数乘法（虽然最后还会是内存瓶颈）。使用专业的openvdb稀疏存储，那么每个儿子块都是8x8x8，可以用高效的>>3和&7来实现三维索引线性化。openvdb最重要的一点是他可以把全部为0的8x8x8区块用一个空指针表示，对于你这种只有一小块有东西，大部分是0的体素能够节省不少空间，可能你是要伺候第三方库，只提供了稠密的固定CXYZ的接口？ — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>