Dgssvd Zgssvd

Dgssvd() Sub Dgssvd ( Jobv As  String, Jobu As  String, M As  Long, N As  Long, Val() As  Double, Ptr() As  Long, Ind() As  Long, Which As  String, Nev As  Long, Ncv As  Long, S() As  Double, V() As  Double, U() As  Double, Optional Info As  Long, Optional Nconv As  Long, Optional Niter As  Long, Optional Base As  Long = -1, Optional Format As  Long = 0, Optional Maxiter As  Long = 300, Optional Tol As  Double = 0  ) 特異値分解 (SVD) (一般疎行列) (リスタート付きランチョス法 (IRLM)) (Arpack) (ドライバ) 目的 行列 A の特異値(s)と対応する右特異ベクトル(v)を, 次の対称固有値問題を解くことによって求める. (A^T*A)*v = s^2*v A は M x N 実一般疎行列である. 必要であれば, 左特異ベクトルも u = A*v/s によって求められる. 疎行列 A は CSC または CSR 形式で格納されているものとする. 対称固有値問題は, Arpack ルーチン Dsaupd および Dseupd を使用し, リスタート付きランチョス法 (Implicitly restarted Lanczos method (IRLM)) により解く. 引数 [in] Jobv = 'N': 右特異ベクトルを計算しない. = 'V': 右特異ベクトルを計算する. [in] Jobu = 'N': 左特異ベクトルを計算しない. = 'V': 左特異ベクトルを計算する. (Jobv も 'V' でなければならない.) [in] M 行列の行数. (M >= 0) (M = 0 の場合, 処理を行わずに戻る) [in] N 行列の列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) 入力行列の非ゼロ要素の値. (Nnz は非ゼロ要素数) [in] Ptr() 配列 Ptr(LPtr - 1) (LPtr >= N + 1 入力行列の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind() 配列 Ind(LInd - 1) (LInd >= Nnz) 入力行列の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). (Nnz は非ゼロ要素数) [in] Which = "LA": 最大のものから Nev 個の特異値を求める. = "SA": 最小のものから Nev 個の特異値を求める. = "LM": 絶対値最大のものから Nev 個の特異値を求める. = "SM": 絶対値最小のものから Nev 個の特異値を求める. = "BE": スペクトルの両端から半分づつ, 合計 Nev 個の特異値を求める. Nev が奇数の場合, 大きい側を1個多く求める. [in] Nev 求める特異値の数. (0 < Nev < N) [in] Ncv 各反復において生成されるランチョス基底ベクトルの数. (Nev < Ncv <= N) (Ncv >= 2*Nev が推奨される) [out] S() 配列 S(LS - 1) (LS >= Nev) A の特異値を返す. 値は昇順に格納される. [out] V() 配列 V(LV1 - 1, LV2 - 1) (LV1 >= N, LV2 >= Nev) Nev本の右特異ベクトル. Jobv = 'N' の場合, 参照されない. [out] U() 配列 U(LU1 - 1, LU2 - 1) (LU1 >= M, LU2 >= Nev) Nev本の左特異ベクトル. Jobu = 'N' の場合, 参照されない. [out] Info リターンコード. = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: 最大反復回数に達した. = 3: リスタート付きランチョス反復中にシフト値が適用できなかった. Ncv を Nev に比較して大きくすると改善するかもしれない (Ncv >= 2*Nev が推奨される). = 11: 初期残差ベクトルが 0 である. = 12: ランチョス分解に失敗した. = 13: LAPACK ルーチンで固有値の計算に失敗した. = 14: Dsaupd は十分な精度の固有値を得られなかった. = 15 〜 17: 内部処理エラー. [out] Nconv (省略可) 収束条件を満たしたリッツ値の数. [out] Niter (省略可) ランチョス反復回数. [in] Base (省略可) Ptr() および Ind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Ptr(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Format (省略可) 行列の格納形式. (省略時 = 0) = 0: CSR 形式. = 1: CSC 形式. [in] Maxiter (省略可) ランチョス反復の最大回数. (MaxIter > 0) (省略時 = 300) [in] Tol (省略可) 停止基準: リッツ値の最大相対誤差. (省略時 = 0) Tol <= 0 の場合, マシン精度とみなす. 使用例 行列 A の大きい方から2個の特異値およびその特異ベクトルを求める. ただし, ( 1 6 11 ) ( 2 7 12 ) A = ( 3 8 13 ) ( 4 9 14 ) ( 5 10 15 ) とする. Sub Ex_Dgssvd() Const M = 5, N = 3, Nnz = M * N, Ncv = 3, Nev = 2 Dim A(Nnz - 1) As Double, Ia(M) As Long, Ja(Nnz - 1) As Long Dim S(Nev - 1) As Double, V(N - 1, Nev - 1) As Double, U(M - 1, Nev - 1) As Double Dim Nconv As Long, Niter As Long, Info As Long, J As Long A(0) = 1: A(1) = 6: A(2) = 11: A(3) = 2: A(4) = 7: A(5) = 12: A(6) = 3: A(7) = 8: A(8) = 13: A(9) = 4: A(10) = 9: A(11) = 14: A(12) = 5: A(13) = 10: A(14) = 15 Ia(0) = 0: Ia(1) = 3: Ia(2) = 6: Ia(3) = 9: Ia(4) = 12: Ia(5) = 15 Ja(0) = 0: Ja(1) = 1: Ja(2) = 2: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2: Ja(6) = 0: Ja(7) = 1: Ja(8) = 2: Ja(9) = 0: Ja(10) = 1: Ja(11) = 2: Ja(12) = 0: Ja(13) = 1 Ja(14) = 2 Call Dgssvd("V", "V", M, N, A(), Ia(), Ja(), "LM", Nev, Ncv, S(), V(), U(), Info, Nconv) Debug.Print "S =" Debug.Print S(0), S(1) Debug.Print "U =" Debug.Print U(0, 0), U(0, 1) Debug.Print U(1, 0), U(1, 1) Debug.Print U(2, 0), U(2, 1) Debug.Print U(3, 0), U(3, 1) Debug.Print U(4, 0), U(4, 1) Debug.Print "V =" Debug.Print V(0, 0), V(0, 1) Debug.Print V(1, 0), V(1, 1) Debug.Print V(2, 0), V(2, 1) Debug.Print "Nconv =" + Str(Nconv) + ", Info =" + Str(Info) End Sub 実行結果 S = 2.46539669691652 35.1272233335747 U = 0.688686643768252 0.354557057037681 0.375554529395872 0.398696369998832 6.24224150234914E-02 0.442835682959984 -0.250709699348889 0.486974995921135 -0.563841813721269 0.531114308882286 V = -0.890317132783019 0.201664911192694 -0.257331626824051 0.516830501392305 0.375653879134918 0.831996091591915 Nconv = 2, Info = 0