◆ Zgttrf()

係数行列のLU分解 (複素3重対角行列)

目的: 本ルーチンはピボットの部分選択および行交換による消去法を用いて複素3重対角行列AのLU分解を計算する. 分解は次の形式である.
A = L * U

ここで, Lは置換行列と対角成分が1の下2重対角行列の積, そして, Uは対角要素と第1・第2上副対角要素のみが0でない上三角行列である.

引数

[in]	N	行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in,out]	Dl()	配列 Dl(LDl - 1) (LDl >= N - 1) [in] AのN-1個の下副対角要素. [out] AのLU分解結果の行列Lを定義するN-1個の乗数.
[in,out]	D()	配列 D(LD - 1) (LD >= N) [in] A の対角要素. [out] AのLU分解結果の上三角行列UのN個の対角要素.
[in,out]	Du()	配列 Du(LDu - 1) (LDu >= N - 1) [in] A のN-1個の上副対角要素. [out] U のN-1個の上副対角要素.
[out]	Du2()	配列 Du2(LDu2 - 1) (LDu2 >= N - 2) UのN-2個の第2上副対角要素.
[out]	IPiv()	配列 IPiv(LIPiv - 1) (LIPiv >= N) ピボットインデックス. 1 <= i <= N について, 行列の第i行は第IPiv(i-1)行と交換されたことを表す. IPiv(i-1)は常にiまたはi+1である. IPiv(i-1) = i は行の交換が不要であったことを示す.
[out]	Info	= 0: 正常終了. = -1: パラメータ N の誤り. (N < 0) = -2: パラメータ Dl() の誤り. = -3: パラメータ D() の誤り. = -4: パラメータ Du() の誤り. = -5: パラメータ Du2() の誤り. = -6: パラメータ IPiv() の誤り. = i > 0: Uのi番目の対角要素が0である. 分解を完了したがUが特異であり, 連立方程式の解の計算に使用すると0による除算が発生する.

使用例: 連立一次方程式 Ax = B を解き, 同時にAの条件数の逆数の推定値(RCond)を求める. ただし,
( 0.81+0.37i 0.20-0.11i 0 )

A = ( 0.64+0.51i -0.80-0.92i -0.93-0.32i )

( 0 0.71+0.59i -0.29+0.86i )

( -0.0484+0.2644i )

B = ( -0.2644-1.0228i )

( -0.5299+1.5025i )

とする.
Sub Ex_Zgttrf()

Const N = 3

Dim Dl(N - 2) As Complex, D(N - 1) As Complex, Du(N - 2) As Complex

Dim Du2(N - 3) As Complex, IPiv(N - 1) As Long

Dim B(N - 1) As Complex, ANorm As Double, RCond As Double, Info As Long

Dl(0) = Cmplx(0.64, 0.51): Dl(1) = Cmplx(0.71, 0.59)

D(0) = Cmplx(0.81, 0.37): D(1) = Cmplx(-0.8, -0.92): D(2) = Cmplx(-0.29, 0.86)

Du(0) = Cmplx(0.2, -0.11): Du(1) = Cmplx(-0.93, -0.32)

B(0) = Cmplx(-0.0484, 0.2644): B(1) = Cmplx(-0.2644, -1.0228): B(2) = Cmplx(-0.5299, 1.5025)

ANorm = Zlangt("1", N, Dl(), D(), Du())

Call Zgttrf(N, Dl(), D(), Du(), Du2(), IPiv(), Info)

If Info = 0 Then Call Zgttrs("N", N, Dl(), D(), Du(), Du2(), IPiv(), B(), Info)

If Info = 0 Then Call Zgtcon("1", N, Dl(), D(), Du(), Du2(), IPiv(), ANorm, RCond, Info)

Debug.Print "X =",

Debug.Print Creal(B(0)), Cimag(B(0)), Creal(B(1)), Cimag(B(1)), Creal(B(2)), Cimag(B(2))

Debug.Print "RCond =", RCond

Debug.Print "Info =", Info

End Sub

実行結果: X = -0.15 0.19 0.2 0.94 0.79 -0.13

RCond = 0.187722560135325

Info = 0