Zcgesv Zgbsv Zgbsvx Zgesv Zgesvx Zgtsv Zgtsvx Zspsv Zspsvx Zsysv Zsysvx

Zgbsvx() Sub Zgbsvx ( Fact As  String, Trans As  String, N As  Long, Kl As  Long, Ku As  Long, Ab() As  Complex, Afb() As  Complex, IPiv() As  Long, Equed As  String, R() As  Double, C() As  Double, B() As  Complex, X() As  Complex, RCond As  Double, FErr() As  Double, BErr() As  Double, Info As  Long, Optional Nrhs As  Long = 1, Optional RPiv As  Double  ) (エキスパートドライバ) 連立一次方程式 AX = B の解 (複素帯行列) 目的 本ルーチンはLU分解を用いて次の複素連立一次方程式を解く. A * X = B または A^T * X = B ここで, Aは下帯幅kl, 上帯幅kuのn×n帯行列, また, XおよびBはn×nrhs行列である. 解の誤差限界および推定条件数も求められる. 解説 本ルーチンでは以下の手順で計算が行われる. Fact = "E"の場合, 連立方程式を均衡化するためにスケーリング係数を計算する. Trans = "N": diag(R)*A*diag(C)*inv(diag(C))*X = diag(R)*B Trans = "T"または"C": (diag(R)*A*diag(C))^T *inv(diag(R))*X = diag(C)*B 連立方程式を均衡化するかどうかは行列Aのスケーリングに依存する. しかし, 均衡化する場合, Aは diag(R)*A*diag(C) により, Bは diag(R)*B (Trans = "N"の場合) または diag(C)*B (Trans = "T"または"C"の場合) により上書きされる. Fact = "N"または"E"の場合, (Fact = "E"の場合には均衡化の後で)LU分解を用いてAを次のように分解する. A = L * U ここで, Lは置換行列と下帯幅klで対角要素が1の下三角行列の積, そして, Uは上帯幅kl+kuの上三角行列である. Uの第i対角要素が0であればUは特異であり, info = iで戻る. そうでなければ, 行列Aの分解形を用いて条件数を推定する. 条件数の逆数がマシンイプシロンより小さければ警告としてinfo = n+1を返すが, 引き続き下記のように解Xを求め誤差限界を計算する. Aの分解形を用いて連立方程式を解きXを求める. 計算された解行列に反復改良を適用して精度向上を図り, その誤差限界および後退誤差推定値を計算する. 均衡化が行われた場合, 均衡化前の元の連立方程式の解を求めるために, 行列Xに左から diag(C) (Trans = "N"の場合) または diag(R) (Trans = "T"または"C"の場合) を乗ずる. 引数 [in] Fact 行列Aの分解形を入力するかどうかを指定. 入力しない場合, 行列Aを分解前に均衡化するかどうかを指定する. = "F": Afb()とIPiv()にAの分解形を入力する. Equed <> "N"であれば, R()およびC()により与えられるスケーリング係数で行列Aが均衡化済であることを示す. Ab(), Afb()およびIPiv()は変更されない. = "N": 行列AをAfb()にコピーしてから分解する. = "E": 行列Aを必要に応じて均衡化し, 次にAfb()にコピーしてから分解する. [in] Trans 連立方程式の形を指定. = "N": A * X = B. (転置なし) = "T": A^T * X = B. (転置あり) = "C": A^H * X = B. (共役転置あり) [in] N 連立方程式の数, すなわち, 行列Aの行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Kl Aの下帯幅 (Kl >= 0) [in] Ku Aの上帯幅 (Ku >= 0) [in,out] Ab() 配列 Ab(LAb1 - 1, LAb2 - 1) (LAb1 >= Kl + Ku + 1, LAb2 >= N) [in] 帯行列形式の行列 A (行1 〜 Kl+Ku+1).   Fact = "F" かつ Equed <> "N"の場合, R()および／またはC()のスケーリング係数によりAが均衡化済でなければならない. [out] Fact = "F"または"N"の場合, および, Fact = "E" かつ Equedの出力 = "N" の場合, Aは変更されない.   Fact = "E" かつ Equedの出力 <> "N" の場合, 次のように均衡化される.     Equed = "R": diag(R)*A,     Equed = "C": A*diag(C), または     Equed = "B": diag(R)*A*diag(C). [in,out] Afb() 配列 Afb(LAfb1 - 1, LAfb2 - 1) (LAfb1 >= 2Kl + Ku + 1, LAfb2 >= N) [in] Fact = "F"の場合, Dgbtrfにより計算された帯行列AのLU分解結果. Uは上帯幅Kl+Kuの上三角行列として第1〜Kl+Ku+1行に格納される. また, 分解中に使われた乗数が第Kl+Ku+2〜2Kl+Ku+1行に格納される. Equed <> "N"の場合, Afb()は均衡化済の行列Aの分解形であること. [out] Fact = "N"の場合, AのLU分解結果を返す.   Fact = "E"の場合, 均衡化済の行列AのLU分解結果を返す(均衡化の形式についてはAb()の説明を参照せよ). [in,out] IPiv() 配列 IPiv(LIPiv - 1) (LIPiv >= N) [in] Fact = "F"の場合, Dgbtrfにより計算された分解 A = L*U のピボットインデックス. 行列の第i行が第IPiv(i-1)行と交換されたことを表す. [out] Fact = "N"の場合, 行列Aの分解 A = L*U のピボットインデックス.   Fact = "E"の場合, 均衡化後の行列Aの分解 A = L*U のピボットインデックス. [in,out] Equed 均衡化の形式の指定. = "N": 均衡化なし. = "R": 行の均衡化. すなわち, diag(R)をAに左から掛ける. = "C": 列の均衡化. すなわち, diag(C)をAに右から掛ける. = "B": 行と列の均衡化. すなわち, A = diag(R)*A*diag(C) とする. [in] Fact = "F"の場合, Afb()に入力された行列の均衡化の形式を表す. [out] Fact <> "F"の場合, 適用された均衡化の形式を返す ("N", "R", "C"または"B").   Fact = "N"であれば常に"N"を返す. [in,out] R() 配列 R(LR - 1) (LR >= N) Aの行の均衡化係数diag(R). Equed = "R"または"B"の場合, diag(R)をAに左から掛ける. Fact = "N"の場合, R()はアクセスされない. [in] Fact = "F"の場合, Afb()に入力された行列の行の均衡化係数. (各要素 > 0) [out] Fact <> "F"の場合, 適用されたAの行の均衡化係数. [in,out] C() 配列 C(LC - 1) (LC >= N) Aの列の均衡化係数diag(C). Equed = "C"または"B"の場合, diag(C)をAに右から掛ける. Fact = "N"の場合, C()はアクセスされない. [in] Fact = "F"の場合, Afb()に入力された行列の列の均衡化係数. (各要素 > 0) [out] Fact <> "F"の場合, 適用された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] Equed = "N"の場合, B()は変更されない.   Trans = "N" かつ Equed = "R"または"B"の場合, B()はdiag(R)*Bで上書きされる.   Trans = "T"または"C" かつ Equed = "C"または"B"の場合, B()はdiag(C)*Bで上書きされる. [out] X() 配列 X(LX1 - 1, LX2 - 1) (LX1 >= max(1, N), LX2 >= Nrhs) (2次元配列) または X(LX - 1) (LX >= max(1, N), Nrhs = 1) (1次元配列) Info = 0 または Info = N+1 の場合, 元の連立方程式のN×Nrhs解行列 X. Equed <> "N"の場合, AおよびBは変更される. また, 均衡化された連立方程式の解は,   Trans = "N" かつ Equed = "C"または"B"であれば inv(diag(C))*X である.   Trans = "T"または"C" かつ Equed = "R"または"B"であれば inv(diag(R))*X である. [out] RCond (均衡化した場合には均衡化後の)行列Aの(1/条件数)の推定値. RCondがマシンイプシロンより小さければ(特に RCond = 0 であれば), 実用精度において行列は特異である. これは Info > 0 を返すことにより通知される. [out] FErr() 配列 FErr(LFErr - 1) (LFErr >= Nrhs) 各解ベクトルX(j)(解行列Xの第j列)の前進誤差限界. X(j)に対応する真の解をXtrueとするとき, FErr(j-1)は, (X(j) - Xtrue)の要素の最大絶対値をX(j)の要素の最大絶対値で割った値の上限の推定値である. この推定値はRCondの推定値と同程度の信頼性があり, ほぼ常に真の誤差よりも大きめに推定される. [out] BErr() 配列 BErr(LBErr - 1) (LBErr >= Nrhs) 各解ベクトルX(j)の要素に関する後退相対誤差. (すなわち, X(j)を真の解にするためのAまたはBの任意の要素の相対変化の最小値) [out] Info = 0: 正常終了. = -1: パラメータ Fact の誤り. (Fact <> "F", "N"および"E") = -2: パラメータ Trans の誤り. (Trans <> "N", "T"および"C") = -3: パラメータ N の誤り. (N < 0) = -4: パラメータ Kl の誤り. (Kl < 0) = -5: パラメータ Ku の誤り. (Ku < 0) = -6: パラメータ Ab() の誤り. = -7: パラメータ Afb() の誤り. = -8: パラメータ IPiv() の誤り. = -9: パラメータ Equed の誤り. (Fact = "F" かつ Equed <> "N", "R", "C"および "B") = -10: パラメータ R() の誤り. (R(i) <= 0 (Fact = "F" かつ Equed = "R"または"B" のときに)) = -11: パラメータ C() の誤り. (C(i) <= 0 (Fact = "F" かつ Equed = "C"または"B" のときに)) = -12: パラメータ B() の誤り. = -13: パラメータ X() の誤り. = -15: パラメータ FErr() の誤り. = -16: パラメータ BErr() の誤り. = -18: パラメータ Nrhs の誤り. (Nrhs < 0) = i (0 < i <= N): Uのi番目の対角要素が0である. 分解は完了したが, Uが特異であり, 解と誤差限界は計算できなかった. RCond = 0 を返す. = N+1: Uは非特異であるが, RCondがマシンイプシロンより小さく, これは実用精度において行列が特異であることを意味する. しかしながら, RCondの値が示すよりも計算値の精度が良いことがあるため, 解と誤差限界は計算される. [in] Nrhs (省略可) 右辺の数, すなわち, 行列Bの列数. (Nrhs >= 0) (Nrhs = 0 の場合, 処理を行わずに戻る) (省略時 = 1) [out] RPiv (省略可) ピボット成長係数の逆数 norm(A)/norm(U). ノルムは"要素最大絶対値"が使われる. これが1よりもずっと小さい場合, (均衡化した)行列AのLU分解の安定性がよくない可能性がある. 解 X, 条件数の推定値 RCond, および, 前進誤差限界 FErr も信用できない可能性があることを意味する. 0 < Info <= N で分解が失敗した場合, ピボット成長係数の逆数はAの先頭Info列分のものである. 出典 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_Zgbsvx() Const N = 3, Kl = 1, Ku = 1 Dim Ab(Kl + Ku, N - 1) As Complex, Afb(2 * Kl + Ku, N - 1) As Complex, IPiv(N - 1) As Long Dim B(N - 1) As Complex, X(N - 1) As Complex Dim R(N - 1) As Double, C(N - 1) As Double, Equed As String Dim FErr(0) As Double, BErr(0) As Double Dim RCond As Double, Info As Long Ab(0, 1) = Cmplx(0.2, -0.11): Ab(0, 2) = Cmplx(-0.93, -0.32) Ab(1, 0) = Cmplx(0.81, 0.37): Ab(1, 1) = Cmplx(-0.8, -0.92): Ab(1, 2) = Cmplx(-0.29, 0.86) Ab(2, 0) = Cmplx(0.64, 0.51): Ab(2, 1) = Cmplx(0.71, 0.59) B(0) = Cmplx(-0.0484, 0.2644): B(1) = Cmplx(-0.2644, -1.0228): B(2) = Cmplx(-0.5299, 1.5025) Call Zgbsvx("N", "N", N, Kl, Ku, Ab(), Afb(), IPiv(), Equed, R(), C(), B(), X(), RCond, FErr(), BErr(), Info) Debug.Print "X =", Debug.Print Creal(X(0)), Cimag(X(0)), Creal(X(1)), Cimag(X(1)), Creal(X(2)), Cimag(X(2)) Debug.Print "RCond =", RCond, "Equed = ", Equed Debug.Print "Info =", Info End Sub 実行結果 X = -0.15 0.19 0.2 0.94 0.79 -0.13 RCond = 0.187722560135325 Equed = N Info = 0