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) 適当な初期値x₀を選ぶ
(2) r₀ = b - A x₀
(3) (r₀^*, r₀) ≠ 0 となる初期値r₀^*を選ぶ
(4) p₀ = r₀
(5) do k = 0, 1, ...
(6) 　α = (r₀^*, r_k) / (r₀^*, A p_k)
(7) 　t_k = r_k - αA p_k
(8) 　ζ = (A^c t_k^c, t_k) / (A^c t_k^c, A t_k)
(9) 　x_k+1 = x_k + αp_k + ζt_k
(10) 　r_k+1 = t_k - ζAt_k
(11) 　||r_k+1|| < ε ならば終了する
(12) 　β = (α/ζ) (r₀^*, r_k+1) / (r₀^*, r_k)
(13) 　p_k+1 = r_k+1 + β(p_k - ζAp_k)
(14) end do

図2-9-1 Bi-CGSTAB法

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

(1) 適当な初期値x₀を選ぶ
(2) q = A x₀
(3) r = b - q
(4) (r₀^*, r) ≠ 0 となる初期値r₀^*を選ぶ
(5) p = r
(6) c = (r₀^*, r)
(7) do k = 0, 1, ...
(8) 　q = A p
(9) 　α = c / (r₀^*, q)
(10) 　t = r - αq
(11) 　u = A t
(12) 　ζ = (u^c, t) / (u^c, u)
(13) 　x = x + αp + ζt
(14) 　r = t - ζu
(15) 　||r|| < ε ならば終了する
(16) 　d = (r₀^*, 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は使用されませんので他のベクトル（例えばr₀^*）と共用することができます。
Bi-CGSTAB法では計算精度の観点から倍精度を使用する必要があります。

2.9.2 対角スケーリング

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

do i = 1 to N
　d = a_ii
　do j = 1 to N
　　a_ij = a_ij / d
　end do
　b_i = b_i / 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	内積(x^c,y)
Zaxpy	y=ax+y
Dznrm2	2乗ノルム\|\|x\|\|₂