« 経済成長と減税を可能にする魔法の税制---累進所得税(2) | トップページ | 大気イオン地表面濃度の異常と地電流・地電位の関係---アイデア(1) »

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

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

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

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

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

計算時間ですが、以下に述べるようなn次正方行列Aとベクトルbに対して、うちの非力なPC(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万要素くらいの定常問題なら、こうしたPCでもなんとかなるでしょうか。

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

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


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

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

Time
図1 (クリックで拡大)

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

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


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

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

Ilu_sor_cp
図2 (クリックで拡大)

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

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

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

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


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

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

Store_new
図3 (クリックで拡大)

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


■モジュール(BiCGStabSOR.f90)の使用法

1次方程式 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) 第1行目、第2行目、...の順に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に格納することもできます。


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

末尾のサンプルプログラム 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
<<<<<


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

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

Uzu
図4 (クリックで拡大)

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

V10
図5 (クリックで拡大)

V90
図6 (クリックで拡大)

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

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

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


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

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

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

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

|

« 経済成長と減税を可能にする魔法の税制---累進所得税(2) | トップページ | 大気イオン地表面濃度の異常と地電流・地電位の関係---アイデア(1) »

「雑感」カテゴリの記事

コメント

「BiCGStab法の記事を拝見して」
 大変勉強になりました。当方では有限要素法プログラムを開発販売していますが、やはり大規模構造におけるiLUの計算時間が非常に長いと感じています。本記事の計算手法を当方でも試してみました結果、残念ながら失敗に終りました。上記計算手法がうまく行く場合とうまく行かない場合があるように思われます。

 当方で試してみました解析種類は、3次元流れの定常解析です(非圧縮で非ニュートン流)。ナビア・ストークス方程式を解くもので、テストモデルはFEM行列の次元数2815、フルバンド幅1099、保存インデックス数437153です。この種のFEMは反復法での収束性があまり良くありません。当方のiLU付きBiCGStab(2)法ソルバーでは解が得られています。本記事の計算手法をテストして調べましたら、内部反復SORソルバーの繰り返し計算で解が全く良化していませんでした。これがうまく行かない原因ではないかと考えています。また、計算時間は、内部反復の計算時間との兼ね合いがありますので、一概にiLUより計算時間が短くなるとは言えないようです。
 注)構造強度解析では、反復法で収束しない境界条件があることはわかっています。

 以上、うまく行かなかった報告で申し訳ありませんが、今後も本記事のような新手法に挑戦され、結果を公開していただければ幸いです。

投稿: 黒田 英夫 | 2009.10.14 22:26

黒田 英夫さま

たいへん参考になるコメントをありがとうございます。

NS方程式の定常問題で、iLU付きBiCGStab(2)では収束するが、SOR-BiCGStabでは収束しないケースがあるとのこと、勉強になりました。 乱流に近いケースでは、ゆらぎが遠くまで伝わるので、非対角項が重要になってくる。 こうした場合には、対角優位を前提にしたSOR内部反復はうまくいかない、ということでしょうか。

WSはいま、大気イオンの接地境界層付近での乱流拡散を考えておりますが、やはり定常問題の収束性の問題にちょっと悩んでいます。境界条件の問題かと思っていましたが、そうではないのかも。 イオンどうしのクーロン相互作用が長距離に及ぶために、係数行列の非対角項が効いてくることが原因なのかも知れません。

一見、時間がかかるように思えても、ILU分解のような前処理をしたほうがよいケースがあるのですね。

若葉マークなので路肩に突っ込んだりするかも知れませんが(笑)、今後とも、アドバイスなどいただけましたら幸いです。

投稿: Wave of sound | 2009.10.16 23:42

WS様
 貴方のご返答を拝見して、直接の回答にはなりませんが、当方コメントを追加させていただきます。

 内部反復SORがうまく行くには、そもそもSORソルバー単独で解を良化できることが前提になるのではないかと考えています。すなわち、SORソルバーで解を良化できないFEM行列には、内部反復SORを適用できないと思います。

 実は、本日、別のFEM対称行列(構造強度解析用)ですが、対角スケーリング前処理のCG法を試してみました。この場合も、CG法反復の解はほとんど良化しませんでした。もちろん、ICCG法なら問題なく収束解が得られます。

 従って、iLU付きのCG法(BiCGStab法などを含む)なら高い確率で収束解が得られますが、内部反復SORのように前処理性能が劣る場合は、その前処理自体で解が良化するかどうか注意する必要があります。その前処理で解が良化するFEM行列であれば、十分実用できると考えます。ただし、応用範囲に制限が付きます。

 以上は私見ですので、異なるご見解などがございましたらお聞かせください。

投稿: 黒田 英夫 | 2009.10.21 17:23

黒田 英夫さま

補足をありがとうございます。 前処理自体で解が良化するかどうか注意する必要があるとのこと、今後の参考にさせていただきます。

いまWSの頭は、モデル化のレベル(乱流モデルの選択やイオン間の相互作用)にあり、解きたい問題のプログラムを走らせるのは、もう少し先になりますが、そのときにいただいたアドバイスが役にたちそうです。

投稿: Wave of sound | 2009.10.26 01:19

コメントを書く



(ウェブ上には掲載しません)


コメントは記事投稿者が公開するまで表示されません。



トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/52086/45864263

この記事へのトラックバック一覧です: BiCGSTAB法による疎行列1次方程式のソルバ:

« 経済成長と減税を可能にする魔法の税制---累進所得税(2) | トップページ | 大気イオン地表面濃度の異常と地電流・地電位の関係---アイデア(1) »