2.9 連立一次方程式の解法

前節までの方法によって電界に関する偏微分方程式(2-1-6)は次式の連立一次方程式になります。

A x = b (2-9-1)
ここで行列Aは複素数非対称行列であり、1行の非ゼロ要素が13個以下の疎行列です。
以下では行列の大きさをNで表します。

2.9.1 Bi-CGSTAB法

式(2-9-1)の連立一次方程式をBi-CGSTAB法(BiConjugate Gradient STABilized, 安定化双共役勾配法)を用いて解きます。 Bi-CGSTAB法のアルゴリズムは以下の通りです[6]-[10]。


(1) 適当な初期値x0を選ぶ
(2) r0 = b - A x0
(3) (r0*, r0) ≠ 0 となる初期値r0*を選ぶ
(4) p0 = r0
(5) do k = 0, 1, ...
(6)  α = (r0*, rk) / (r0*, A pk)
(7)  tk = rk - αA pk
(8)  ζ = (Ac tkc, tk) / (Ac tkc, A tk)
(9)  xk+1 = xk + αpk + ζtk
(10)  rk+1 = tk - ζAtk
(11)  ||rk+1|| < ε ならば終了する
(12)  β = (α/ζ) (r0*, rk+1) / (r0*, rk )
(13)  pk+1 = rk+1 + β(pk - ζApk)
(14) end do

図2-9-1 Bi-CGSTAB法

ここでxcは複素共役を表します。 また(x,y)はベクトルの内積(Σxiyi)を表します。
これをプログラムに即して書くと以下のようになります。


(1) 適当な初期値x0を選ぶ
(2) q = A x0
(3) r = b - q
(4) (r0*, r) ≠ 0 となる初期値r0*を選ぶ
(5) p = r
(6) c = (r0*, r)
(7) do k = 0, 1, ...
(8)  q = A p
(9)  α = c / (r0*, q)
(10)  t = r - αq
(11)  u = A t
(12)  ζ = (uc, t) / (uc, u)
(13)  x = x + αp + ζt
(14)  r = t - ζu
(15)  ||r|| < ε ならば終了する
(16)  d = (r0*, r)
(17)  β = (α/ζ) (d/c)
(18)  c = d
(19)  p = r + β(p - ζq)
(20) end do

図2-9-2 Bi-CGSTAB法(実装版)

作業ベクトルq,t,uを追加しています。 ベクトルは毎回上書きされますので添え字kを省略しています。
演算量の多い行列とベクトルの積は反復回数あたり2回((8)と(11))です。
なお、(3)以降ベクトルbは使用されませんので他のベクトル (例えばr0*)と共用することができます。
Bi-CGSTAB法では計算精度の観点から倍精度を使用する必要があります。

2.9.2 対角スケーリング

共役勾配法では対角スケーリングを行うことによって収束を速め計算時間を短縮することができます。
対角スケーリングは反復計算の前に行います。対角スケーリングによって使用メモリーが増えることはありません。


do i = 1 to N
 d = aii
 do j = 1 to N
  aij = aij / d
 end do
 bi = bi / d
end do

図2-9-3 対角スケーリング

図2-9-4に対角スケーリングがないときとあるときの収束状況の一例を示します。 対角スケーリングによって収束が速くなることがわかります。


図2-9-4 対角スケーリングと収束状況(Nx=Ny=Nz=100)

2.9.3 BLAS Level-1

図2-9-2ではベクトル同士の演算は以下の複素数倍精度のBLAS Level-1の関数で表現することができます。

表2-9-1 BLAS Level-1関数(複素数倍精度)
関数名処理
Zcopyコピー(y=x)
Zdotu内積(x,y)
Zdotc内積(xc,y)
Zaxpyy=ax+y
Dznrm22乗ノルム||x||2