Zgejsv Zgesdd Zgesvd Zgesvdq Zgesvdx Zggsvd3

Zgesvdq() Sub Zgesvdq ( Joba As  String, Jobp As  String, Jobr As  String, Jobu As  String, Jobv As  String, M As  Long, N As  Long, A() As  Complex, S() As  Double, U() As  Complex, V() As  Complex, Numrank As  Long, Info As  Long, Optional Cond As  Double  ) 特異値分解 (SVD) (前処理付きQR法) (複素行列) 目的 本ルーチンは M x N 複素行列 A の特異値分解(SVD)を求める. ただし, M >= N である. A の SVD は次のように表される. [++] [xx] [x0] [xx] A = U * SIGMA * V^H, [++] = [xx] * [ox] * [xx] [++] [xx] ただし, SIGMA は N x N 対角行列, U は M x N 正規直交行列, V は N x N ユニタリ行列である. SIGMA の対角要素は A の特異値である. U および V の列は A のそれぞれ右および左特異ベクトルである. 引数 [in] Joba 計算された SVD の精度のレベルを指定する. = "A": 要求精度は || delta A ||_F <= f(M, N) * EPS * || A ||_F で上限が示される後退誤差を持つことに対応する. ただし, EPS はマシンイプシロン (= Dlamch("E")) である. これにより, Rank revealing QR分解で計算された三角行列が EPS * ||A||_F のオーダーの基準値以下であれば, このルーチンがそれを切り捨てることができるようにする (積極的な切り捨てレベル). = "M": "A" と同様であるが, 切り捨てはより穏やかで, QR分解の三角行列の対角要素にドロップがある場合にのみ許される (中間的な切り捨てレベル). = "H": 高精度を要求する. Rank revealing QR分解に基づくランクの決定を行わない. = "E": "H" と同様であるが, 追加で列スケーリングした A の条件数の推定を行い Cond に返す. N^(-1/4)*Cond <= ||pinv(A_scaled)||_2 <= N^(1/4)*Cond. [in] Jobp = "P": A の行は ||A(i,:)||_infinity に関して降順に並べ替えられる. これにより, 余分なデータ移動分のコストをもって数値精度を改善することができる. 数値的安定性のために推奨される. = "N": 行の並べ替えを行わない. [in] Jobr = "T": 初期並べ替え付きQR分解の後, 求められた三角行列 R の転置 R^H に Zgesvd を適用する. 余分なデータ移動(行列の転置)を含む. 実験, 研究, 開発に有用である. = "N": 三角行列 R を Zgesvd に入力する. 余分なデータ移動が少ない. [in] Jobu = "A": M 本の左特異ベクトルすべてを計算し, 配列 U()に返す. U() の項を参照せよ. = "S" または "U": min(M, N) (= N) 本の左特異ベクトルを計算し, 配列 U()に返す. U() の項を参照せよ. = "R": ランク数 Numrank を決定し Numrank 本の左特異ベクトルだけを計算し, 配列 U()に返す. = "N": 左特異ベクトルを計算しない. 注 - Jobu = "F" はサポートされない. [in] Jobv = "A" または "V": N 本の右特異ベクトルすべてを計算し, 配列 V()に返す. = "R": ランク数 Numrank を決定し Numrank 本の右特異ベクトルだけを計算し, 配列 V()に返す. このオプションは Jobu = "R" または Jobu = "N" の場合にだけ許される. それ以外であれば不正である. = "N": 右特異ベクトルを計算しない. [in] M 入力行列 A の行数. (M >= 0) (M = 0 の場合, 処理を行わずに戻る) [in] N 入力行列 A の列数. (M >= N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in,out] A() 配列 A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] 入力行列 A. [out] Jobu <> "N" または Jobv <> "N" の場合, A() の下三角部分にハウスホルダーベクトルが Zgeqp3 によって格納される. Jobu = "F" の場合, これらのハウスホルダーベクトルと work[0:N-1] を初期の A の並べ替え付きQR分解の Q を復元するために使用することができる. U() の項を参照せよ. [out] S() 配列 S(LS - 1) (LS >= N) A の特異値 (S(i) >= S(i+1) となるように並んでいる). [out] U() 配列 U(LU1 - 1, LU2 - 1) (LU1 >= M (Jobu = "A", "S", "U" または "R" の場合)) (LU2 >= M (Jobu = "A" の場合), LU2 >= N (Jobu = "S", "U" または "R" の場合)) Jobu = "A": U() に M 本の左特異ベクトルを返す. Jobu = "S", "U" または "R": U() に先頭から N 本または Numrank 本の左特異ベクトルを返す. Jobu = "N": U() は参照されない. [out] V() 配列 V(LV1 - 1, LV2 - 1) (LV1 >= N (Jobv = "A", "V" または "R", または Joba = "E" の場合)) (LV2 >= N (Jobv = "A", "V" または "R", または Joba = "E" の場合)) Jobv = "A" または "V": V() に N x N ユニタリ行列 V^H を返す. Jobv = "R": V() に V^H (大きい方 Numrank 個の特異値の右特異ベクトルで, 行ごとに格納される) の最初の Numrank 行を返す. Jobv = "N" かつ Joba = "E": V() は作業領域として使用される. Jobv = "N" かつ Joba <> "E": V() は参照されない. [out] Numrank Joba の指定に従って Rank revealing QR分解の後, 最初に決定されたランク数. Jobv = "R" かつ Jobu = "R" の場合, 先頭から Numrank 個の特異値およびベクトルのみが Zgesvd の呼び出しにおいて要求される. 計算された特異値に 0 のものがあれば Numrank の最終的な値は減るかもしれない. [out] Info = 0: 正常終了 = -1: パラメータ Joba の誤り (Joba <> "A", "M", "H" および "E") = -2: パラメータ Jobp の誤り (Jobp <> "P" および "N") = -3: パラメータ Jobr の誤り (Jobr <> "T" および "N") = -4: パラメータ Jobu の誤り (Jobu <> "A", "S", "R" および "N") = -5: パラメータ Jobv の誤り (Jobv <> "A", "V", "R" および "N") = -6: パラメータ M の誤り (M < 0) = -7: パラメータ N の誤り (N < 0 or N > M) = -8: パラメータ A() の誤り. = -9: パラメータ S() の誤り. = -10: パラメータ U() の誤り. = -11: パラメータ V() の誤り. > 0: Zbdsqr が収束しなかった. 中間結果の二重対角形 B (Zgesvd で計算される) の上副対角要素のうち Info 個が 0 に収束しなかった. [out] Cond (省略可) Joba = "E" の場合, 列スケーリングした A の条件数の推定値を返す. 出典 LAPACK 使用例 行列Aの特異値, 左および右特異ベクトルを求める. ただし, ( 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 ) である. Sub Ex_Zgesvdq() Const M = 3, N = 3 Dim A(M - 1, N - 1) As Complex, V(N - 1, N - 1) As Complex, U(M - 1, N - 1) As Complex Dim S(N - 1) As Double, Numrank 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) Call Zgesvdq("A", "P", "N", "S", "A", M, N, A(), S(), U(), V(), Numrank, Info) Debug.Print "Singular values =", S(0), S(1), S(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^H =" 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 "Numrank =", Numrank, "Info =", Info End Sub 実行結果 Singular values = 2.07084030821889 1.23513084760163 0.901483337149816 U = 0.322093837376809 -0.402732043434635 9.95958553097357E-02 0.309610800223195 -0.255722271568769 -0.689932737910445 -0.631621088952191 0.037032935734802 0.438560869464762 -1.80488745428177E-02 -7.00863751874239E-02 0.699280401316292 -0.71181974068539 0.348707521966053 9.32902099996814E-02 -0.222663951646943 0.551156386049191 9.82856061041462E-02 V^H = 0.603023132575377 0 -0.36655126147564 -0.394522570687146 0.66381846531382 -0 0.134961467338043 -0.149629139353428 0.44238913491109 -0 0.297134277147366 0.762299060909463 -0.16075433957414 0.566138903287951 0.366911928609536 -0.61977189734833 -0.331437452593181 0.15827959884907 Numrank = 3 Info = 0