セル番号を混ぜただけでLU-SGSが2.8倍遅くなった — RCM再順序付けと開始セルが決めるスイープ方向

LU-SGSではセル番号が、1回のスイープで情報がどこまで届くかを決めます。RCMはバンド幅を24まで縮めますが、反復9回と21回を分けるのは開始セルと流れの向きの関係です。

セル番号を混ぜると何が変わるのか

同じ格子です。セル番号だけをランダムに入れ替えました。行列は行と列の順番が変わるだけです。収束解はそのままです。それでもLU-SGS(Lower-Upper Symmetric Gauss–Seidel、下三角と上三角の2回のスイープで陰的システムを近似する手法)の反復は9回から25回になりました。差は「1回のスイープで情報がどこまで進むか」にあります。RCM(Reverse Cuthill–McKee、幅優先探索で行列のバンド幅を縮める番号付け)が直すものと直せないものを、576セルの出力で分けます。

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セルの1次風上対流拡散で色の広がり方が変わります。

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

ランダム番号では前進スイープが1回通っても色がまだらに残ります。RCM SWでは流れの矢印に沿って、1回のスイープでほぼ全体が埋まります。RCM start をSEにすると bw はそのままで、反復回数だけ伸びます。

LU-SGSの2回のスイープを式で書く#

定常線形システム $A\phi = b$ を擬似時間(pseudo-time)前進で解くとします。行列を対角と非対角に分けます。

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$ は擬似時間刻みです。

この近似分解は2回のスイープで解けます。前進スイープは番号順に進み、番号が小さい隣接セルの新しい値だけを使います。後退スイープは逆順に進み、番号が大きい隣接セルの寄与を足します。

\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$ が行列全体を含みます。前進スイープ1回が厳密解になります。

Cuthill–McKeeは幅優先探索です#

CuthillとMcKeeが1969年に発表したアルゴリズムは単純です。開始セルを1つ選びます。幅優先探索(BFS、近い隣接セルから層ごとに訪問する方法)で番号を付けます。同じ層の中では隣接数が少ないセルから付けます。こうすると、隣接セル同士の番号差は1層の幅を超えません。

この番号差の最大値がバンド幅です。

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

$p(i)$ はセル $i$ の新しい番号で、最大値は面を共有するすべてのセル対について取ります。

1971年、Georgeはこの順序を反転するだけで直接分解のフィルイン(fill-in)が減ることを示しました。これがRCMのRです。開始セルは通常、グラフの「端」から選びます。BFSを繰り返して最も遠いセルへ移っていく擬似周辺節点(pseudo-peripheral node)探索が標準です。

コード:576セル、4通りの番号付け#

24×24格子に2つの問題を設定します。1つは30°方向の流れに小さな拡散を加えた対流支配問題です(セルPéclet数100)。もう1つは純粋な拡散です。格子生成器が番号を無秩序に付けた状況をまねるため、ラベルを混ぜます。その後RCMで番号を付け直します。対流問題では、RCMの開始セルを4つの隅に変えて試します。

import math, random
from collections import deque
 
def build_grid(n, ux, uy, diff):
    # セル (i,j) -> 隣接リスト、係数 a_nb、対角 a_P (1次風上 + 拡散)
    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$ が後)では、どちらも先の順番になります。そのため前進スイープ1回で、南西の隅から北東の隅まで情報が届きます。残りの9回は、小さな拡散項と擬似時間項を片付ける分です。

ランダム番号では、上流セルが先の順番になる確率は半々です。情報は1回のスイープで数セル進んで途切れます。24セルを横切るには何度も反復が要ります。その結果が25回、約2.8倍です。

バンド幅も553から24に減りました。この数字は反復回数とは別の話です。実際のコードでは、隣接セルのデータがメモリ上でどれだけ離れているかを表します。バンド幅が小さいとキャッシュヒットが良くなり、1回の反復が速くなります。つまりRCMは反復回数と1回あたりのコストを同時に減らします。ただしPythonの実験ではキャッシュ効果を測れないので、ここでは反復回数だけを比べます。

同じバンド幅24で、反復は9回と21回#

4つの隅の実験のほうが興味深い結果です。どの隅から始めてもバンド幅は24です。それでもSWとNEは9回、SEとNWは21回です。

RCMのBFSの層は、開始した隅を中心とする斜めの帯になります。SWから始めると帯が流れに垂直に並び、順番が流れに沿って増えます。NEから始めても、順番が流れに逆らって増えるだけです。その場合は後退スイープが流れの向きになります。SGSは前進と後退の両方を行うので、2つは対称です。同じ理由で、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はその床まで縮めて、そこで止まっています。

4つの隅の $\beta$ はすべて24です。反復だけがSW 9、NE 9、SE 21、NW 21です。21/9は2.33で、ランダムの25よりは良く、流れに沿った9の2倍を超えます。SEとNWが同じ21なのは、スイープが2回、互いに逆向きだからです。一方の対角で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$ に上流を入れる仕事は、同じ再順序付けの中でも別です。拡散には上流がないので、2番目はほとんど空です。

では、なぜ「Reverse」なのか#

実験ではCMとRCMの反復回数は同じでした。反転する理由はなさそうに見えます。効く場所は別です。

Georgeが示したのは直接分解のフィルインです。帯行列をCholeskyやLUで分解すると、各行の最初の非ゼロ要素から対角までが埋まります。この「プロファイル(profile)」の合計は、RCMではほぼ常にCM以下になります。そのためILU前処理や直接法を使うつもりならRCMが標準です。LU-SGSだけを使うならCMでもRCMでも差はありません。対称な2回のスイープが、順序の向きを打ち消します。

1Dでは左から右へ番号を振れば終わりです。番号の問題は2D・3Dの非構造格子で起きます。格子生成器はセルを作られた順に番号付けします。点の挿入順やフロントの前進順に従います。ソルバーから見ればランダムな番号に近いものです。

バンド幅24でも21回で止まるなら、開始セルを捨てる#

再順序付けは初期化で1回です。計算中に番号を付け直す理由はありません。面の接続から隣接グラフを作り、セル配列と面の左右インデックス、境界リストを新しい番号へ移します。四面体でもプリズムでも、見るのはグラフだけです。コストはBFS 1回で、セル数に比例します。

対流が強いときは、開始セルがバンド幅より先です。今回の出力では4つの隅の $\beta$ がすべて24なのに、SWとNEは9回、SEとNWは21回です。擬似周辺セルはグラフ上で最も遠いセルにすぎず、流入は知りません。9回は、開始が流れの軸上の隅に当たった場合です。当たる保証はありません。再順序付けの前後で解が機械精度で変わったら、インデックス変換が漏れています。番号は解を変えてはいけません。1Dを左から番号付けすれば終わっていた分解は、Burgers例のLU-SGSにあります。その番号の問題は、ここから始まります。