前回係数行列の仕上がり具合の比較検証の実学で終わりましたが、引続き進めるに当たって最低限の基礎知識の共有が不可欠だと判断し、座学を交えた実用的な講義に変更します。俺なりにポイントを伝授すので絶対飽きさせない自信があります。ご辛抱してお付き合いくださいね。先ず一般的なILU前処理編から始めます。Last time the practical exercise is conducted in view of the comparison and verification purpose for the incidence matrices. Now I would change to a practical lecture course as it is essential to share the basic knowledge for ensuring the meaningful continuity. I shall make the best effort to provide with a valid points on my eyes in confidence so as not to get you bored. Please be patient together with me for the practical ILU preconditioning course.

Gaussian Elimination

疎な係数行列\(A\)をもつ連立一次方程式

$$\sum_{j=1}^{n}a_{ij}x_{j} = b_{i} \hspace{20pt} i=1,\ldots,n \tag{1.0}$$

の求解を求める。

式 (1.0) で \(a_{11}\neq0\) と仮定する。第 1 式に \(m_{i}^{(1)}=a_{i1}/a_{11}\) をかけ，\(x_{1}\) の係数である \(a_{21}\) から\(a_{n1}\) までの成分を\(i≥2, j≥2\)について

$$a_{ij}^{(1)} = a_{ij} − m_{i}^{(1)}a_{1j}$$

$$b_{i}^{(1)}= b_{i} − m_{i}^{(1)}b_{1}$$

の過程で消去すると、方程式は式(1.0)と同値な

の形になる。
一般形式で表記すると、第\(k\) 段目の消去は演算は\(a_{kk}^{(k-1)}\neq0\)と仮定するとき、\(i≥k+1, j≥k+1\) に対し

$$m_{i}^{(k)} = a_{ik}^{(k-1)}/a_{kk}^{(k-1)}$$

$$a_{ij}^{(k)} = a_{ij}^{(k−1)} – m_{i}^{(k)}a_{kj}^{(k-1)}$$

$$b_{i}^{(k)} = b_{i}^{(k-1)} – m_{i}^{(k)}b_{k}^{(k-1)}$$

で計算する。

この消去を第\(1\)式から第\(n−1\)式まで行なうことで、最終的に次式が得られる。

$$ \begin{equation}
\begin{pmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
& a_{22}^{(1)} & \cdots & a_{2n}^{(1)} \\
& & \cdots & \vdots \\
0 & & & a_{nn}^{(n-1)}
\end{pmatrix}
\begin{pmatrix}
x_{1} \\ x_{2} \\ \vdots \\ x_{n}
\end{pmatrix}
=
\begin{pmatrix}
b_{1} \\ b_{2}^{(1)} \\ \vdots \\ b_{n}^{(n-1)}
\end{pmatrix}
\end{equation} $$

上式は式 (1.0) と同値である。従って、\(a_{nn}^{n-1}\neq0\)ならば、直ちに\(n\)番目の未知数\(x_{n}\)が求まり、これを直前の式に代入することで\(x_{n-1}\)が求まる。
以下、順番に\(k=n−1,\ldots,1\)に対して

$$x_{k} = {(b_{k}^{(k-1)} – \sum_{j=k+1}^{n}a_{kj}^{(k-1)}x_{j})}/a_{kk}^{(k-1)}$$

で解が求まる。
以上が ガウス消去法の手順である。

次式を用いて、順次\(i=1,\ldots,n\)を計算すると、

$$z_{i} = b_{i} − \sum_{k=1}^{i-1} l_{ik}z_{k} \tag{1.1}$$

\(z\)が求まる。求めた\(z\)を用いて、\(i=n\)から順次\(n-1,n-2,\ldots,1\)の順に次の計算行うと、

$$x_{i} = {(z_{i} – \sum_{k=i+1}^{n}u_{ik}x_{k})}/u_{ii} \tag{1.2}$$

\(x\)が求まる。前者を前進消去(Forward Elimination)、後者を後退代入(Backward Substitution)、総称してガウスの消去法(Gaussian Elimination)と呼ぶ。

また、以上の式を行列形式で解法すると、

$$Ax = b$$

に於いて、与えられた係数行列\(A\)に対して、\(A=LU\)を満たす下三角行列\(L\)と上三角行列\(U\)を見い出すことを、\(A\)を\(LU\)分解するという。\(A\)が\(LU\)分解されている時、作業用のベクトル\(z\)を用意し、

$$Lz = b$$

を\(z\)について解いた後、

$$Ux = z$$

を\(x\)について解くと、最終的な解が得られる。

依って、ガウスの消去法は下記アルゴリズム1.0に総括される。

ALGORITHM 1.0 Gaussian Elimination
1.   For i=2,...,n Do:
2.      For k=1,...,i−1 Do:
3.         aik := aik / akk
4.         For j=k+1,...,n Do:
5.            aij := aij − aik ∗ akj
6.         EndDo
7.      EndDo
8.      bi := bi + aik * bk
9.   EndDo

この前進消去過程では、右辺\(b\)の値も\(A\)と一緒に変形する。実際の数値計算では、同じ\(A\)に対して何度も異なる\(b\)を使って連立 1 次方程式を解く事になるので、このアルゴリズムをそのままCodingすると、\(b\)に応じて幾度も同じ\(A\)を変形する事になって効率が悪い。そこで、予め行列\(A\)を\(LU\)分解し、得られた三角行列\(L\)と\(U\)を用いて\(b\)を変形するコードに書き直す。こうすると、異なる\(b\)が与えられた場合でも、先に求めた\(LU\)を利用して効率的に求解できる。

Ordering

\({(i,j,k)}\)軸方向に応じて順々に順序付けしていく事(Ordering)を自然な順序付け(Natural Ordering)と云う。その他の有効な順序付けとしてRed-Black順序付け(Red-Black Ordering)がある｡過去には、D4順序付けとよばれていた方法であるが､解くべき変数を約半分に削減できるので大幅な記憶領域、演算量の削減を図る事ができる。

Red-Black順序付けは次の手順にて行う。
先ず､格子全てにBlackラベルをつけ､１番目の格子にRedラベルをつける｡次に、自然な順序付けに従い２番目の格子を調べる｡もし､この格子がBlackラベル格子にのみ囲まれていればこの格子にRedレベルを付ける。もし､そうでなければBlackラベルのままとする｡このプロセスを格子全てにつき繰り返す。即ち、格子点の座標を\(i,j\)とするとき、先ずは\(i+j\)が奇数となる格子点を番号付けし、その後で\(i+j\)が偶数となる格子点を番号付けする事で得られる。

ここで､例題４-３-２にRed-Black順序付けを適用してみる｡
Red-Black順序付けを適用した格子システムを図2.0に示す。尚、本例では坑井はないと仮定している。

係数行列の構築方法を理解する為､格子B2について考察してみよう｡
格子B2は上流格子R3に繋がっているが､Red-Black順序付けを適用した場合、得られる縮小係数行列ではRed格子を除きBlack格子のみで構成されるので、結果的にR3の上流格子に該当するB3とB6がR3の代わりにB2とコネクションを形成する事になる。

得られるRed-Black係数行列を図2.1-1、Black-Black縮小係数行列を図2.1-2に示す。

次に､Red-Black順序付けの理論的意義について下記に述べる｡Red-Black順序付けを適用した係数行列は次の形式で示される。

$$ \begin{equation}
\begin{pmatrix}
J_{RR} & J_{RB} \\
J_{BR} & J_{BB}
\end{pmatrix}
\begin{pmatrix}
X_{R} \\ X_{B}
\end{pmatrix}
=
\begin{pmatrix}
R_{R} \\ R_{B}
\end{pmatrix}
\end{equation} \tag{2.0}$$

ここで、
\(J_{RR}\): 係数行列におけるRed-Red格子
\(J_{BR}\): 係数行列におけるBlack-Red格子
\(J_{RB}\): 係数行列におけるRed-Black格子
\(J_{BB}\): 係数行列におけるBlack-Black格子
\(X_{R}\): Red格子における解
\(X_{B}\): Black格子における解
\(R_{R}\): Red格子における残差
\(R_{B}\): Black格子における残差

対角行列で\(J_{BR}\)を消去すると､先の係数行列は次のような縮小行列となる。

$$ \begin{equation}
\begin{pmatrix}
J_{RR} & J_{RB} \\
0 & J^{’}_{BB}
\end{pmatrix}
\begin{pmatrix}
X_{R} \\ X_{B}
\end{pmatrix}
=
\begin{pmatrix}
R_{R} \\ R^{’}_{B}
\end{pmatrix}
\end{equation} \tag{2.1}$$

ここで、

$$\begin{align}J_{BB}^{’} &= J_{BB} – J_{BR} J_{RR}^{-1} J_{RB} \\\tag{2.2}\\
R_{B}^{’} &= R_{B} – J_{BR} J_{RR}^{-1} R_{R} \end{align}$$

Black格子解は次式から求められる。

$$X_{B} = J^{’}_{BR}R_{B} \tag{2.3}$$

又､Red格子解は次式から逐次代入法によって求められる。

$$X_{R} = J^{-1}_{RR}(R_{R} – J_{RB}X_{R}) \tag{2.4}$$

Symbolic Factorization

ガウス消去法の\(LU\)分解の過程に於いて元々の係数行列で要素の値が零である部分に零以外の値が出現することになる。これをフィルイン(Fill-in)という。\(LU\)分解によって掃き出し成分が生じる位置をSymbolicallyに決定するプロセスをシンボリック分解(Symbolic Factorization)といい、ワッツ(1981年Watts)により提唱された｡Symbolic Factorizationは係数行列又は非零成分変更の都度必要となる｡具体的には、坑井が開閉した場合とか､\(ILU\)分解のレベルを変えた場合などに再度算定する必要がある。

係数行列を\(LU\)分解する際､行列其々の非零成分の構造を調べ､その位置を記憶する配列を確保しておくが、このフィルインの位置(レベル)によって収束性如いては演算速度が左右される。フィルインをなるべく少なくするように必要な記憶領域と演算量を削減することが普通である。

係数行列の元の非零成分\(a_{ij}\)のレベルを０レベルとすると、係数行列\(a_{ij}\)の全成分のフィルインのレベル\(lev_{ij}\)は次式で示される。

$$lev_{ij} =
\begin{cases}
0 & if ~a_{ij}\neq0,~or~i=j \\
∞ & otherwise
\end{cases} \tag{3.0}$$

係数行列\(a_{ij}\)に対してアルゴリズムに従ってガウス消去法を実行すると､掃き出し成分のレベルはそれを掃き出す成分らのレベルの総和であるから､第\(k\)段目の部分軸における掃き出し成分\(a_{ij}\)のレベルは一般に次式で与えられ、

$$lev_{ij}=min\lbrace{lev_{ij}, lev_{sum}=lev_{ik}+lev_{kj}+1}\rbrace \tag{3.1}$$

フィルインレベル(p)のILUアルゴリズムは3.0に更新される。

ALGORITHM 3.0 ILU Factorization(p)
1.   Define levij = 0 for all nonzero elements aij
2.   For i=2,...,n Do:
3.      For each k=1,...,i−1 and for levik ≤ p Do:
4.         aik := aik / akk
5.         ai∗ := ai∗ − aik * ak∗
6.         Update levels of fill of the nonzero aij's using (1.4)
7.      EndDo
8.      Replace any element in row i with levij > p by zero
9.   EndDo

ここで､例題４-３-２係数行列に対して完全にシンボリック分解を行ってみよう。

例えば、図3.0に示す左下三角行列に於ける第４行目の掃き出し成分のレベルは､式(3.1)より図3.1のように算定される。図中、xはフィルインの可能性を示す。

i,j	Non Zero Location	levij
(4,1)	(1,2) (1,4) (1,8)	min{ levij=(0) (1) (0), levsum=(2) (1) (2) }
(4,3)	(3,4) (3,10)	min{ levij=(1) (0), levsum=(1) (2) }
(4,2)	(2,5) (2,8) (2,9) (2,4)	min{ levij=(1) (2) (0) (1), levsum=(1) (2) (3) (1) }

以降同様に算定し、ガウス消去法が終了した段階で全ての掃き出し成分のレベルが求まる。最終的な結果を図3.2に示す。ここで、1wは坑井のよるフィルインを示す。

以上、ILU前処理は、オーダリング、シンボリック分解、\(LU\)分解の三段階の処理を順次行うことにより実行される。

以上
つづく