Shuffling Cell Numbers Made LU-SGS 2.8× Slower — RCM Reordering and the Start Cell That Sets Sweep Direction

What changes when you shuffle cell numbers?#

Same mesh. Only the cell numbers got shuffled. The matrix is a row-and-column permutation. The converged solution stays put. LU-SGS (Lower-Upper Symmetric Gauss–Seidel, which approximates an implicit system with one lower-triangular and one upper-triangular sweep) still went from 9 iterations to 25. The gap is how far information travels in one sweep. The 576-cell printout splits what RCM (Reverse Cuthill–McKee, a breadth-first renumbering that shrinks matrix bandwidth) fixes from what it cannot.

The question in a 2017 note: which neighbor is Lower?#

In late 2017 I added LU-SGS to an unstructured 3D code and left a note. The gist: LU-SGS splits the work into a Lower part and an Upper part. On an unstructured grid it is hard to tell which neighbors are Lower and which are Upper. So I found an open-source RCM routine and reordered the cells during initialization.

Looking back, the first half of that question is only half right. The grid does not define Lower and Upper. The numbering does. Neighbors with a smaller number are Lower, larger ones are Upper. Any numbering at all lets LU-SGS run. The real question is which numbering is good.

The note also says, "this RCM code supposedly only works for tetrahedral meshes, but I think it should work for others." That part is correct. RCM never looks at cell shape. It only sees the list of face neighbors, i.e. the graph.

Switching random labels and RCM changes how far the color spreads on a 20×20 first-order upwind convection–diffusion problem.

Random numbers leave only scattered patches after one forward sweep. RCM SW fills in along the flow arrow in that same sweep. Set RCM start to SE and bw stays put while the iteration count grows.

The two LU-SGS sweeps as equations#

Say we solve the steady linear system $A\phi = b$ by pseudo-time marching. Split the matrix into diagonal and off-diagonal parts.

$D$ is the diagonal $a_P$, $L$ holds the coefficients $a_{nb}$ of lower-numbered neighbors, $U$ those of higher-numbered ones, $r = b - A\phi$ is the residual, and $\Delta\tau$ is the pseudo-time step.

This approximate factorization is solved with two sweeps. The forward sweep walks in numbering order and uses only the fresh values of lower-numbered neighbors. The backward sweep walks in reverse and adds the higher-numbered contributions.

$k<c$ means "a neighbor ranked before $c$", and $d = a_P + 1/\Delta\tau$.

Here is the crux. If cell $c$'s upstream neighbor is ranked earlier, the correction that neighbor just received flows straight into $c$. If it is ranked later, it waits for the next iteration. For pure upwind convection numbered along the flow, $L$ contains the whole matrix. One forward sweep is the exact solve.

Cuthill–McKee is a breadth-first search#

The algorithm Cuthill and McKee published in 1969 is simple. Pick one start cell. Number cells by breadth-first search (BFS, visiting neighbors layer by layer outward). Within a layer, number cells with fewer neighbors first. Neighbor number gaps then never exceed the width of one layer.

The largest such gap is the bandwidth.

$p(i)$ is the new number of cell $i$, and the maximum runs over all cell pairs sharing a face.

In 1971 George showed that simply reversing this order reduces fill-in in direct factorization. That is the R in RCM. The start cell is usually taken from the "edge" of the graph. The standard approach is a pseudo-peripheral node search, which repeats BFS and hops to the farthest cell.

Code: 576 cells, four numberings#

We set up two problems on a 24×24 grid. One is convection-dominated, with flow at 30° and small diffusion (cell Péclet number 100). The other is pure diffusion. To mimic a mesh generator's arbitrary numbering, we shuffle the labels, then renumber with RCM. For the convection case we also try all four corners as the RCM start.

import math, random
from collections import deque
 
def build_grid(n, ux, uy, diff):
    # cell (i,j) -> neighbor list, coefficients a_nb, diagonal a_P (1st-order upwind + diffusion)
    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: lower-rank neighbors only (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: higher-rank neighbors only (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]
    # natural order uses the unshuffled grid as is
    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

Same 576 cells, same matrix, same tolerance ($10^{-8}$). The only thing that changed is the order in which cells are visited.

Upstream neighbor rank separates 25 from 9#

With flow at 30°, each cell's upstream neighbors are to the west and south. In natural order ($i$ first, then $j$), both come earlier. So a single forward sweep carries information from the southwest corner all the way to the northeast. The remaining 9 iterations clear the small diffusion and pseudo-time terms.

With random numbering, an upstream neighbor comes earlier only half the time. Information travels a few cells per sweep and then stalls. Crossing 24 cells takes many iterations. The result is 25 iterations, about 2.8×.

Bandwidth also dropped from 553 to 24. That number is a separate story from iteration count. In a real code it measures how far apart neighbor data sits in memory. Small bandwidth means better cache hits and a cheaper iteration. So RCM cuts both the iteration count and the cost per iteration. The Python test cannot measure cache effects, so here we only compare iterations.

Same bandwidth of 24, yet 9 vs. 21 iterations#

The four-corner test is more interesting. Every corner gives bandwidth 24. Yet SW and NE take 9 iterations, while SE and NW take 21.

RCM's BFS layers are diagonal bands centered on the start corner. Starting from SW, the bands sit perpendicular to the flow and rank increases along it. Starting from NE, rank simply increases against the flow. The backward sweep then runs with the flow. SGS does both sweeps, so the two cases are symmetric. For the same reason, Cuthill–McKee and reverse CM need the same number of iterations.

Starting from SE, rank grows toward the northwest. That opposes the flow's $x$ component and follows its $y$ component. Each sweep receives only one of the two upstream neighbors on time. Hence 21. Bandwidth only describes the matrix shape. It knows nothing about flow direction.

25/9 is 2.78. The title's 2.8× is that rounding. Same 576 cells, same $10^{-8}$. Bandwidth falls from 553 to 24. 553 sits next to the cell count 576, so neighbor indices span almost the whole matrix. 24 is the floor natural order already hits: on a 24×24 grid a vertical neighbor is $n=24$ away. RCM reached that floor and stopped.

All four corners have $\beta=24$. Iterations are SW 9, NE 9, SE 21, NW 21. 21/9 is 2.33, better than the random 25 and still more than double the flow-aligned 9. SE and NW share 21 because the two sweeps run opposite ways. A 9 on one diagonal means a 9 on the other.

Diffusion only differs by 13%#

Pure diffusion tells a different story. Natural order and RCM need 1549 iterations, random numbering 1747. That is a 13% gap.

Diffusion has no upstream. Information spreads equally in every direction. Gauss–Seidel is slow here because of the long wave, not the numbering. Ordering will not shrink that wave. It belongs with a Krylov solver. Ordering changes iteration counts sharply only when convection dominates. Strongly hyperbolic problems such as compressible Euler are that case.

The gap between 1549 and 1747 is 198 iterations. 198/1549 is 0.128, the 13% above. Bandwidth can stretch to 561 and the iteration count barely moves. That is the opposite of convection, where 553 produced 25. Shrinking $\beta$ and putting the upstream neighbor into $k \lt c$ are separate jobs inside one reordering. Diffusion has no upstream, so the second job is nearly empty.

So why "Reverse"?#

In the test, CM and RCM needed the same iterations. Reversing looks pointless. The payoff shows up somewhere else.

George's result concerns fill-in during direct factorization. When a banded matrix is factored with Cholesky or LU, each row fills from its first nonzero up to the diagonal. The sum of those spans, the profile, is almost always no larger for RCM than for CM. If the plan is ILU preconditioning or a direct solver, RCM is the default. If the code only runs LU-SGS, CM and RCM behave the same. The two symmetric sweeps erase the order's direction.

On a 1D line, numbering left to right is the whole job. The numbering problem starts in 2D and 3D unstructured grids. A mesh generator stamps cells in the order it creates them. That order follows point insertion or front advance. From the solver it looks close to random.

When bandwidth 24 still takes 21 iterations, drop the start cell#

Reorder once, at initialization. There is no reason to renumber during the run. Build the adjacency graph from face connectivity, then move the cell arrays, each face's left and right indices, and the boundary lists onto the new numbers. Tets, prisms, hexes: the graph is enough. The cost is one BFS, linear in the cell count.

When convection is strong, the start cell matters before the bandwidth does. In this printout every corner has $\beta=24$, yet SW and NE take 9 iterations and SE and NW take 21. A pseudo-peripheral cell is only the farthest cell on the graph. It does not know the inflow. The 9 showed up because the start happened to land on a corner along the flow axis. Nothing promises that landing. If the solution changes at machine precision after the reorder, an index map was missed. Numbering must not change the answer. The factorization that finishes once you number a 1D line from the left is in LU-SGS on a Burgers example. The numbering problem starts here.