Zgecov Zgecovs Zgecovy Zgels Zgelsd Zgelss Zgelsy Zgetsls

Zgelsy() Sub Zgelsy ( M As  Long, N As  Long, A() As  Complex, B() As  Complex, Jpvt() As  Long, RCond As  Double, Rank As  Long, Info As  Long, Optional Nrhs As  Long = 1  ) 優決定または劣決定系連立一次方程式 Ax = b の解 (完全直交分解) (複素行列) 目的 本ルーチンは複素線形最小二乗問題の最小ノルム解を, Aの完全直交分解を使って求める. || A * X - B || を最小化する. AはM×N行列で, ランク落ちしていてもよい. いくつかの右辺ベクトル b および解ベクトル x を1回の呼び出しで扱うことができる. これらのベクトルは, M×Nrhs右辺行列BおよびN×Nrhs解行列Xの列として格納される. まず最初に列のピボット選択付きQR分解を求める. A * P = Q * [ R11 R12 ] [ 0 R22 ] R11は最大の主小行列を示し, その推定条件数は1/rcondより小さい. R11の次数(rank)はAの有効ランク数である. そして, R22は無視してよいと考えられ, R12は右からの直交変換により消え, 次の完全直交分解になる. A * P = Q * [ T11 0 ] * Z [ 0 0 ] これより, 最小ノルム解は次のように求められる. X = P * Z^H [ T11^(-1)*Q1^H*B ] [ 0 ] ここで, Q1はQのはじめのrank列よりなる. 引数 [in] M 行列 A の行数. (M >= 0) (M = 0 の場合, Rank = 0 として戻る) [in] N 行列 A の列数. (N >= 0) (N = 0 の場合, Rank = 0 として戻る) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M×N行列 A. [out] A()は完全直交分解により上書きされる. [in,out] B() 配列 B(LB1 - 1, LB2 - 1) (LB1 >= max(M, N), LB2 >= Nrhs) (2次元配列) または B(LB - 1) (LB >= max(M, N), Nrhs = 1) (1次元配列) [in] M×Nrhs右辺行列 B. [out] N×Nrhs解行列 X. [in,out] Jpvt() 配列 Jpvt(LJpvt - 1) (LJpvt >= N) [in] Jpvt(i-1) <> 0 のとき, Aのi列をA*Pの前の方に持ってくる. 0のとき, i列はフリーな列である. [out] Jpvt(i) = k であれば, A*Pのi列はAのk列だったことを示す. [in] RCond RCondはAの有効ランク数を決めるために使われる. 有効ランク数は, Aのピボット選択付きQR分解の最大の主小行列R11の次数である. R11の推定条件数 < 1/RCond である. [out] Rank Aの有効ランク数, すなわち, 部分行列R11の次数. これはAの完全直交分解における部分行列T11の次数に等しい. [out] Info = 0: 正常終了. = -1: パラメータ M の誤り. (M < 0) = -2: パラメータ N の誤り. (N < 0) = -3: パラメータ A() の誤り. = -4: パラメータ B() の誤り. = -5: パラメータ Jpvt() の誤り. = -9: パラメータ Nrhs の誤り. (Nrhs < 0) [in] Nrhs (省略可) 右辺の数, すなわち, 行列BおよびXの列数. (Nrhs >= 0) (Nrhs = 0 の場合, Rank = 0 として戻る) (省略時 = 1) 出典 LAPACK 使用例 優決定系連立1次方程式 Ax = B の最小二乗解を求める. また, 分散を求める. ただし, ( -0.82+0.83i 0.18-0.94i -0.18-0.12i ) A = ( -0.76-0.24i 0.57-0.16i -0.08-0.27i ) ( 1.90+0.26i -0.98+0.54i 0.21+0.28i ) ( 0.50-0.30i -0.31+0.37i 0.22+0.19i ) ( 1.7126-0.6648i ) B = ( 0.8697+0.7604i ) ( -2.1048-1.6171i ) ( -0.9297+0.1252i ) とする. Sub Ex_Zgelsy() Const M = 4, N = 3 Dim A(M - 1, N - 1) As Complex, B(M - 1) As Complex, Ci(N - 1) As Complex Dim Jpvt(N - 1) As Long, RCond As Double, Rank As Long, Info As Long, I As Long A(0, 0) = Cmplx(-0.82, 0.83): A(0, 1) = Cmplx(0.18, -0.94): A(0, 2) = Cmplx(-0.18, -0.12) A(1, 0) = Cmplx(-0.76, -0.24): A(1, 1) = Cmplx(0.57, -0.16): A(1, 2) = Cmplx(-0.08, -0.27) A(2, 0) = Cmplx(1.9, 0.26): A(2, 1) = Cmplx(-0.98, 0.54): A(2, 2) = Cmplx(0.21, 0.28) A(3, 0) = Cmplx(0.5, -0.3): A(3, 1) = Cmplx(-0.31, 0.37): A(3, 2) = Cmplx(0.22, 0.19) B(0) = Cmplx(1.7126, -0.6648): B(1) = Cmplx(0.8697, 0.7604) B(2) = Cmplx(-2.1048, -1.6171): B(3) = Cmplx(-0.9297, 0.1252) For I = 0 To N - 1 Jpvt(I) = 0 Next RCond = 0.0001 Call Zgelsy(M, N, A(), B(), Jpvt(), RCond, Rank, Info) If Info <> 0 Then Debug.Print "Error in Zgelsy: Info =", Info Exit Sub End If Debug.Print "X =" Debug.Print Creal(B(0)), Cimag(B(0)), Creal(B(1)), Cimag(B(1)) Debug.Print Creal(B(2)), Cimag(B(2)) Call Zgecovy(0, N, A(), Jpvt(), Ci(), Info) Debug.Print "Var =" Debug.Print Creal(Ci(0)), Creal(Ci(1)), Creal(Ci(2)) Debug.Print "Rank =", Rank, "Info =", Info End Sub 実行結果 X = -0.819999999999999 -0.939999999999999 0.740000000000001 0.200000000000002 0.480000000000002 0.21 Var = 5.46169501938982 15.1880504061464 21.4241290120714 Rank = 3 Info = 0