2

我正在尝试使用 C++ AMP 优化算法(Lattice Boltzmann)以进行并行计算。并寻找一些优化内存布局的建议,刚刚发现从结构中删除一个参数到另一个向量(阻塞向量)会增加大约 10%。

任何人都有任何可以进一步改善这一点的提示,或者我应该考虑什么?下面是每个时间步执行的最耗时的函数,以及用于布局的结构。

struct grid_cell {
//  int blocked;    // Define if blocked
    float n;        // North
    float ne;       // North-East
    float e;        // East
    float se;       // South-East
    float s;
    float sw;
    float w;
    float nw;
    float c;        // Center
};

int collision(const struct st_parameters param, vector<struct grid_cell> &node, vector<struct grid_cell> &tmp_node, vector<int> &obstacle) {
    int x,y;
    int i = 0;
    float c_sq = 1.0f/3.0f;     // Square of speed of sound
    float w0 = 4.0f/9.0f;       // Weighting factors
    float w1 = 1.0f/9.0f;
    float w2 = 1.0f/36.0f;

    int chunk = param.ny/20;
    float total_density = 0;

    float u_x,u_y;              // Avrage velocities in x and y direction
    float u[9];                 // Directional velocities
    float d_equ[9];             // Equalibrium densities
    float u_sq;                 // Squared velocity
    float local_density;        // Sum of densities in a particular node

    for(y=0;y<param.ny;y++) {
        for(x=0;x<param.nx;x++) {
            i = y*param.nx + x; // Node index
            // Dont consider blocked cells
            if (obstacle[i] == 0) {
                // Calculate local density
                local_density = 0.0;
                local_density += tmp_node[i].n;
                local_density += tmp_node[i].e;
                local_density += tmp_node[i].s;
                local_density += tmp_node[i].w;
                local_density += tmp_node[i].ne;
                local_density += tmp_node[i].se;
                local_density += tmp_node[i].sw;
                local_density += tmp_node[i].nw;
                local_density += tmp_node[i].c;

                // Calculate x velocity component
                u_x = (tmp_node[i].e + tmp_node[i].ne + tmp_node[i].se - 
                      (tmp_node[i].w + tmp_node[i].nw + tmp_node[i].sw)) 
                      / local_density;
                // Calculate y velocity component
                u_y = (tmp_node[i].n + tmp_node[i].ne + tmp_node[i].nw - 
                      (tmp_node[i].s + tmp_node[i].sw + tmp_node[i].se)) 
                      / local_density;
                // Velocity squared
                u_sq = u_x*u_x +u_y*u_y;

                // Directional velocity components;
                u[1] =  u_x;        // East
                u[2] =        u_y;  // North
                u[3] = -u_x;        // West
                u[4] =      - u_y;  // South
                u[5] =  u_x + u_y;  // North-East
                u[6] = -u_x + u_y;  // North-West
                u[7] = -u_x - u_y;  // South-West
                u[8] =  u_x - u_y;  // South-East

                // Equalibrium densities
                // Zero velocity density: weight w0
                d_equ[0] = w0 * local_density * (1.0f - u_sq / (2.0f * c_sq));
                // Axis speeds: weight w1
                d_equ[1] = w1 * local_density * (1.0f + u[1] / c_sq
                                 + (u[1] * u[1]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[2] = w1 * local_density * (1.0f + u[2] / c_sq
                                 + (u[2] * u[2]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[3] = w1 * local_density * (1.0f + u[3] / c_sq
                                 + (u[3] * u[3]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[4] = w1 * local_density * (1.0f + u[4] / c_sq
                                 + (u[4] * u[4]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                // Diagonal speeds: weight w2
                d_equ[5] = w2 * local_density * (1.0f + u[5] / c_sq
                                 + (u[5] * u[5]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[6] = w2 * local_density * (1.0f + u[6] / c_sq
                                 + (u[6] * u[6]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[7] = w2 * local_density * (1.0f + u[7] / c_sq
                                 + (u[7] * u[7]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[8] = w2 * local_density * (1.0f + u[8] / c_sq
                                 + (u[8] * u[8]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));

                // Relaxation step
                node[i].c = (tmp_node[i].c + param.omega * (d_equ[0] - tmp_node[i].c));
                node[i].e = (tmp_node[i].e + param.omega * (d_equ[1] - tmp_node[i].e));
                node[i].n = (tmp_node[i].n + param.omega * (d_equ[2] - tmp_node[i].n));
                node[i].w = (tmp_node[i].w + param.omega * (d_equ[3] - tmp_node[i].w));
                node[i].s = (tmp_node[i].s + param.omega * (d_equ[4] - tmp_node[i].s));
                node[i].ne = (tmp_node[i].ne + param.omega * (d_equ[5] - tmp_node[i].ne));
                node[i].nw = (tmp_node[i].nw + param.omega * (d_equ[6] - tmp_node[i].nw));
                node[i].sw = (tmp_node[i].sw + param.omega * (d_equ[7] - tmp_node[i].sw));
                node[i].se = (tmp_node[i].se + param.omega * (d_equ[8] - tmp_node[i].se));
            }
        }
    }
    return 1;
}
4

4 回答 4

7

通常,您应该确保在不同 cpu 上使用的数据不共享(容易)并且不在同一缓存行上(错误共享,请参见此处的示例:错误共享是没有乐趣的)。同一个 cpu 使用的数据应该靠近在一起才能从缓存中受益。

于 2012-03-13T13:24:51.400 回答
6

众所周知,当前的 GPU 依赖于内存布局。如果没有关于您的应用程序的更多详细信息,我建议您探索以下内容:

  1. 单位步长访问非常重要,因此 GPU 更喜欢“数组结构”而不是“结构数组”。正如您将“blocked”字段移动到向量“obstacle”一样,将“grid_cell”的所有字段转换为单独的向量应该是有利的。对于不访问所有字段的循环,这应该显示对 CPU 的好处。

  2. 如果“障碍物”非常稀疏(我猜这不太可能),那么将其移动到自己的向量中就特别有价值。像 CPU 这样的 GPU 将以缓存行或其他形式从内存系统加载多个字,当您不需要某些数据时,您会浪费带宽。对于许多系统而言,内存带宽是瓶颈资源,因此任何减少带宽的方法都有帮助。

  3. 这更具推测性,但现在您正在写入所有输出向量,内存子系统可能正在避免读取“节点”中的值,这些值将被简单地覆盖

  4. 在某些系统上,片上存储器被拆分为组,并且在您的结构中具有奇数个字段可能有助于消除组冲突。

  5. 一些系统还将“矢量化”加载和存储,因此再次从结构中删除“阻塞”可能会启用更多矢量化。向数组结构的转变减轻了这种担忧。

感谢您对 C++ AMP 的关注。

大卫卡拉汉

http://blogs.msdn.com/b/nativeconcurrency/ C++ AMP 团队博客

于 2012-03-16T18:08:21.953 回答
1

一些小的通用上衣:

  • 跨多个处理器共享的任何数据结构都应该是只读的。

  • 任何需要修改的数据结构对于处理器来说都是唯一的,并且不会与另一个处理器所需的数据共享内存位置。

  • 确保您的内存被安排好,以便您的代码可以连续扫描它(不要走大步或跳来跳去)。

于 2012-03-13T15:47:22.613 回答
0

对于任何研究这个主题的人来说,一些提示。Lattice-Boltzmann 通常是带宽受限的。这意味着它的性能主要取决于可以从内存加载和写入内存的数据量。

  • 使用高效的编译编程语言:C 或 C++是基于 CPU 实现的不错选择。

  • 选择具有高带宽的架构。对于 CPU,这意味着高时钟 RAM 和大量内存通道(四通道或更多)。

  • 这使得使用有效利用缓存预取的适当线性内存布局变得至关重要:数据在内存中以小部分排列,即所谓的缓存行。每当处理器访问一个元素时,它所在的整个缓存行(在现代架构上为 64 字节)都会被加载。这意味着一次加载 8 个双精度或 16 个浮点值!虽然我没有发现这对于多核处理器来说是个问题,因为它们共享 L3 缓存,但这应该会导致多核架构出现问题,因为对同一缓存行的更改必须保持同步,并且当其他处理器出现问题时正在处理另一个处理器正在处理的数据(错误共享)。这可以通过引入填充来避免,这意味着您添加不会用于填充其余缓存行的元素。假设您要更新具有 D3Q27 格的 27 速度的离散化单元,那么在双精度数(8 字节)的情况下,数据位于 4 个不同的缓存行上。您应该添加 5 双填充以匹配 32 字节(4*8 字节)。

unsigned int const PAD          = (64 - sizeof(double)*D3Q27.SPEEDS % 64); ///< padding: number of doubles
size_t       const MEM_SIZE_POP = sizeof(double)*NZ*NY*NX*(SPEEDS+PAD);    ///< amount of memory to be allocated

大多数编译器自然地将数组的开头与缓存行对齐,因此您不必照顾它。

  • 线性索引不方便访问。因此,您应该尽可能高效地设计访问。你可以写一个包装类。在任何情况下内联这些函数,这意味着每个调用都被它们在代码中的定义替换。
inline size_t const D3Q27_PopIndex(unsigned int const x, unsigned int const y, unsigned int const z, unsigned int const d)
{
    return (D3Q27.SPEEDS + D3Q27.PAD)*(NX*(NY*z + y) + x) + D3Q27.SPEEDS*p + d;
}
  • 此外,可以通过最大化计算和通信之间的比率来增加缓存局部性,例如使用三维空间循环阻塞OpenMP 的缩放问题),这意味着每个代码都在一个单元格而不是单个单元格上工作。

  • 通常,实现使用两个不同的种群 A 和 B,并执行从一个实现到另一个实现的冲突和流式传输。这意味着内存中的每个值都存在两次,一次是碰撞前,一次是碰撞后。存在用于重组步骤和存储的不同策略,您只需在内存中保留一个总体副本。例如,参见P. Bailey 等人提出的AA 模式。-“使用图形处理器加速格子玻尔兹曼流体流动模拟”(https://www2.cs.arizona.edu/people/pbailey/Accelerating_GPU_LBM.pdf)或深奥的扭曲作者:M. Geier 和 M. Schönherr -“Esoteric Twist: An Efficient in-Place Streaming Algorithmus for the Lattice Boltzmann Method on the Massively Parallel Hardware”(https://pdfs.semanticscholar.org/ea64/3d63667900b60e6ff49f2746211700e63802.pdf)。我已经使用宏实现了第一个,这意味着人口的每次访问都调用以下形式的宏:

#define O_E(a,b)       a*odd + b*(!odd)

#define  READ_f_0      D3Q27_PopIndex(x,           y, z,           0, p)
#define  READ_f_1      D3Q27_PopIndex(O_E(x_m, x), y, z, O_E( 1,  2), p)
#define  READ_f_2      D3Q27_PopIndex(O_E(x_p, x), y, z, O_E( 2,  1), p)
...

#define  WRITE_f_0     D3Q27_PopIndex(x,           y, z,           0, p)
#define  WRITE_f_1     D3Q27_PopIndex(O_E(x_p, x), y, z, O_E( 2,  1), p)
#define  WRITE_f_2     D3Q27_PopIndex(O_E(x_m, x), y, z, O_E( 1,  2), p)
...
  • 如果您有多个交互群体,请使用网格合并。在内存中线性排列索引并将两个不同的群体并排放置。人口 p 的访问工作如下:
inline size_t const D3Q27_PopIndex(unsigned int const x, unsigned int const y, unsigned int const z, unsigned int const d, unsigned int const p = 0)
{
    return (D3Q27.SPEEDS*D3Q27.NPOP + D3Q27.PAD)*(NX*(NY*z + y) + x) + D3Q27.SPEEDS*p + d;
}
  • 对于规则网格,使算法尽可能可预测。让每个单元执行碰撞和流式处理,然后反向执行边界条件。如果您有许多对算法没有直接贡献的单元格,请使用您也可以存储在填充中的逻辑掩码省略它们!

  • 在编译时让编译器知道所有内容:例如边界条件模板,其中包含一个处理索引更改的函数,因此您不必重写每个边界条件。

  • 现代架构具有可以执行SIMD操作的寄存器,因此对多个数据执行相同的指令。一些处理器 (AVX-512) 最多可以处理 512 位数据,因此 32 位的速度几乎是单个数字的两倍。这似乎对 LBM 非常有吸引力,特别是自从引入了收集和散射(https://en.wikipedia.org/wiki/Gather-scatter_(vector_addressing))但由于当前的带宽限制(也许值得使用 DDR5 和内核很少的处理器),我认为这不值得麻烦:单核性能和并行扩展更好(M. Wittmann 等人 - “Lattice Boltzmann Benchmark Kernels as a Testbed for Performance Analysis” - https: //arxiv.org/abs/1711.11468) 但由于带宽有限,整体算法的性能并没有更好。因此,它仅适用于受计算能力而非带宽限制的架构。在 Xeon Phi 架构上,结果似乎是显着的 Robertsen 等人。-“至强融核处理器上格子玻尔兹曼方法的高性能 SIMD 实现”(https://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.5072)。

在我看来,对于简单的 2D 实现来说,大部分这些都不值得。在那里进行简单的优化,循环阻塞,线性内存布局,但忘记更复杂的访问模式。在 3D 中,效果可能是巨大的:我在现代 12 核处理器上使用 D3Q19 实现了高达 95% 的并行可扩展性和超过 150 Mlups 的整体性能。要获得更高的性能,请切换到更合适的架构,例如针对带宽进行了优化的具有CUDA C的GPU

于 2019-10-20T12:19:50.053 回答