Dgsev Dgsgv Dnaupd Dneupd Dsaupd Dseupd Dssev Dssgv

Dssgv() Sub Dssgv ( Jobz As  String, Uplo As  String, Uplo2 As  String, N As  Long, Val() As  Double, Ptr() As  Long, Ind() As  Long, Val2() As  Double, Ptr2() As  Long, Ind2() As  Long, Which As  String, Nev As  Long, Ncv As  Long, D() As  Double, Z() As  Double, Optional Info As  Long, Optional Nconv As  Long, Optional Niter As  Long, Optional Base As  Long = -1, Optional Base2 As  Long = -1, Optional Format As  Long = 0, Optional Mode As  Long = 2, Optional Sigma As  Double = 0, Optional Maxiter As  Long = 300, Optional Tol As  Double = 0  ) 一般化固有値問題 (対称疎行列) (リスタート付きランチョス法 (IRLM)) (Arpack) (ドライバ) 目的 一般化固有値問題 A*z = λ*B*z において, λは固有値, z はその固有ベクトルを表す. A は対称疎行列, B は半正定値対称疎行列である. 本ルーチンは, Arpack ルーチン Dsaupd および Dseupd を使用して, リスタート付きランチョス法 (Implicitly restarted Lanczos method (IRLM)) により上の問題の固有値 λ および固有ベクトル z を求める. 対称疎行列 A および B は上三角部分のみ, 下三角部分のみあるいは両方(行列全体)が CSC 形式または CSR 形式で格納されているものとする. 本ルーチンは以下の2つのモードで動作する. (1) 通常インバースモード Mode = 2 であれば, 通常インバースモードで計算を行う. 大きい方から, あるいは, 小さい方から数個の固有値およびその固有ベクトルを求めることができる. Sigma は参照されない. (2) シフト＆インバートモード Mode = 3 であれば, シフト＆インバートモードで計算を行う. 指定した σ (シフト) の近くの固有値を求めることができる. OP = (A - σ*M)^(-1)*M, M = B として, 別の固有値問題 OP*x = ν*M*x の固有値 ν を求めたとすると, 元の問題の固有値 λ は ν = 1/(λ - σ) より求めることができる. すなわち, 最大の ν を求めれば σ に最も近い λ を求められる. パラメータ Which は求める固有値 ν を指定する. 一般的には, "LM"とすれば σ の最も近くの固有値 λ を求めることができる. 引数 [in] Jobz = 'N': 固有値のみ求める. = 'V': 固有値と固有ベクトルを求める. [in] Uplo 行列 A の上三角部分あるいは下三角部分のどちらが格納されているかを指定. = 'U': 上三角部分. = 'L': 下三角部分. = 'B': 行列全体(上および下三角部分). 指定された以外の三角部分(対角成分をを除く)は無視される. [in] Uplo2 行列 B の上三角部分あるいは下三角部分のどちらが格納されているかを指定. = 'U': 上三角部分. = 'L': 下三角部分. = 'B': 行列全体(上および下三角部分). 指定された以外の三角部分(対角成分をを除く)は無視される. [in] N 行列の行および列数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) 行列 A の非ゼロ要素の値. (Nnz は非ゼロ要素数) [in] Ptr() 配列 Ptr(LPtr - 1) (LPtr >= N + 1 行列 A の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind() 配列 Ind(LInd - 1) (LInd >= Nnz) 行列 A の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). (Nnz は非ゼロ要素数) [in] Val2() 配列 Val2(LVal2 - 1) (LVal2 >= Nnz2) 行列 B の非ゼロ要素の値. (Nnz2 は非ゼロ要素数) [in] Ptr2() 配列 Ptr2(LPtr2 - 1) (LPtr2 >= N + 1 行列 B の列ポインタ(CSC形式の場合)または行ポインタ(CSR形式の場合). [in] Ind2() 配列 Ind2(LInd2 - 1) (LInd2 >= Nnz2) 行列 B の行インデクス(CSC形式の場合)または列インデクス(CSR形式の場合). (Nnz2 は非ゼロ要素数) [in] Which = "LA": 最大のものから Nev 個の固有値を求める. = "SA": 最小のものから Nev 個の固有値を求める. = "LM": 絶対値最大のものから Nev 個の固有値を求める. = "SM": 絶対値最小のものから Nev 個の固有値を求める. = "BE": スペクトルの両端から半分づつ, 合計 Nev 個の固有値を求める. Nev が奇数の場合, 大きい側を1個多く求める. [in] Nev 求める固有値の数. (0 < Nev < N) [in] Ncv 行列 V の列数. (Nev < Ncv <= N) 各反復において生成されるランチョス基底ベクトルの数 (Ncv >= 2*Nev が推奨される). [out] D() 配列 D(LD - 1) (LD >= Nev) A*z = λ*B*z の固有値を近似するリッツ値を返す. 値は昇順に格納される. [out] Z() 配列 Z(Ldz - 1, LZ - 1) (Ldz >= N, LZ >= Nev) 固有値問題 A*z = λ*B*z のリッツ値に対応するリッツベクトル(近似固有ベクトル)を返す. Jobz = "N" の場合, Z() は参照されない. [out] Info リターンコード. = 0: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: 最大反復回数に達した. = 3: リスタート付きランチョス反復中にシフト値が適用できなかった. Ncv を Nev に比較して大きくすると改善するかもしれない (Ncv >= 2*Nev が推奨される). = 11: 初期残差ベクトルが 0 である. = 12: ランチョス分解に失敗した. = 13: LAPACK ルーチンで固有値の計算に失敗した. = 14: Dsaupd は十分な精度の固有値を得られなかった. = 15 〜 19: 内部処理エラー. [out] Nconv (省略可) 収束条件を満たしたリッツ値の数. [out] Niter (省略可) ランチョス反復回数. [in] Base (省略可) Ptr() および Ind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Ptr(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Base2 (省略可) Ptr2() および Ind2() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Ptr2(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Format (省略可) 行列の格納形式. (省略時 = 0) = 0: CSR 形式. = 1: CSC 形式. [in] Mode (省略可) 動作モード. (省略時 = 1) = 2: 通常インバースモード. = 3: シフト＆インバートモード. [in] Sigma (省略可) シフト (σ). Mode = 2 のときには参照されない. (省略時 = 0) [in] Maxiter (省略可) ランチョス反復の最大回数. (MaxIter > 0) (省略時 = 300) [in] Tol (省略可) 停止基準: リッツ値の最大相対誤差. (省略時 = 0) Tol <= 0 の場合, マシン精度とみなす. 使用例 対称行列 A および B について大きい方から2個の一般化固有値とその固有ベクトルを求める. ただし, ( 0.54 -0.90 -0.94 ) ( 1.18 0.54 -1.22 ) A = ( -0.90 0.70 1.04 ) B = ( 0.54 0.60 -0.71 ) ( -0.94 1.04 1.65 ) ( -1.22 -0.71 1.66 ) とする. Sub Ex_Dssgv() Const N = 3, Nnz = 6, Ncv = N, Nev = Ncv - 1 Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(Nnz - 1) As Double, Ib(N) As Long, Jb(Nnz - 1) As Long Dim D(Nev - 1) As Double, Z(N - 1, Nev - 1) As Double Dim Nconv As Long, Info As Long A(0) = 0.54: A(1) = -0.9: A(2) = 0.7: A(3) = -0.94: A(4) = 1.04: A(5) = 1.65 Ia(0) = 0: Ia(1) = 1: Ia(2) = 3: Ia(3) = 6 Ja(0) = 0: Ja(1) = 0: Ja(2) = 1: Ja(3) = 0: Ja(4) = 1: Ja(5) = 2 B(0) = 1.18: B(1) = 0.54: B(2) = 0.6: B(3) = -1.22: B(4) = -0.71: B(5) = 1.66 Ib(0) = 0: Ib(1) = 1: Ib(2) = 3: Ib(3) = 6 Jb(0) = 0: Jb(1) = 0: Jb(2) = 1: Jb(3) = 0: Jb(4) = 1: Jb(5) = 2 Call Dssgv("V", "L", "L", N, A(), Ia(), Ja(), B(), Ib(), Jb(), "LA", Nev, Ncv, D(), Z(), Info, Nconv) Debug.Print "D =", D(0), D(1) Debug.Print "Z =" Debug.Print Z(0, 0), Z(0, 1) Debug.Print Z(1, 0), Z(1, 1) Debug.Print Z(2, 0), Z(2, 1) Debug.Print "Nconv =" + Str(Nconv) + ", Info =" + Str(Info) End Sub 実行結果 D = 0.510423243055614 7.37370278804149 Z = 1.46384155786189 -9.13803098094185E-02 -0.442269445774234 1.71612038711754 1.35130854246605 0.947378600715302 Nconv = 2, Info = 0