Dgejsv Dgesdd Dgesvd Dgesvdq Dgesvdx Dggsvd3

Dggsvd3() Sub Dggsvd3 ( Jobu As  String, Jobv As  String, Jobq As  String, M As  Long, N As  Long, P As  Long, K As  Long, L As  Long, A() As  Double, B() As  Double, Alpha() As  Double, Beta() As  Double, U() As  Double, V() As  Double, Q() As  Double, Info As  Long  ) 一般化特異値分解 (GSVD) 目的 本ルーチンはM×N実行列 A および P×N実行列 B の一般化特異値分解(GSVD)を求める. U^T*A*Q = D1*(0 R), V^T*B*Q = D2*(0 R) ここで, U, VおよびQは直交行列である. K+L = (A^T, B^T)^T の有効ランク数 とすると, R は (K+L)×(K+L)上三角行列, D1およびD2はそれぞれ M×(K+L) および P×(K+L)「対角」行列となり, 次のような構造をもつ. M-K-L >= 0 の場合: K L D1 = K (I 0) L (0 C) M-K-L (0 0) K L D2 = L (0 S) P-L (0 0) N-K-L K L (0 R) = K ( 0 R11 R12) L ( 0 0 R22) ただし, C = diag(Alpha(K), ... , Alpha(K+L-1)), S = diag(Beta(K), ... , Beta(K+L-1)), C^2 + S^2 = I. 終了時, R は A(0〜K+L-1, N-K-L〜N-1) に格納される. M-K-L < 0 の場合: K M-K K+L-M D1 = K (I 0 0) M-K (0 C 0) K M-K K+L-M D2 = M-K (0 S 0) K+L-M (0 0 I) P-L (0 0 0) N-K-L K M-K K+L-M (0 R) = K ( 0 R11 R12 R13) M-K ( 0 0 R22 R23) K+L-M ( 0 0 0 R33) ただし, C = diag(Alpha(K), ... , Alpha(M-1)), S = diag(Beta(K), ... , Beta(M-1)), C^2 + S^2 = I 終了時, (R11 R12 R13) は A(0〜M-1, N-K-L〜N-1) に格納される. また, R33 は ( 0 R22 R23) B(M-K〜L-1, N+M-K-L〜N-1) に格納される. C, S, R および必要ならば直交変換行列 U, V, Q が求められる. 行列 B がN×N行列の場合, AとBのGSVDは暗黙のうちに A*inv(B) のSVDを与える. A*inv(B) = U*(D1*inv(D2))*V^T (A^T,B^T)^T の列が正規直交であれば, AとBのGSVDはAとBのCS分解に等しい. さらに, GSVDを固有値問題を解くために使うことができる. A^T*A*X = λB^T*B*X 文献によってはAとBのGSVDを次のように表していることがある. U^T*A*X = (0 D1), V^T*B*X = (0 D2) ここで, UとVは直交行列で, Xは特異でないものとする. D1とD2は「対角」行列である. 非特異行列Xを次のようにとることにより, GSVDの前者の形式を後者の形式に変換することができる. X = Q*(I 0 ) (0 inv(R)) 引数 [in] Jobu = 'U': 直交行列 U を求める. = 'N': U を求めない. [in] Jobv ='V': 直交行列 V を求める. = 'N': V を求めない. [in] Jobq = "Q": 直交行列 Q を求める. = 'N': Q を求めない. [in] M 行列 A の行数. (M >= 0) [in] N 行列 A および B の列数. (N >= 0) [in] P 行列 B の行数. (P >= 0) [out] K [out] L KおよびLは, 「目的」に説明されているサブブロックの次数を表す. (K + L = (A^T, B^T)^T の有効ランク数) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M×N行列 A. [out] A()に三角行列R, あるいは, Rの一部が入る. 詳細は「目的」を参照せよ. [in,out] B() 配列 B(LB1 - 1, LB2 - 1) (LB1 >= P, LB2 >= N) [in] P×N行列 B. [out] M-K-L < 0 の場合, B()に三角行列Rが入る. 詳細は「目的」を参照せよ. [out] Alpha() 配列 Alpha(LAlpha - 1) (LAlpha >= N) [out] Beta() 配列 Beta(LBeta - 1) (LBeta >= N) Alpha()およびBeta()にAおよびBの一般化特異値が入る.   Alpha(0 to K-1) = 1,   Beta(0 to K-1) = 0, M-K-L >= 0 の場合,   Alpha(K to K+L-1) = C,   Beta(K to K+L-1) = S, M-K-L < 0 の場合,   Alpha(K to M-1)= C, Alpha(M to K+L-1) = 0,   Beta(K to M-1) = S, Beta(M to K+L-1) = 1, また   Alpha(K+L to N-1) = 0,   Beta(K+L to N-1) = 0. [out] U() 配列 U(LU1 - 1, LU2 - 1) (LU1 >= M, LU2 >= M) Jobu = 'U': U()にM×M直交行列Uが入る. Jobu = 'N': U()は参照されない. [out] V() 配列 V(LV1 - 1, LV2 - 1) (LV1 >= P, LV2 >= P) Jobv = 'V': V()にP×P直交行列Vが入る. Jobv = 'N': V()は参照されない. [out] Q() 配列 Q(LQ1 - 1, LQ2 - 1) (LQ1 >= N, LQ2 >= N) Jobq = 'Q': Q()にN×N直交行列Qが入る. Jobq = 'N': Q()は参照されない. [out] Info = 0: 正常終了. = -1: パラメータ Jobu の誤り. (Jobu <> "U" および "N") = -2: パラメータ Jobv の誤り. (Jobv <> "V" および "N") = -3: パラメータ Jobq の誤り. (Jobq <> "Q" および "N") = -4: パラメータ M の誤り. (M < 0) = -5: パラメータ N の誤り. (N < 0) = -6: パラメータ P の誤り. (P < 0) = -9: パラメータ A() の誤り. = -10: パラメータ B() の誤り. = -11: パラメータ Alpha() の誤り. = -12: パラメータ Beta() の誤り. = -13: パラメータ U() の誤り. = -14: パラメータ V() の誤り. = -15: パラメータ Q() の誤り. = 1: ヤコビ型手続きが収束しなかった. 出典 LAPACK 使用例 行列 A および 行列 B の一般化特異値分解(GSVD)を求める. ただし, ( 1 6 11 ) ( 1 5 9 ) ( 2 7 12 ) ( 2 6 10 ) A = ( 3 8 13 ), B = ( 3 7 11 ) ( 4 9 14 ) ( 4 8 12 ) ( 5 10 15 ) である. Sub Ex_Dggsvd3() Const M = 5, N = 3, P = 4 Dim A(M - 1, N - 1) As Double, B(P - 1, N - 1) As Double Dim Alpha(N - 1) As Double, Beta(N - 1) As Double Dim U(M - 1, M - 1) As Double, V(P - 1, P - 1) As Double, Q(N - 1, N - 1) As Double Dim K As Long, L As Long, Info As Long A(0, 0) = 1: A(0, 1) = 6: A(0, 2) = 11 A(1, 0) = 2: A(1, 1) = 7: A(1, 2) = 12 A(2, 0) = 3: A(2, 1) = 8: A(2, 2) = 13 A(3, 0) = 4: A(3, 1) = 9: A(3, 2) = 14 A(4, 0) = 5: A(4, 1) = 10: A(4, 2) = 15 B(0, 0) = 1: B(0, 1) = 5: B(0, 2) = 9 B(1, 0) = 2: B(1, 1) = 6: B(1, 2) = 10 B(2, 0) = 3: B(2, 1) = 7: B(2, 2) = 11 B(3, 0) = 4: B(3, 1) = 8: B(3, 2) = 12 Call Dggsvd3("U", "V", "Q", M, N, P, K, L, A(), B(), Alpha(), Beta(), U(), V(), Q(), Info) Debug.Print "Alpha =", Alpha(0), Alpha(1), Alpha(2) Debug.Print "Beta =", Beta(0), Beta(1), Beta(2) Debug.Print "R =" Debug.Print A(0, 0), A(0, 1), A(0, 2) Debug.Print A(1, 0), A(1, 1), A(1, 2) Debug.Print A(2, 0), A(2, 1), A(2, 2) Debug.Print "U =" Debug.Print U(0, 0), U(0, 1), U(0, 2), U(0, 3), U(0, 4) Debug.Print U(1, 0), U(1, 1), U(1, 2), U(1, 3), U(1, 4) Debug.Print U(2, 0), U(2, 1), U(2, 2), U(2, 3), U(2, 4) Debug.Print U(3, 0), U(3, 1), U(3, 2), U(3, 3), U(3, 4) Debug.Print U(4, 0), U(4, 1), U(4, 2), U(4, 3), U(4, 4) Debug.Print "V =" Debug.Print V(0, 0), V(0, 1), V(0, 2), V(0, 3) Debug.Print V(1, 0), V(1, 1), V(1, 2), V(1, 3) Debug.Print V(2, 0), V(2, 1), V(2, 2), V(2, 3) Debug.Print V(3, 0), V(3, 1), V(3, 2), V(3, 3) Debug.Print "Q =" Debug.Print Q(0, 0), Q(0, 1), Q(0, 2) Debug.Print Q(1, 0), Q(1, 1), Q(1, 2) Debug.Print Q(2, 0), Q(2, 1), Q(2, 2) Debug.Print "K =", K, "L =", L, "Info =", Info End Sub 実行結果 Alpha = 0.802331750723485 0.826878800377287 0 Beta = 0.596878347555838 0.56238016455652 0 R = 0 -5.40364104666407 -35.7382954334068 0 0 24.1572950255889 0 0 0 U = -9.72172871075144E-02 -0.768471729530407 -0.369047538143976 -0.45043765685718 -0.246799173164912 0.116318998962025 -0.535228820674366 1.34135766249015E-02 0.516457650749435 0.658096931268109 0.329855285031565 -0.301985911818325 0.857835814511989 -0.145130339105955 -0.207496746990269 0.543391571101105 -6.87430029622829E-02 -0.279722206322779 0.542638353392326 -0.572100607164143 0.756927857170645 0.164499905893758 -0.222479646670135 -0.463528008178625 0.368299596051215 V = -4.71369494280112E-02 -0.835331136734781 0.537498761722548 0.105333191097239 0.239901343964157 -0.492389424301732 -0.795175740823531 0.260183668214875 0.526939637356325 -0.149447711868684 -2.21448035205818E-02 -0.836366909721466 0.813977930748493 0.193494000564365 0.279821782621565 0.470850050409352 Q = 0.408248290463863 -0.906218690095711 -0.1100046319686 -0.816496580927726 -0.308596432505139 -0.487956530009011 0.408248290463863 0.289025825085432 -0.865908428049421 K = 0 L = 2 Info = 0