如何将菱形分块（Diamond Tiling）推广至高维空间？

核心结论 #

菱形分块完全可以推广到2D/3D空间，核心思路是基于依赖图的低维投影而非笛卡尔积——把高维的计算依赖映射到低维的"波前"维度，再在每个波前内划分菱形分块，从根源上避免梯形分块的冗余内存访问问题。

2D场景手动实现分步流程（以5点模板为例） #

假设2D PDE的5点模板计算逻辑为：u[i][j] = (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1])/4 + f[i][j]，计算域为[0, N-1] x [0, M-1]，边界条件已预处理完成。

步骤1：定义低维投影的波前维度 #

2D计算的依赖是上下左右相邻点，我们用对角线波前作为低维投影维度：定义波前索引k = i + j，每个波前包含所有满足i+j=k的网格点。

• 波前范围：k从0到(N-1)+(M-1) = N+M-2
• 每个波前内的i取值范围：i ∈ [max(0, k-(M-1)), min(N-1, k)]，对应j = k - i

步骤2：在波前内划分菱形分块 #

对每个波前k，将其包含的点划分为大小为B的菱形块（块内点数量不超过B）：

• 块起始i索引：i_start = max(0, k-(M-1))，每次步进B
• 块内i范围：i ∈ [i_start, min(i_start+B-1, min(N-1, k))]

步骤3：可直接编码的循环调度伪代码 #

// 预处理边界条件
init_boundary(u);

// 遍历所有波前
for (int k = 0; k <= N+M-2; k++) {
    // 计算当前波前的i范围
    int i_min = max(0, k - (M-1));
    int i_max = min(N-1, k);
    
    // 划分菱形块并处理
    for (int i_start = i_min; i_start <= i_max; i_start += B) {
        int i_end = min(i_start + B - 1, i_max);
        
        // 处理块内每个点
        for (int i = i_start; i <= i_end; i++) {
            int j = k - i;
            // 边界过滤（理论上已通过i_min/i_max限制，此处为双重保险）
            if (j < 0 || j >= M) continue;
            // 执行5点模板计算
            u[i][j] = (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1])/4.0 + f[i][j];
        }
    }
}

低维投影说明 #

这里的波前k=i+j是把2D计算的依赖关系投影到1D维度，同一波前内的点之间没有依赖（5点模板的依赖点不在同一对角线），因此块内的点可以并行处理；同时菱形块的划分确保每个网格点仅被加载一次，彻底消除了梯形分块的冗余内存访问。

3D场景推广思路（以7点模板为例） #

针对3D PDE的7点模板（依赖上下、左右、前后相邻点），低维投影采用空间对角线波前：k = i + j + l，每个波前包含所有满足i+j+l=k的网格点。

分步调度逻辑 #

1. 波前维度定义：k从0到(N-1)+(M-1)+(P-1) = N+M+P-3（N,M,P为三维计算域的尺寸）
2. 块划分：在每个波前k内，先固定i的范围并划分块，再在每个i对应的子空间内划分j的块，块大小设为Bx（i方向）、By（j方向）
3. 可直接编码的循环框架

// 预处理3D边界条件
init_boundary_3d(u);

// 遍历所有空间对角线波前
for (int k = 0; k <= N+M+P-3; k++) {
    // 计算当前波前的i范围
    int i_min = max(0, k - (M-1) - (P-1));
    int i_max = min(N-1, k);
    
    // 划分i方向的块
    for (int i_start = i_min; i_start <= i_max; i_start += Bx) {
        int i_end = min(i_start + Bx - 1, i_max);
        
        for (int i = i_start; i <= i_end; i++) {
            // 计算当前i对应的j范围
            int j_min = max(0, k - i - (P-1));
            int j_max = min(M-1, k - i);
            
            // 划分j方向的块
            for (int j_start = j_min; j_start <= j_max; j_start += By) {
                int j_end = min(j_start + By - 1, j_max);
                
                // 处理块内每个点
                for (int j = j_start; j <= j_end; j++) {
                    int l = k - i - j;
                    if (l < 0 || l >= P) continue;
                    // 执行7点模板计算
                    u[i][j][l] = (u[i-1][j][l] + u[i+1][j][l] + 
                                  u[i][j-1][l] + u[i][j+1][l] + 
                                  u[i][j][l-1] + u[i][j][l+1])/6.0 + f[i][j][l];
                }
            }
        }
    }
}

低维投影说明 #

3D场景中k=i+j+l将三维依赖关系投影到1D波前，同一波前内的点无依赖关系，块划分在波前的二维子空间内进行，同样保证了内存访问无冗余，完美适配HPC场景的内存带宽瓶颈。

关键调优注意事项 #

• 块大小B（或Bx,By）需根据硬件缓存大小调优，一般设置为缓存能容纳的最大块尺寸，避免缓存失效。
• 边界点需提前预处理，循环中可跳过边界计算（或初始化时直接赋值）。
• 并行化可在波前级别或块级别进行，同一波前内的块/点可并行执行。

内容的提问来源于stack exchange，提问作者比尔盖子