◆ Dgttrf()

係数行列の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.58 1.14 0 ) ( -0.5960 )

A = ( -1.56 2.21 0.16 ), B = ( -0.6798 )

( 0 0.24 0.25 ) ( -0.0406 )

とする.
Sub Ex_Dgttrf()

Const N = 3

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

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

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

Dl(0) = -1.56: Dl(1) = 0.24

D(0) = -0.58: D(1) = 2.21: D(2) = 0.25

Du(0) = 1.14: Du(1) = 0.16

B(0) = -0.596: B(1) = -0.6798: B(2) = -0.0406

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

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

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

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

Debug.Print "X =", B(0), B(1), B(2)

Debug.Print "RCond =", RCond

Debug.Print "Info =", Info

End Sub