只是打乱单元编号，LU-SGS 就慢了 2.8 倍 — RCM 重排序与决定扫描方向的起始单元

在 LU-SGS 中，单元编号决定信息在一次扫描中能传多远。RCM 把带宽降到 24，但区分 9 次迭代和 21 次迭代的，是起始单元与流动方向的关系。

打乱单元编号，会改变什么？

同一套网格。只把单元编号随机打乱。矩阵不过是行列换序。收敛解不动。可是 LU-SGS（Lower-Upper Symmetric Gauss–Seidel，用一次下三角扫描和一次上三角扫描近似求解隐式系统的方法）的迭代从 9 次变成了 25 次。差距在于一次扫描里信息能走多远。576 个单元的输出把 RCM（Reverse Cuthill–McKee，用广度优先搜索缩小矩阵带宽的重新编号方法）修得了的和修不了的分开。

2017 年笔记里的问题：哪个邻居是 Lower？#

2017 年底给非结构 3D 代码加入 LU-SGS 时，留下过一段笔记。大意如下。LU-SGS 把计算分成 Lower 部分和 Upper 部分。但在非结构网格上，很难判断周围单元哪个是 Lower、哪个是 Upper。于是找来一份 RCM 开源代码，在初始化阶段对单元重新排序。

现在回头看，这个问题的前半句只对了一半。Lower 和 Upper 不由网格决定，而由单元的编号决定。编号比自己小的邻居是 Lower，比自己大的是 Upper。只要有编号，LU-SGS 就能运行。真正的问题是"哪种编号更好"。

笔记里还有一句："据说这份 RCM 代码只适用于四面体网格，但我觉得其他网格也能用。"这句是对的。RCM 不看单元形状。它只看共享面的邻居列表，也就是图。

在 random labels 和 RCM 之间切换，20×20 个单元的一阶迎风对流–扩散里，颜色铺开的距离会变。

RCM start:

flow angle = 30° diffusion D = 0.010 (cell Pe = 100) cells/frame = 16

sweeping| iter = 0 | residual = 1.0e+0natural: 9 it, bw 20random: 21 it, bw 392RCM SW: 9 it, bw 20RCM SE: 18 it, bw 20RCM NW: 18 it, bw 20RCM NE: 9 it, bw 20

随机编号时，一次前向扫描过后颜色仍是零星的。RCM SW 顺着流动箭头，一次扫描就几乎铺满。把 RCM start 换成 SE，bw 不动，迭代次数变长。

把 LU-SGS 的两次扫描写成公式#

设用伪时间（pseudo-time）推进求解定常线性系统 $A\phi = b$ 。把矩阵拆成对角和非对角部分。

A = D - L - U, \qquad (D + \tfrac{1}{\Delta\tau}I - L)\,(D + \tfrac{1}{\Delta\tau}I)^{-1}\,(D + \tfrac{1}{\Delta\tau}I - U)\,\Delta\phi = r

$D$ 是对角元 $a_P$ ， $L$ 是编号较小邻居的系数 $a_{nb}$ ， $U$ 是编号较大邻居的系数， $r = b - A\phi$ 是残差， $\Delta\tau$ 是伪时间步长。

这个近似分解通过两次扫描求解。前向扫描按编号顺序前进，只使用编号较小邻居的新值。后向扫描按逆序前进，加上编号较大邻居的贡献。

\Delta\phi^{*}_c = \frac{r_c + \sum_{k<c} a_{ck}\,\Delta\phi^{*}_k}{d}, \qquad \Delta\phi_c = \Delta\phi^{*}_c + \frac{\sum_{k>c} a_{ck}\,\Delta\phi_k}{d}

$k<c$ 表示"顺序排在 $c$ 之前的邻居"， $d = a_P + 1/\Delta\tau$ 。

关键就在这里。如果单元 $c$ 的上游邻居排在前面，这个邻居本次扫描得到的修正会立即传给 $c$ 。如果排在后面，就要等到下一次迭代。对纯迎风对流而言，若沿流动方向编号， $L$ 就包含了整个矩阵。一次前向扫描就是精确解。

Cuthill–McKee 就是广度优先搜索#

Cuthill 和 McKee 在 1969 年发表的算法很简单。选一个起始单元。用广度优先搜索（BFS，从近到远逐层访问邻居）编号。同一层内，邻居数少的单元先编号。这样，相邻单元之间的编号差不会超过一层的宽度。

这个编号差的最大值就是带宽。

\beta = \max_{(i,j)\,\in\,\text{faces}} \lvert p(i) - p(j) \rvert

$p(i)$ 是单元 $i$ 的新编号，最大值取遍所有共享面的单元对。

1971 年，George 证明只要把这个顺序倒过来，直接分解中的填充（fill-in）就会减少。这就是 RCM 里的 R。起始单元通常取在图的"边缘"。标准做法是伪外围节点（pseudo-peripheral node）搜索：反复做 BFS，跳到最远的单元。

代码：576 个单元，四种编号#

在 24×24 网格上设置两个问题。一个是 30° 方向流动加小扩散的对流主导问题（单元 Péclet 数 100）。另一个是纯扩散。为了模拟网格生成器随意编号的情况，先打乱标签，再用 RCM 重新编号。对流问题中，还把 RCM 起始单元换到四个角分别试一遍。

import math, random
from collections import deque
 
def build_grid(n, ux, uy, diff):
    # 单元 (i,j) -> 邻居列表、系数 a_nb、对角元 a_P（一阶迎风 + 扩散）
    nbr, coef = [], []
    for j in range(n):
        for i in range(n):
            lst, cf = [], []
            for di, dj, nx_, ny_ in ((1,0,1,0),(-1,0,-1,0),(0,1,0,1),(0,-1,0,-1)):
                ii, jj = i+di, j+dj
                if 0 <= ii < n and 0 <= jj < n:
                    lst.append(jj*n+ii)
                    cf.append(diff + max(-(ux*nx_+uy*ny_), 0.0))
            nbr.append(lst); coef.append(cf)
    ap = 4*diff + abs(ux) + abs(uy)
    return nbr, coef, ap
 
def shuffle_labels(nbr, coef, rng):
    m = len(nbr)
    perm = list(range(m)); rng.shuffle(perm)   # old -> new
    inv = [0]*m
    for o, nw in enumerate(perm): inv[nw] = o
    nb2 = [[int(perm[k]) for k in nbr[inv[c]]] for c in range(m)]
    cf2 = [coef[inv[c]] for c in range(m)]
    return nb2, cf2, perm
 
def bfs_levels(nbr, s):
    dist = [-1]*len(nbr); dist[s] = 0; q = deque([s]); last = s
    while q:
        c = q.popleft(); last = c
        for k in nbr[c]:
            if dist[k] < 0:
                dist[k] = dist[c]+1; q.append(k)
    return dist, last
 
def peripheral_start(nbr):
    s, ecc = 0, -1
    while True:
        dist, far = bfs_levels(nbr, s)
        if dist[far] <= ecc: return s
        ecc, s = dist[far], far
 
def cuthill_mckee(nbr, start):
    seen = [False]*len(nbr); seen[start] = True
    order, q = [], deque([start])
    while q:
        c = q.popleft(); order.append(c)
        for k in sorted(nbr[c], key=lambda k: len(nbr[k])):
            if not seen[k]:
                seen[k] = True; q.append(k)
    return order
 
def bandwidth(nbr, order):
    pos = {c: p for p, c in enumerate(order)}
    return max(abs(pos[c]-pos[k]) for c in range(len(nbr)) for k in nbr[c])
 
def lusgs_solve(nbr, coef, ap, order, dtau, tol=1e-8, itmax=4000):
    m = len(nbr); x = [0.0]*m; b = [1.0]*m
    d = ap + 1.0/dtau
    rank = [0]*m
    for p, c in enumerate(order): rank[c] = p
    bnorm = m**0.5
    for it in range(1, itmax+1):
        r = [b[c] - ap*x[c] + sum(a*x[k] for a, k in zip(coef[c], nbr[c])) for c in range(m)]
        if sum(v*v for v in r)**0.5/bnorm < tol: return it-1
        dx = [0.0]*m
        for c in order:                      # forward: 只用 rank 较小的邻居 (L)
            s = r[c] + sum(a*dx[k] for a, k in zip(coef[c], nbr[c]) if rank[k] < rank[c])
            dx[c] = s/d
        for c in reversed(order):            # backward: 只用 rank 较大的邻居 (U)
            s = sum(a*dx[k] for a, k in zip(coef[c], nbr[c]) if rank[k] > rank[c])
            dx[c] += s/d
        for c in range(m): x[c] += dx[c]
    return itmax
 
rng = random.Random(7)
n = 24
th = math.radians(30)
for name, ux, uy, df, dtau in (("convection", math.cos(th), math.sin(th), 0.01, 1e3),
                               ("diffusion", 0.0, 0.0, 1.0, 0.5)):
    nb0, cf0, ap = build_grid(n, ux, uy, df)
    nb, cf, perm = shuffle_labels(nb0, cf0, rng)
    m = len(nb)
    s = peripheral_start(nb)
    cm = cuthill_mckee(nb, s)
    rcm = cm[::-1]
    # 自然顺序直接使用打乱前的网格
    orders = {
        "natural (i,j)": (nb0, cf0, list(range(m))),
        "random labels": (nb, cf, list(range(m))),
        "Cuthill-McKee": (nb, cf, cm),
        "reverse CM": (nb, cf, rcm),
    }
    print(f"[{name}] cells={m}")
    for k, (g, c, o) in orders.items():
        it = lusgs_solve(g, c, ap, o, dtau)
        print(f"  {k:15s} bandwidth={bandwidth(g, o):5d}  LU-SGS iters={it:5d}")
    if name == "convection":
        for corner, old in (("SW", 0), ("SE", n-1), ("NW", n*n-n), ("NE", n*n-1)):
            o = cuthill_mckee(nb, perm[old])[::-1]
            print(f"  RCM start {corner}    bandwidth={bandwidth(nb, o):5d}  LU-SGS iters={lusgs_solve(nb, cf, ap, o, dtau):5d}")

[convection] cells=576
  natural (i,j)   bandwidth=   24  LU-SGS iters=    9
  random labels   bandwidth=  553  LU-SGS iters=   25
  Cuthill-McKee   bandwidth=   24  LU-SGS iters=    9
  reverse CM      bandwidth=   24  LU-SGS iters=    9
  RCM start SW    bandwidth=   24  LU-SGS iters=    9
  RCM start SE    bandwidth=   24  LU-SGS iters=   21
  RCM start NW    bandwidth=   24  LU-SGS iters=   21
  RCM start NE    bandwidth=   24  LU-SGS iters=    9
[diffusion] cells=576
  natural (i,j)   bandwidth=   24  LU-SGS iters= 1549
  random labels   bandwidth=  561  LU-SGS iters= 1747
  Cuthill-McKee   bandwidth=   24  LU-SGS iters= 1549
  reverse CM      bandwidth=   24  LU-SGS iters= 1549

同样 576 个单元，同样的矩阵，同样的收敛标准（ $10^{-8}$ ）。改变的只有访问单元的顺序。

区分 25 次和 9 次的是上游邻居的顺序#

流动方向为 30° 时，每个单元的上游邻居在西侧和南侧。自然顺序（先 $i$ 后 $j$ ）下，两者都排在前面。因此一次前向扫描就能把信息从西南角一直送到东北角。剩下的 9 次迭代用来消化小的扩散项和伪时间项。

随机编号时，上游邻居排在前面的概率只有一半。信息在一次扫描中走几个单元就断了。要穿过 24 个单元，需要多次迭代。结果是 25 次，约 2.8 倍。

带宽也从 553 降到 24。这个数字和迭代次数是两回事。在实际代码中，它表示邻居单元的数据在内存中相隔多远。带宽小，缓存命中率高，每次迭代就更快。所以 RCM 同时减少了迭代次数和每次迭代的开销。不过 Python 实验测不到缓存效应，这里只比较迭代次数。

同样带宽 24，迭代却是 9 次和 21 次#

四个角的实验更有意思。无论从哪个角出发，带宽都是 24。但 SW 和 NE 是 9 次，SE 和 NW 是 21 次。

RCM 的 BFS 层是以起始角为中心的斜向条带。从 SW 出发，条带与流动垂直，顺序沿着流动增加。从 NE 出发，只是顺序逆着流动增加。这时后向扫描顺着流动方向。SGS 前后两次扫描都做，所以两种情况是对称的。同理，Cuthill–McKee 和 reverse CM 的迭代次数也相同。

从 SE 出发，顺序朝西北方向增加。这与流动的 $x$ 分量相反，与 $y$ 分量相同。无论前向还是后向扫描，都只有一个上游邻居能及时传到。所以是 21 次。带宽只看矩阵的形状，不知道流动方向。

25 除以 9 是 2.78。标题里的 2.8 倍就是这个四舍五入。同样 576 个单元，同样 $10^{-8}$ 。带宽从 553 降到 24。553 贴着单元数 576。邻居编号几乎横跨整张矩阵。24 是自然顺序已经碰到的 $\beta$ 下限。24×24 上，纵向邻居在 $j$ 增加 1 时，编号差正好是 $n=24$ 。RCM 收到了这个下限，再往下没有了。

四个角的 $\beta$ 都是 24。迭代只有 SW 9、NE 9、SE 21、NW 21。21/9 是 2.33，比随机的 25 好，仍比顺着流动的 9 多一倍以上。SE 和 NW 同为 21，是因为两次扫描方向相反。一条对角线上出 9，另一条对角线也是 9。这和 CM 与 reverse CM 都是 9 是同一种对称。

扩散问题只差 13%#

纯扩散的结果不同。自然顺序和 RCM 都是 1549 次，随机编号是 1747 次。相差 13%。

扩散没有上游。信息向各个方向以同样速度扩散。这里 Gauss–Seidel 慢，是因为长波，不是因为编号。这个长波靠排序缩不掉，得交给Krylov 求解器。编号能把迭代次数拉开，是在对流主导的时候。像可压缩 Euler 这样双曲性很强的问题，正是这种情况。

1549 和 1747 相差 198 次。198/1549 是 0.128，也就是 13%。带宽拉到 561，迭代也只动这么多。这和对流里 553 引出 25 次正好相反。把 $\beta$ 缩小，和把上游邻居放进 $k \lt c$ ，是同一次重排序里的两件事。扩散没有上游，第二件事几乎是空的。

那为什么要 "Reverse"？#

实验中 CM 和 RCM 的迭代次数相同。倒序看起来没有意义。起作用的地方在别处。

George 证明的是直接分解中的填充。用 Cholesky 或 LU 分解带状矩阵时，每一行从第一个非零元到对角元都会被填满。这些跨度之和称为"轮廓（profile）"，RCM 几乎总是不大于 CM。所以如果打算用 ILU 预条件或直接解法，RCM 是默认选择。如果只用 LU-SGS，CM 和 RCM 没有区别。对称的两次扫描抵消了顺序的方向。

1D 上从左到右编号就结束了。编号问题出在 2D、3D 非结构网格。网格生成器按单元生成的先后编号。这个顺序跟着点的插入或阵面的推进走。从求解器看，它接近随机编号。

带宽已经是 24 却停在 21 次，就丢掉起始单元#

重排序只在初始化做一次。计算中没有理由再改编号。由面的连接做成邻接图，再把单元数组、每个面的左右索引和边界列表换到新编号上。四面体还是三棱柱，看的都是图。代价是一次 BFS，与单元数成正比。

对流强的时候，起始单元比带宽更靠前。这次输出里四个角的 $\beta$ 都是 24，SW 和 NE 却是 9 次，SE 和 NW 是 21 次。伪外围单元只是图上最远的单元，不知道入口在哪。9 次是起始恰好落在流动轴线上的角。没有落在那里的保证。重排序前后解在机器精度上变了，就是漏了索引转换。编号不该改变解。1D 从左边编号就结束的分解，在Burgers 算例里的 LU-SGS里。编号问题从这里开始。