Zgejsv Zgesdd Zgesvd Zgesvdq Zgesvdx Zggsvd3

Zggsvd3() Sub Zggsvd3 ( 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  Complex, B() As  Complex, Alpha() As  Double, Beta() As  Double, U() As  Complex, V() As  Complex, Q() As  Complex, Info As  Long  ) 一般化特異値分解 (GSVD) (複素行列) 目的 本ルーチンはM×N複素行列 A および P×N複素行列 B の一般化特異値分解(GSVD)を求める. U^H*A*Q = D1*(0 R), V^H*B*Q = D2*(0 R) ここで, U, VおよびQはユニタリ行列である. K+L = (A^H, B^H)^H の有効ランク数 とすると, 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^H (A^H,B^H)^H の列が正規直交であれば, AとBのGSVDはAとBのCS分解に等しい. さらに, GSVDを固有値問題を解くために使うことができる. A^H*A*X = λB^H*B*X 文献によってはAとBのGSVDを次のように表していることがある. U^H*A*X = (0 D1), V^H*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^H, B^H)^H の有効ランク数) [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)を求める. ただし, ( 0.20-0.11i -0.93-0.32i 0.81+0.37i ) A = ( -0.80-0.92i -0.29+0.86i 0.64+0.51i ) ( 0.71+0.59i -0.15+0.19i 0.20+0.94i ) ( 0.57-0.91i -0.28-0.45i 0.25+0.91i ) B = ( 0.83+0.46i 0.63-0.19i -0.69+0.09i ) ( 0.24-1.33i -0.56-0.67i 0.90+1.25i ) とする. Sub Ex_Zggsvd3() Const M = 3, N = 3, P = 3 Dim A(M - 1, N - 1) As Complex, B(P - 1, N - 1) As Complex Dim Alpha(N - 1) As Double, Beta(N - 1) As Double Dim U(M - 1, M - 1) As Complex, V(P - 1, P - 1) As Complex, Q(N - 1, N - 1) As Complex Dim K As Long, L As Long, 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, 0) = cmplx(0.57, -0.91): B(0, 1) = cmplx(-0.28, -0.45): B(0, 2) = cmplx(0.25, 0.91) B(1, 0) = cmplx(0.83, 0.46): B(1, 1) = cmplx(0.63, -0.19): B(1, 2) = cmplx(-0.69, 0.09) B(2, 0) = cmplx(0.24, -1.33): B(2, 1) = cmplx(-0.56, -0.67): B(2, 2) = cmplx(0.9, 1.25) Call Zggsvd3("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 Creal(A(0, 0)), Cimag(A(0, 0)), Creal(A(0, 1)), Cimag(A(0, 1)) Debug.Print Creal(A(1, 0)), Cimag(A(1, 0)), Creal(A(1, 1)), Cimag(A(1, 1)) Debug.Print Creal(A(2, 0)), Cimag(A(2, 0)), Creal(A(2, 1)), Cimag(A(2, 1)) Debug.Print Creal(A(0, 2)), Cimag(A(0, 2)) Debug.Print Creal(A(1, 2)), Cimag(A(1, 2)) Debug.Print Creal(A(2, 2)), Cimag(A(2, 2)) Debug.Print "U =" Debug.Print Creal(U(0, 0)), Cimag(U(0, 0)), Creal(U(0, 1)), Cimag(U(0, 1)) Debug.Print Creal(U(1, 0)), Cimag(U(1, 0)), Creal(U(1, 1)), Cimag(U(1, 1)) Debug.Print Creal(U(2, 0)), Cimag(U(2, 0)), Creal(U(2, 1)), Cimag(U(2, 1)) Debug.Print Creal(U(0, 2)), Cimag(U(0, 2)) Debug.Print Creal(U(1, 2)), Cimag(U(1, 2)) Debug.Print Creal(U(2, 2)), Cimag(U(2, 2)) Debug.Print "V =" Debug.Print Creal(V(0, 0)), Cimag(V(0, 0)), Creal(V(0, 1)), Cimag(V(0, 1)) Debug.Print Creal(V(1, 0)), Cimag(V(1, 0)), Creal(V(1, 1)), Cimag(V(1, 1)) Debug.Print Creal(V(2, 0)), Cimag(V(2, 0)), Creal(V(2, 1)), Cimag(V(2, 1)) Debug.Print Creal(V(0, 2)), Cimag(V(0, 2)) Debug.Print Creal(V(1, 2)), Cimag(V(1, 2)) Debug.Print Creal(V(2, 2)), Cimag(V(2, 2)) Debug.Print "Q =" Debug.Print Creal(Q(0, 0)), Cimag(Q(0, 0)), Creal(Q(0, 1)), Cimag(Q(0, 1)) Debug.Print Creal(Q(1, 0)), Cimag(Q(1, 0)), Creal(Q(1, 1)), Cimag(Q(1, 1)) Debug.Print Creal(Q(2, 0)), Cimag(Q(2, 0)), Creal(Q(2, 1)), Cimag(Q(2, 1)) Debug.Print Creal(Q(0, 2)), Cimag(Q(0, 2)) Debug.Print Creal(Q(1, 2)), Cimag(Q(1, 2)) Debug.Print Creal(Q(2, 2)), Cimag(Q(2, 2)) Debug.Print "K =", K, "L =", L, "Info =", Info End Sub 実行結果 Alpha = 0.430704325811232 0.976338298760024 0.991145496564202 Beta = 0.902493093451408 0.216248760399642 0.132780287093004 R = -2.49290056457435 0 0.457594814366194 0.194431879391052 0 0 -1.26470568603387 0 0 0 0 0 1.6180003034571 1.39174354257137 -7.00070473318143E-02 -0.174007882105091 -1.75471215104157 0 U = 0.568161065126274 0.474196630652255 -0.180338037878126 -6.24249232440974E-03 0.201921419044388 0.445329052169583 0.137381224968446 -0.594162404216661 -0.365049262631842 -0.282806243264565 -0.454409702123437 -0.623737510850184 0.622382357104722 0.180027732589302 -0.33819653446823 -0.524051826005947 0.431431460509124 -7.13435328533903E-02 V = -0.259639196635024 0.401093020738518 4.88190028247251E-02 3.25673943950224E-02 -0.400793045759536 -0.226433378198375 -0.714410691722334 0.432401521895175 -0.122502964271075 0.738104165913481 -0.437440072643577 -0.328403014673516 0.331492445838016 -0.811406620889563 -0.275482203370116 -0.121855253273566 2.47600244928659E-02 0.374672988073607 Q = 0.498168838100246 0.295619901645173 0.262656203742932 -0.420386203776387 0.659025789544342 -0.110123110678936 0.347354481147794 -5.02490827933195E-02 -0.285036186254502 -0.369795841286094 0.793403710637756 4.02194525204689E-02 -0.644893096482004 5.32610074710518E-02 0.613682926576038 -0.231885896388499 -6.10928315913905E-02 0.383622238414744 K = 0 L = 3 Info = 0