请教：数组索引连续化的性能问题 #3
Comments
|
我看一下。
只有 64 个元素的连续，怎么好意思叫“很长片段”的？一般预取都可以连续一整页（4KB）。
这是什么？4D 体素切片？咱们不是说好了用 HWC 嘛？ |
|
上一次和你说SOA好，但是对于这种切片操作，最实惠的反而又是AOS了。 优化前（CXYZ） 优化后（XYZC） 原理：本来一次只能连续拷贝 64 个元素的（Z），把 C 排到 Z 后面，就可以一次连续拷贝 64*9 个元素了（ZC）。 |
|
无法使用 stream 或 prefetch，因为 numpy 不提供让数组对齐到 64 字节（缓存行）的方式。 |
|
可以了 #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;
};
} |
|
不并行： 4线程并行： 达到 17.2 GB/s 接近 20.8 GB/s 理论上限。 |
|
哦还有你这个是不是knn之类的，START_X是任意数的，如果你能保证START_X都是64的整数倍，可以把src数组block化，这样就完全没有跳跃访问了。 |
|
感谢解答。 我说一下我用这个是做什么。我提取切片用来做深度学习，切片尺寸为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的技巧。 再次感谢解答。 |
我学到了。切片时，让连续的片段尽可能长，使用（XYZC）会更快。 |
震惊！！原来 std::memcpy 也是可以被优化的吗？因为它没有考虑内存对齐问题？ Numpy 数组也有对齐，可怕。
从一开始 numpy 的CXYZ 的 387 it, 通过 XYZC + stream + prefetch 提升到 775it。提升到 2x 速率，太生猛了。 |
START_X的确是任意数。 我的所有操作都是内存 I/O 密集的任务。没有计算。我记得小彭老师介绍cuda编程的时候说过N卡的显存速度快与主存。我输出的数据最后也要搬到显存上做深度学习，那我是不是可以直接在cuda上做切片。。。小彭老师能预估一下这个性能提升大概有多大？ |
memcpy 只有在数据量大于 4096 时才会尝试用 stream 指令，否则依然会使用 store 指令。而你这里都是一小段一小段分别 memcpy，他是不知道你是一个大数据量拷贝的。 |
显存确实比内存快，但是从内存到显存，却比两者都慢！因为需要经过 PCIe 总线。 如果说：CPU=湖泊，PCIe=河流，GPU=大海。 你的问题就变成：我有很多货物，但是有很大一部分是要扔掉的。是提前在湖泊里挑选出需要的货物留下，在运河里开精简的货物，还是先把所有的货物一次性运到大海里再进行挑选。显然河流是瓶颈，在湖泊里提前挑选好可以尽可能降低河流的负担。 设显存到显存的时间成本为1，内存到内存的时间成本为2，内存到显存的时间成本为3。 所以如果你是一次性切片的话，那么两种方案的速度比较如下：
显然在 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 体素，多次切片的话，那么两种方案的速度比较如下：
也就是每一个集装箱都有用，但是要派发给不同的国家，这时就是要先全部挪到大海里，再统一重新分配一下高效，不提前全部通过河流，每挑选一次就需要通过河流一次，反而不划算了。 |
|
哦对了，你这个果真是uint64的数据？如果不需要，可以改用uint32的，还有你这个图片，看起来像是可以变成稀疏化体素的，不知道是自己写还是因为别的第三方库就要求这种XYZC的稠密格式？有没有办法改用专业的openvdb？ |
|
还有一件事，你知不知道你这个切片其实还是有一定计算的，那就是整数索引的乘法计算，因为你的数组大小都不保证是2的幂次方，不仅没有对齐保障也需要一直计算整数乘法（虽然最后还会是内存瓶颈）。使用专业的openvdb稀疏存储，那么每个儿子块都是8x8x8，可以用高效的 |
|
这是因为我们的体像素本身就是要稠密的，后续的深度学习也需要用稠密数组。
虽然你看到图像有很多黑像素，但是在 那个立方形感兴趣框里面，大部分是有效像素。说体像素，一般是说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: ***@***.***>
|
我又来打扰了。我记得小彭老师说过，连续的内存访问速度比较快。可是我的高维数组进行抽提时，低维也有很长片段的连续，为啥还是很慢。需要用什么 stream 或者 prefetch 来优化吗？
这个 4.9 GByte/s 的内存访问速度就很离谱，我还没有做任何计算。
https://github.com/chenxinfeng4/hpc_issue_array_continguous
The text was updated successfully, but these errors were encountered: