◆ Zgesv()

Sub Zgesv	(	N As	Long,
		A() As	Complex,
		IPiv() As	Long,
		B() As	Complex,
		Info As	Long,
		Optional Nrhs As	Long = `1`
	)

(シンプルドライバ) 連立一次方程式 AX = B の解 (複素行列)

目的: 本ルーチンは次の複素連立一次方程式を解く.
A * X = B

ここで, AはN×N行列, また, XおよびBはN×Nrhs行列である.

まず, 行交換によるピボットの部分選択を行うLU分解を用いて, 次のようにAを分解する.
A = P * L * U

ここで, Pは置換行列, Lは対角要素が1の下三角行列, そして, Uは上三角行列である. 次に, 分解されたAを用いて連立方程式 A * X = B の解を求める.

引数

[in]	N	連立方程式の数, すなわち, 行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in,out]	A()	配列 A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= N) [in] N×N係数行列 A. [out] 分解 A = PLU のLおよびU. Lの対角要素(= 1)は格納されない.
[out]	IPiv()	配列 IPiv(LIPiv - 1) (LIPiv >= N) 置換行列Pを定義するピボットインデックス. 第i行が第IPiv(i-1)行と交換されたことを表す.
[in,out]	B()	配列 B(LB1 - 1, LB2 - 1) (LB1 >= max(1, N), LB2 >= Nrhs) (2次元配列) または B(LB - 1) (LB >= max(1, N), Nrhs = 1) (1次元配列) [in] N×Nrhs右辺行列 B. [out] Info = 0 の場合, N×Nrhs解行列 X.
[out]	Info	= 0: 正常終了. = -1: パラメータ N の誤り. (N < 0) = -2: パラメータ A() の誤り. = -3: パラメータ IPiv() の誤り. = -4: パラメータ B() の誤り. = -6: パラメータ Nrhs の誤り. (Nrhs <= 0, または, Nrhs <> 1 かつ B()が1次元配列) = i > 0: Uのi番目の対角要素が0である. 分解を完了したが, Uが特異であるため解を計算できなかった.
[in]	Nrhs	(省略可) 右辺の数, すなわち, 行列Bの列数. (Nrhs >= 0) (Nrhs = 0 の場合, 処理を行わずに戻る) (省略時 = 1)

出典: LAPACK

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

A = ( -0.8-0.92i -0.29+0.86i 0.64+0.51i )

( 0.71+0.59i -0.15+0.19i 0.2+0.94i )

( -0.5853-0.9457i )

B = ( -2.1697-1.0006i )

( 0.0116-0.5094i )

とする.
Sub Ex_Zgesv()

Const N = 3

Dim A(N - 1, N - 1) As Complex, B(N - 1) As Complex, IPiv(N - 1) As Long

Dim ANorm As Double, RCond As Double, Info As Long

A(0, 0) = Cmplx(0.2, -0.11): A(0, 1) = Cmplx(-0.93, -0.32): A(0, 2) = Cmplx(0.81, 0.37)

A(1, 0) = Cmplx(-0.8, -0.92): A(1, 1) = Cmplx(-0.29, 0.86): A(1, 2) = Cmplx(0.64, 0.51)

A(2, 0) = Cmplx(0.71, 0.59): A(2, 1) = Cmplx(-0.15, 0.19): A(2, 2) = Cmplx(0.2, 0.94)

B(0) = Cmplx(-0.5853, -0.9457): B(1) = Cmplx(-2.1697, -1.0006): B(2) = Cmplx(0.0116, -0.5094)

ANorm = Zlange("1", N, N, A())

Call Zgesv(N, A(), IPiv(), B(), Info)

If Info = 0 Then Call Zgecon("1", N, A(), ANorm, RCond, Info)

Debug.Print "X =", 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
複素数の実数部

Zlange
Function Zlange(Norm As String, M As Long, N As Long, A() As Complex, Optional Info As Long) As Double
行列の1ノルム, フロベニウスノルム, 無限ノルム, または, 要素の最大絶対値 (複素行列)

Zgecon
Sub Zgecon(Norm As String, N As Long, A() As Complex, ANorm As Double, RCond As Double, Info As Long)
行列の条件数 (複素行列)

Zgesv
Sub Zgesv(N As Long, A() As Complex, IPiv() As Long, B() As Complex, Info As Long, Optional Nrhs As Long=1)
(シンプルドライバ) 連立一次方程式 AX = B の解 (複素行列)

実行結果: X = 0.79 -0.13 0.13 0.75 -0.91 0.3

RCond = 0.250214147937285

Info = 0