◆ Zpbsv()

Sub Zpbsv	(	Uplo As	String,
		N As	Long,
		Kd As	Long,
		Ab() As	Complex,
		B() As	Complex,
		Info As	Long,
		Optional Nrhs As	Long = `1`
	)

(シンプルドライバ) 連立一次方程式 AX = B の解 (正定値エルミート帯行列)

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

ここで, AはN×N正定値エルミート帯行列, また, XおよびBはN×Nrhs行列である.

まず, コレスキー分解によりAを次のように分解する.
A = U^H * U (Uplo = "U"の場合)

A = L * L^H (Uplo = "L"の場合)

ここで, UはAと同じ上帯幅の上三角行列, LはAと同じ下帯幅の下三角行列である. 次に, 分解されたAを用いて連立一次方程式 A * X = B の解を求める.

引数

[in]	Uplo	= "U": Aの上三角部分を格納. = "L": Aの下三角部分を格納.
[in]	N	連立方程式の数, すなわち, 行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[in]	Kd	行列Aの上帯幅(Uplo = "U"の場合)あるいは下帯幅(Uplo = "L"の場合). (Kd >= 0)
[in,out]	Ab()	配列 Ab(LAb1 - 1, LAb2 - 1) (LAb1 >= Kd + 1, LAb2 >= N) [in] Kd+1×N対称帯行列形式のN×N正定値エルミート帯行列 A. Uploに従って上または下三角部分を格納. 詳細は下記を参照のこと. [out] Info = 0の場合, 帯行列Aのコレスキー分解 A = U^HU または A = LL^H の三角行列UまたはLをAと同様の形式で格納.
[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: パラメータ Uplo の誤り. (Uplo <> "U"および"L") = -2: パラメータ N の誤り. (N < 0) = -3: パラメータ Kd の誤り. (Kd < 0) = -4: パラメータ Ab() の誤り. = -5: パラメータ B() の誤り. = -7: パラメータ Nrhs の誤り. (nrhs < 0) = i > 0: Aのi×i首座小行列が正定値でないため分解を完了できなかった. また, 解は計算されなかった.
[in]	Nrhs	(省略可) 右辺の数, すなわち, 行列Bの列数. (Nrhs >= 0) (Nrhs = 0 の場合, 処理を行わずに戻る) (省略時 = 1)

詳細

次の例は, n = 6, kd = 2, Uplo = "U" の場合の対称帯行列形式を表す.

入力:

   *    *   a13  a24  a35  a46
   *   a12  a23  a34  a45  a56
  a11  a22  a33  a44  a55  a66

出力:

   *    *   u13  u24  u35  u46
   *   u12  u23  u34  u45  u56
  u11  u22  u33  u44  u55  u66

同様に uplo = "L" の場合, Aの形式は次のようになる.

入力:

  a11  a22  a33  a44  a55  a66
  a21  a32  a43  a54  a65   *
  a31  a42  a53  a64   *    *

出力:

  l11  l22  l33  l44  l55  l66
  l21  l32  l43  l54  l65   *
  l31  l42  l53  l64   *    *

*で示された配列要素は使用されない.

出典: LAPACK

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

A = ( 0.29+0.44i 0.62 -0.01-0.02i )

( 0 -0.01+0.02i 0.46 )

( 1.6236-0.7300i )

B = ( 0.1581+0.1537i )

( 0.1132-0.2290i )

とする.
Sub Ex_Zpbsv()

Const N As Long = 3, Kd = 1

Dim Ab(Kd, N - 1) As Complex, B(N - 1) As Complex

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

Ab(0, 0) = Cmplx(2.2, 0): Ab(0, 1) = Cmplx(2.11, 0): Ab(0, 2) = Cmplx(2.93, 0)

Ab(1, 0) = Cmplx(-0.32, 0.81): Ab(1, 1) = Cmplx(0.37, -0.8)

B(0) = Cmplx(-1.7808, -1.5394): B(1) = Cmplx(1.9646, 1.3687): B(2) = Cmplx(2.3245, 1.6291)

ANorm = Zlanhb("1", "L", N, Kd, Ab())

Call Zpbsv("L", N, Kd, Ab(), B(), Info)

If Info = 0 Then Call Zpbcon("L", N, Kd, Ab(), 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.59 -0.28 -0.2 -0.04 0.24 -0.49

RCond = 0.124521368143895

Info = 0