多面体编译模型深度工程实战：从依赖多面体、仿射调度到 Tiling 与 GPU 代码生成的全链路解析

多面体编译模型深度工程实战：从依赖多面体、仿射调度到 Tiling 与 GPU 代码生成的全链路解析

执行摘要

绝大多数人对编译器循环优化的理解停留在「循环交换、循环展开、向量化」这一串模式匹配式的 pass 上。但当你面对一个三层嵌套的 stencil、一个带非完美嵌套的矩阵乘、或者一次需要同时满足数据局部性与并行度的卷积时，基于语法树的启发式 pass 会迅速失效——因为每一个变换都要重新做一遍合法性检查，而变换之间还会互相破坏彼此的前提。

多面体模型（Polyhedral Model）给出的是另一条路：把「循环嵌套」翻译成几何对象（整数点集与仿射映射），把「变换是否合法」翻译成依赖多面体上的线性约束求解，把「找一个好的变换」翻译成整数线性规划（ILP）。在这一框架里，循环交换、fission、fusion、skewing、tiling、并行化不再是各自独立的 pass，而只是同一个调度系数矩阵的不同解。

本文沿着真实工程链路拆解多面体编译：迭代域与访问函数的仿射表示、依赖分析的精确化路径（GCD → Banerjee → Farkas 引理 + ILP）、Feautrier/PLuTo 的仿射调度搜索、Tiling 与 AST 代码生成，以及它在 MLIR Affine Dialect 与深度学习编译器里的工业落地形态。每一节都附带可对照的表示与代码。

一、先建立坐标：为什么循环优化需要「几何化」

传统编译器把循环表示成 AST + CFG。这种表示对控制流友好，但对数据流极不友好。考虑下面这个最简单的例子：

for (int i = 1; i < N; i++)
  for (int j = 1; j < N; j++)
    A[i][j] = A[i-1][j] + A[i][j-1];

这是典型的二维前缀和 / 波前（wavefront）计算。要让它并行，正确做法是对 i 和 j 做 skewing：令新循环变量 t = i + j，同一条反对角线上的点互不依赖，可并行。

问题在于：编译器怎么「知道」可以 skew？在 AST 上，你需要为 skew 单独写一个合法性证明；如果循环里还有一个 if (j > i) 的三角边界，或者数组下标是 A[2*i+3][j]，这个证明就得重写。组合爆炸。

多面体模型的做法是把程序换一套坐标系：

• 每个动态执行实例（statement instance）用迭代向量 x = (i, j) 表示；
• 循环边界写成一组仿射不等式 D = { x | B·x + b ≥ 0 }，这就是迭代域多面体；
• 每次内存访问写成仿射函数 f(x) = F·x + f，称为访问关系。

于是「程序的动态执行」= 多面体中的整数点集，「数据依赖」= 两个整数点之间的二元关系。优化就变成了在这个几何表示上做变换与求解，最后再把结果翻译回 AST。这个「翻译回去」的步骤叫 code generation（CLooG / isl AST 生成器）。

二、依赖分析：从 GCD 判定到 Farkas 引理

变换合法性的唯一判据是：保持所有依赖的方向。所以依赖分析必须精确。

2.1 依赖的本质

语句 S1 的迭代点 x 写 A[f(x)]，语句 S2 的迭代点 y 读 A[g(y)]。若存在 x, y 使 f(x) = g(y) 且 x 在 y 之前执行（按原程序顺序），则存在依赖，依赖向量为 d = y - x。

2.2 三级精度递进

现代多面体工具（isl、Omega+、Polly）全部走第三条路。其核心是 Farkas 引理：一个多面体非空，等价于不存在一组非负线性组合能导出矛盾。用它把「多面体非空」这一存在性命题转成一组线性不等式，再交给 ILP 求解器（PIP、GLPK、isl 自带的 Gomory cut 算法）得出精确依赖。

import islpy as isl

# 迭代域: 0 <= i < N, 0 <= j < N
ctx = isl.Context()
space = isl.Space.create_from_names(ctx, set=["i", "j"], params=["N"])
domain = (isl.BasicSet.universe(space)
          .add_constraint(isl.Constraint.ineq_from_names(space, {"i": 1}))
          .add_constraint(isl.Constraint.ineq_from_names(space, {"i": -1, "N": 1, 1: -1}))
          .add_constraint(isl.Constraint.ineq_from_names(space, {"j": 1}))
          .add_constraint(isl.Constraint.ineq_from_names(space, {"j": -1, "N": 1, 1: -1})))

# 依赖: A[i][j] 写, A[i-1][j] 读  →  存在 x 使 f(x)=g(y)
# f(x) = (i, j)     g(y) = (i'+1, j')
# 依赖关系: i = i'+1, j = j'  且 (i',j') 在 (i,j) 之后
# 依赖向量 d = (1, 0)

上面这段用 islpy 手工构造的过程，在真正的编译器里由前端自动完成：静态控制部分（SCoP）提取器扫描 CFG，把不含 while、不含不可预测分支、数组下标全为循环变量与参数仿射组合的循环区域抽出来，逐个构造迭代域与访问关系。

2.3 SCoP：多面体模型的适用边界

SCoP（Static Control Part）是这套方法的硬边界。一个区域能进入多面体框架，必须满足：

1. 循环边界是循环变量与外层参数的仿射函数；
2. 数组下标同上（不能出现 A[B[i]] 这类间接寻址）；
3. 无 while、无 break/goto、无条件依赖数据的分支（可保留 if，但条件需仿射并可转为 guard 约束）；
4. 所有函数调用可内联或无副作用注解。

现实中的数值代码（stencil、BLAS、卷积、物理仿真）大多满足；通用业务代码大多不满足。这也解释了为什么多面体编译在 HPC 与 AI 编译器里是主力，而在通用编译器（GCC/LLVM 主线）里始终只是可选 pass。

三、仿射调度：把「变换」变成「解方程」

有了依赖多面体，接下来要为每个 statement 的每个迭代点 x 计算一个多维时间戳 θ(x) = (θ1(x), θ2(x), ...)，使得新执行顺序依然满足所有依赖。这就是调度（schedule）。

3.1 合法性的线性表达

若存在依赖 S1(x) → S2(y)，则必须有 θ_{S1}(x) ≺ θ_{S2}(y)（字典序严格小于）。当 θ 被限制为仿射函数 θ_S(x) = T_S · x + c_S 时，这个约束对系数 T_S 是完全线性的。于是：

• 所有合法变换的集合 = 一个关于 T 的线性不等式系统的解集（一个多面体）；
• 循环交换、skew、shift、fusion 全部只是这个解集中的不同点。

3.2 目标函数：PLuTo 的做法

光合法不够，还要「好」。PLuTo 的经典做法是最小化依赖距离：对每个依赖，最小化 T·d 在字典序下的第一个非零分量。这直接对应「让 producer 和 consumer 在时间上尽量靠近」，即提升数据局部性；同时把强连通分量内部的依赖距离设为 0，从而暴露出可并行的外层循环。

# 概念性伪代码：PLuTo 的 ILP 目标
# 变量: 每个 statement S 的调度系数 T_S (行向量 per dimension)
# 约束: 对每条依赖 (S1,x) -> (S2,y), 依赖向量 d:
#       T_S2 · (x + d) - T_S1 · x >= 0     (字典序, 由 Farkas 线性化)
# 目标: 字典序最小化  sum_over_dependencies ( T · d )
#       —— 按依赖的"承载面"分组, 逐层贪心求解

实际实现中，PLuTo 会把依赖按其可达基（carried basis）分组，从最内层维开始逐维求解 ILP，每求出一维就固定它再求下一维。这是把「找最优多维调度」这个 NP-hard 问题变成可工程化的贪心分层求解的关键工程妥协。

3.3 一个具体结果

回到开篇的波前例子。原程序依赖向量集是 {(1,0), (0,1)}，强连通（SCC 包含 S1 自身），传统并行化分析会判定「不可并行」。但多面体调度求出的解是：

θ(S(i,j)) = (i + j,  i)

第一维 t = i + j 上依赖距离恒为 0（因为 (1,0) 与 (0,1) 在 i+j 上的投影都是 1... 注意这里依赖方向为 -），于是新循环的最外层 t 是完全并行的，内层 i 保持串行。这正是手写的 skew 变换，只不过它是解出来的，不是匹配出来的。

四、Tiling 与代码生成

调度只给出顺序，要真正拿到性能还需要 Tiling（分块）：把迭代空间切成小方块，让每块的数据能装进 L1/L2 或 shared memory。

4.1 Tiling 的几何解释

Tiling 本质是在调度上再叠加一层变换：把 θ 拆成 (tile_id, intra_tile_id)。若 tile 大小为 b，则 θ_tile = floor(θ / b)，θ_point = θ mod b。这两层各自成为一层循环，外层遍历 tile，内层遍历块内点。

关键工程点：Tiling 的合法性必须由 tile 间的依赖方向保证。若某个依赖在某个维度上距离为负，则该维度不可分块（FCO——Full Column Rank 条件不满足时需先做 skewing 修正）。这也是为什么「先调度、后分块」是多面体工具链的固定顺序。

// isl AST 生成器输出（示意）：波前 + tiling 后的代码
for (int t = 0; t < 2*N-1; t += 32)          // tile 外层，可并行
  for (int ii = max(1, t); ii < min(N, t+32); ii += 32)
    for (int tt = t; tt < min(t+32, 2*N-1); tt++)      // 波前维
      #pragma omp parallel for
      for (int i = max(1, tt-N+1, ii); i <= min(N-1, tt-1, ii+31); i++) {
        int j = tt - i;
        A[i][j] = A[i-1][j] + A[i][j-1];
      }

注意生成代码里那些 max/min——它们不是人写的，而是 isl 对多面体做整数投影（existential elimination）后重新构造出的边界不等式。这是多面体代码生成最"魔法"也最容易出 bug 的地方：投影算法的复杂度会随维数指数上升，是超大循环嵌套编译时间爆炸的主因。

4.2 从调度树到 GPU

在 PPCG / MLIR GPU 路径上，多面体调度还能直接映射硬件层级：

同样的调度框架，只改目标函数与映射规则，就能从 OpenMP 代码变成 CUDA 代码。这是多面体模型相对传统 pass 最本质的工程优势。

五、工业落地形态与工程边界

5.1 三条主要落地路径

1. Polly（LLVM）：作为 LLVM 的 opt pass 存在，抽取 SCoP → isl 调度 → 生成 OpenMP/向量化代码 → 回填 LLVM IR。默认在 -O3 的某些发行版里开启，但对非数值代码几乎无收益。
2. PPCG / Pluto 独立工具：面向 C 数值内核，产出 CUDA/OpenCL，在 stencil 与稠密线性代数上收益显著。
3. MLIR Affine Dialect：把多面体表示作为一等公民放进 IR 层——affine.for、affine.if、affine.load/store 直接携带仿射约束，调度变换由 dialect 转换完成，最后 lowering 到 SCF + Vector。TVM、IREE、Polygeist 都在这条路径上。

5.2 必须清醒的三个工程代价

• 编译时间：ILP 与投影算法复杂度高。实测中，六层以上嵌套、或迭代域约束超过几十条时，isl 求解可能从毫秒级劣化到分钟级。生产里必须设 SCoP 规模上限（语句数、维数、约束数）并做 fallback。
• 数值参数的黑盒：当 array size 是运行时值时，多面体分析全部走参数化（parametric）模式，约束求解难度上一个量级。常见做法是生成多版本代码 + 运行时 dispatch（versioning），对 N 较小或整除性特殊的情形走专门路径。
• 与手写 intrinsics 的冲突：多面体生成的是「干净的仿射循环」，一旦你手工写了 SIMD intrinsic 或 register blocking，SCoP 提取就断了。二者是竞争而非叠加关系。

5.3 生产建议

对数值密集型团队，我的判断是：

1. 不要指望通用编译器自动生效。GCC/LLVM 的多面体 pass 覆盖面有限，收益常常为 0。
2. 把多面体用在「内核生成器」而非「编译器」上——即为 stencil、卷积、matmul 这类模式化内核单独跑 PPCG/MLIR 生成代码，作为库交付，而不是让全量代码过一遍多面体 pass。
3. 先保证 SCoP 的可提取性：这是最大的工程杠杆。把间接寻址、while 循环、数据相关分支从热点内核里剥离出去，收益远大于调任何 isl 参数。
4. 把调度结果固化：ILP 求解有不稳定风险，一旦得到好的调度，用 schedule 文件（isl 可序列化）固化下来，避免编译器版本升级导致性能回退。

六、结论

多面体编译的核心思想可以用一句话概括：把「程序变换的合法性」从语法层面的模式匹配，提升为几何层面的约束求解。

这条链路串起来是：

1. SCoP 提取把可分析的循环区域从 CFG 中隔离出来，代价是严格的适用边界；
2. 仿射表示把迭代域、访问函数、依赖全部写成线性不等式，代价是只能处理仿射下标；
3. Farkas + ILP给出精确依赖与合法调度的完整解空间，代价是 NP-hard 需用贪心分层逼近；
4. 调度树 → AST 生成把几何解翻译回循环代码，代价是投影算法的复杂度爆炸；
5. 层级映射把调度维度映射到 OpenMP/CUDA/warp，代价是必须与手写优化二选一。

真正值得带走的工程原则是：当优化的搜索空间可以用数学对象完整刻画时，就该用求解器去搜索，而不是写一千条启发式规则。多面体模型、寄存器分配的图着色、指令调度的约束规划，本质上是同一个思路在不同坐标上的投影——理解了这一点，再看 MLIR 的 dialect 体系、TVM 的 schedule 原语，会发现它们都是同一套思想在编译器工程里的具体落地。