BiCGSTAB法による疎行列1次方程式のソルバ

大気電場の成因（雷雲内での電荷分離のメカニズム？）とか、開口端補正に対する粘性の影響など、いま考えている問題はいくつかあるのですが、手計算によるWSの解析的なアプローチはこのところ進展していません。

そんな折、etopirica5さんが円筒管の開口端補正についてのレビンとシュヴィンガーの理論値（半径×0.61）を、放射境界条件を用いた数値計算（有限要素法による音場解析）で再現しておられるのを発見。 アウトプットの質、量、スピードともにすごい方ですw

触発されて、上のような（WSにとっての）未解決問題に、数値計算でアプローチしてみよう、と思い立ちました。

まずは基本的な道具づくりからと、非対称疎行列を係数に持つ１次方程式Ax=bのソルバを作ってみました。 BiCGStab法で、前処理内部反復にSORを用いています。 記事の末尾に、ソルバ（f90モジュール）とサンプルプログラムを置いておきます。 ご自由にお使い下さい。 ただし、もし、バグなどあっても責任はとれませんので、ご利用は自己責任でお願いします。 今回は、開発のてんまつとモジュールの使い方を書きます。

計算時間ですが、以下に述べるようなn次正方行列Aとベクトルbに対して、うちの非力なＰＣ（1.5GHz PowerPC G4）で、行列要素の値のセットと高速化のためのインデックス作成を含めて
n=10000のとき t=10.2秒
n=100000のとき t=152.5秒=2.54分
n=1000000のとき t=2300秒=38.35分
となりました。
計算時間はおよそnの1.18乗に比例しています。 2次元の100万要素くらいの定常問題なら、こうしたＰＣでもなんとかなるでしょうか。

（試した行列Aは次のような行列です。A ≡ D + F、ただし、Dは対角行列で、各対角要素は2以上3以下の一様乱数で与える。 Fは各行に-1以上1以下の一様乱数をランダムに10個配置した疎行列。 また、ベクトルbの各要素は2以上3以下の一様乱数で与える。）

開発段階では、行列の格納方法を変えるだけで桁違いの高速化が達成できたりしました。 まだまだ、改良の余地があると思います。 アドバイスをいただけたらうれしいです。

■最初の試み---ILU（不完全LU）分解による前処理つきのBiCGSTAB法

BiCGSTAB法というのは、係数行列が正値対称な場合のCG法（共役勾配法）を非対称行列に拡張したものです。 収束をよくするための前処理として、係数行列Aの不完全LU分解を併用する方法を最初に試しました。

図１　（クリックで拡大）

次数nが大きくなるにつれて、ILU分解に時間がかかり、計算時間はおよそnの２乗に比例して増えます。 Aが疎行列であることを思えば、計算時間はnの１乗に比例してほしいところです。 ILU分解のアルゴリズムの再検討も行ったのですが、高速化する方法がわからなかった(^^;;

そこで、ILU分解を使わない方法を検討して、内部反復にSORを用いる方法を試しました。 対角優位の係数行列の場合にはとくに良さそうです。

■内部反復にSORを用いるBiCGSTAB+SORのハイブリッド

阿部邦美氏・張紹良氏・杉原正顯氏の発表：前処理のための内部反復をもつ一般化共役残差法とその内部解法に関する一考察 を参考に、BiCGSTABとSORを組み合わせた方法を試してみました。

図２　（クリックで拡大）

この方法だと事前にILU分解を行う必要がないので、計算時間が約４割短くなりました。

しかし、依然として、計算時間はｎの１乗ではなく、ｎの２乗に比例して増えます。どうしてなんだろう、と思って調べてみると、１次方程式を解く反復部分ではなくて、行列Aの非零要素を配列にセットし、行圧縮形式(Compressed Sparse Row (CSR) format)(*)に並べ替え、高速化のために対角要素のインデックスを作成する部分で時間を要していることがわかりました。

速く走ることばかり考えていたら、実はその前のシューズをはく部分で時間がかかっていた、みたいな感じです。

*) CSR形式については、以下を参照。
Yousef Saad, Iterative Methods for Sparse Linear Systems (iter1.pdf, chap3.4)

■格納順序を工夫して高速化

もともとは、行列Aの要素をどんな順序でセットしてもよいプログラムだったのですが、上で見たように、それではあまりにも時間を浪費するので、少し使いやすさを犠牲にすることにしました。 行列Aの要素は、まず第１行目の要素たち、つぎに第２行目の要素たち、... のように、上から順番にセットしなければならない、という制限をつけて、そのかわり高速化を可能にしました。 （同一行のなかでは、要素はどんな順にセットしても構いません。 上から順にセットする必要はあるが、左右の順序は任意です。）

図３　（クリックで拡大）

これで計算時間は大幅に短縮し、次数ｎの約1.18乗に比例するようになりました。

■モジュール（BiCGStabSOR.f90）の使用法

１次方程式 Ax=b について

行列Aの次数n=10000
Aの非零要素を格納する配列の大きさ（概算で大きめに宣言）nf0=300000
内部でのSOR反復の回数itn=10
収束条件（残差÷初期残差）eps=1.0d-30
SORの加速係数omega=1.5

の場合で説明します。

----------
1) use BiCGStabSOR ... モジュールの使用を宣言します
2) real(8) :: b(n), x(n) ... 右辺のベクトルb, 解ベクトルxを格納する配列を宣言します
2) call allocate_aa(n, nf0) ... 行列Aの要素を格納する配列を準備します。最初は全ての要素がゼロです
3) 第１行目、第２行目、...の順にAの要素をセットします。 対角要素は非零である必要があります。
　　do i = 1, n
　　　call add_aa(i, j1, val1) ... Aの第(i,j1)要素に値val1を加えます
　　　call add_aa(i, j2, val2) ... Aの第(i,j2)要素に値val2を加えます
　　　...
　　　call make_index_row(i) ... Aの第i行目のインデックスを作成します
　　end do
4) 右辺のベクトルbをセットします。
　　b(1) = u1, b(2) = u2, ....
5) 解ベクトルxの初期値をセットします。（ゼロベクトルでもよい）
　　x(1) = v1, x(2) = v2, ....
6) モジュールのメインルーチンを呼び出します。
　　call BiCGStab(b, x, eps, itn, omega) ... xに解がセットされます
7) ここで解xの表示などを行います。
8) call deallocate_aa() ... 行列Aを格納していた配列を解放します
----------

解は、6)のメインルーチンBiCGStabで求まります。 この解xに、行列Aを左から乗じたベクトルy=Axを計算したければ、以下のようにサブルーチンtrans_Aを呼び出してください。（事前にreal(8) :: y(n)を宣言しておく必要があります）

call trans_A(x, y)

すると積Axがyにセットされます。

同様に、「call trans_A(z, y)」によって、Aと任意のベクトルzの積を計算してyに格納することもできます。

■使用例１　サンプルプログラム

末尾のサンプルプログラム test_BiCG.f90 とモジュール BiCGStabSOR.f90 を同一ディレクトリに置き、そのディレクトリに移動します。 コンパイルして実行します。 以下はg95を用いた場合です。

>>>>>
$ g95 -c BiCGStabSOR.f90
$ g95 test_BiCG.f90 BiCGStabSOR.o
$ time ./a.out
行列Aをセットし、インデックスを作成
連立方程式を解く
BiCGStab : iteration = 1 err = 0.3439110121467133
BiCGStab : iteration = 2 err = 0.1164432737586217
BiCGStab : iteration = 3 err = 0.06451884273381157
BiCGStab : iteration = 4 err = 0.020630954370641726
BiCGStab : iteration = 5 err = 0.006727731895106646
...
BiCGStab : iteration = 62 err = 1.8051163991796378E-28
BiCGStab : iteration = 63 err = 1.1466359963051961E-28
BiCGStab : iteration = 64 err = 1.5807610755819108E-29
BiCGStab : iteration = 65 err = 7.465540104283266E-30
BiCGStab : iteration = 66 err = 4.876834841150964E-31
i x_i (ax)_i b_i (ax-b)_i
1 -0.22718379891456397 2.6356894031632683 2.635689403163269 -8.881784197001252E-16
2 0.7909582558745546 2.0326006598770623 2.032600659877062 4.440892098500626E-16
3 0.22327875208768816 2.7630785314831896 2.7630785314831883 1.3322676295501878E-15
4 -0.4923048950956566 2.9563035743776704 2.956303574377671 -4.440892098500626E-16
5 -0.5758304521678598 2.7402388919144904 2.740238891914487 3.552713678800501E-15
.....
9996 -0.6908738292005479 2.165197746362537 2.165197746362537 0.
9997 3.3107586353027236 2.87137534352951 2.8713753435295075 2.6645352591003757E-15
9998 1.1468878038275612 2.7981291597243425 2.79812915972434 2.6645352591003757E-15
9999 1.1213607440230315 2.2077956653665733 2.207795665366575 -1.7763568394002505E-15
10000 1.2789141395192902 2.5819734537508343 2.581973453750834 4.440892098500626E-16
real 0m12.646s
user 0m10.202s
sys 0m0.317s
<<<<<

■使用例２　移流拡散方程式

２次元対流の中で粒子が拡散する場合の定常濃度分布を計算します。 正方形領域内に図のような対流があるとします。

図４　（クリックで拡大）

拡散スピードにくらべて対流の速度が小さい場合（図５）と、大きい場合（図６）の定常濃度分布を、ｚ軸に濃度をとって示しました。

図５　（クリックで拡大）

図６　（クリックで拡大）

境界条件は、奥の辺（高濃度）と手前の辺（低濃度）ではディリクレ条件、左右の辺ではノイマン条件としています。

拡散スピードにくらべて対流の速度が大きいと、ディリクレ条件を設定した辺の近傍で濃度勾配がとても大きくなります。 メッシュが粗いと濃度変化に追いつきません。 そうした場合にはメッシュを細かくして計算するか、あるいは、Scharfetter-Gummelの離散化スキーム(*)を使うとよいようです。

*) Scharfetter-Gummel Discretization Scheme for Drift-Diffusion Equations

■モジュールとテストプログラム

最後に、上で紹介したモジュール BiCGStabSOR.f90 とテストプログラム test_BiCG.f90 を置いておきます。

ダウンロード BiCGStabSOR.f90 (4.4K)

ダウンロード test_BiCG.f90 (1.7K)