◆ Zgttrf()

Sub Zgttrf	(	N As	Long,
		Dl() As	Complex,
		D() As	Complex,
		Du() As	Complex,
		Du2() As	Complex,
		IPiv() As	Long,
		Info As	Long
	)

係数行列の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による除算が発生する.

出典: LAPACK

使用例: 連立一次方程式 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

Cmplx
Function Cmplx(R As Double, Optional I As Double=0) As Complex
複素数の作成

Cimag
Function Cimag(A As Complex) As Double
複素数の虚数部

Creal
Function Creal(A As Complex) As Double
複素数の実数部

Zlangt
Function Zlangt(Norm As String, N As Long, Dl() As Complex, D() As Complex, Du() As Complex, Optional Info As Long) As Double
行列の1ノルム, フロベニウスノルム, 無限ノルム, または, 要素の最大絶対値 (複素3重対角行列)

Zgttrf
Sub Zgttrf(N As Long, Dl() As Complex, D() As Complex, Du() As Complex, Du2() As Complex, IPiv() As Long, Info As Long)
係数行列のLU分解 (複素3重対角行列)

Zgtcon
Sub Zgtcon(Norm As String, N As Long, Dl() As Complex, D() As Complex, Du() As Complex, Du2() As Complex, IPiv() As Long, ANorm As Double, RCond As Double, Info As Long)
行列の条件数 (複素3重対角行列)

Zgttrs
Sub Zgttrs(Trans As String, N As Long, Dl() As Complex, D() As Complex, Du() As Complex, Du2() As Complex, IPiv() As Long, B() As Complex, Info As Long, Optional Nrhs As Long=1)
LU分解済の連立一次方程式 AX = B, ATX = B または AHX = B の解 (複素3重対角行列)

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

RCond = 0.187722560135325

Info = 0