ZHscIc0 ZHscIcSolve ZHscSsorSolve ZHsrIc0 ZHsrIcSolve ZHsrSsorSolve

ZHsrIc0() Sub ZHsrIc0 ( N As  Long, Val() As  Complex, Rowptr() As  Long, Colind() As  Long, Val2() As  Complex, Idiag() As  Long, Optional Info As  Long, Optional ByVal Uplo As  String = "L", Optional Base As  Long = -1  ) 不完全コレスキー分解(フィルインなし)(IC0)前処理のための初期化 (正定値エルミート行列) (CSR) 目的 連立一次方程式の正定値エルミートで疎な係数行列の不完全コレスキー分解(フィルインなし)を求める. A = L * D * L^H + R または A = U^H * D * U + R ここで, L は下三角行列, U は上三角行列, D は対角行列である. 係数行列 A の要素が 0 であった場所がコレスキー分解後に 0 でなくなることをフィルインが生じるというが, それを避けるために上式の L または U ではそのような場所を強制的に 0 にすることにする. R はこの操作による完全なコレスキー分解との差分を表すが, R が大きくはないものとすると, この不完全分解を使って連立一次方程式を解いてもよい近似解になっていることが期待される. これを前処理行列 M として使用する. M = L * D * L^H または M = U^H * D * U 本ルーチンは Val2() に L または U および D を出力する. また, Idiag() に対角要素のインデックスを出力する. Val2() および Idiag() を ZHsrIcSolve() が使用する. 引数 [in] N 行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) (Nnz は行列 A の非ゼロ要素数) 行列 A の非ゼロ要素の値. (対角要素と下三角要素のみ参照する) [in] Rowptr() 配列 Rowptr(LRowptr - 1) (LRowptr >= N + 1) 行列 A の行ポインタ. [in] Colind() 配列 Colind(LColind - 1) (LColind >= Nnz) 行列 A の列インデクス. [out] Val2() 配列 Val2(LVal2 - 1) (LVal2 >= Nnz) 前処理行列 M の非ゼロ要素の値. (A の対角要素と下三角要素と同じ場所に書き込まれる) [out] Idiag() 配列 Idiag(LIdiag) (LIdiag >= N) 対角要素のインデックス. [out] Info (省略可) = 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = j > 0: 行列が特異である(j番目の対角要素が0). [in] Uplo (省略可) 入力行列 A が上三角部分あるいは下三角部分のどちらを格納しているか指定する. (省略時 = "L") = "U": 上三角部分 (U^H*D*U と分解する). = "L": 下三角部分 (L*D*L^H と分解する). [in] Base (省略可) Rowptr() および Colind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Rowptr(0) = 1 であれば 1, そうでなければ 0 とみなす) 使用例 連立一次方程式 Ax = B を IC(0) 前処理付き CG 法で解く. ただし, ( 1.4 -1.5+0.46i 0.16+0.23i ) A = ( -1.5-0.46i 1.44 -0.12+0.04i ) ( 0.16-0.23i -0.12-0.04i 0.05 ) ( -2.3215-1.1316i ) B = ( 1.7972+2.0692i ) ( -0.4042-0.0049i ) とする. Sub Ex_ZCg_Ic0_Hsr() Const N = 3, Nnz = 6, Tol = 0.0000000001 '1.0e-10 Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(N - 1) As Complex, X(N - 1) As Complex Dim XX(N - 1) As Complex, YY(N - 1) As Complex Dim Iter As Long, Res As Double, IRev As Long, Info As Long A(0) = Cmplx(1.4): A(1) = Cmplx(-1.5, -0.46): A(2) = Cmplx(1.44): A(3) = Cmplx(0.16, -0.23): A(4) = Cmplx(-0.12, -0.04): A(5) = Cmplx(0.05) 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) = Cmplx(-2.3215, -1.1316): B(1) = Cmplx(1.7972, 2.0692): B(2) = Cmplx(-0.4042, -0.0049) Dim M(Nnz - 1) As Complex, Id(N - 1) As Long Call ZHsrIc0(N, A(), Ia(), Ja(), M(), Id(), Info) If Info <> 0 Then Debug.Print "Ic0 Info =" + Str(Info) IRev = 0 Do Call ZCg_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res) If IRev = 1 Then '- Matvec Call HsrZusmv("L", N, Cmplx(1), A(), Ia(), Ja(), XX(), Cmplx(0), YY()) ElseIf IRev = 3 Then '- Psolve Call ZHsrIcSolve(N, M(), Ia(), Ja(), Id(), YY(), XX(), Info) If Info <> 0 Then Debug.Print "Ic0Solve Info =" + Str(Info) ElseIf IRev = 10 Then '- Check convergence If Res < Tol Then IRev = 11 End If Loop While IRev <> 0 Debug.Print "X =" Debug.Print "(" + CStr(Creal(X(0))) + "," + CStr(Cimag(X(0))) + ")" Debug.Print "(" + CStr(Creal(X(1))) + "," + CStr(Cimag(X(1))) + ")" Debug.Print "(" + CStr(Creal(X(2))) + "," + CStr(Cimag(X(2))) + ")" Debug.Print "Iter =" + Str(Iter) + ", Res =" + Str(Res) + ", Info =" + Str(Info) End Sub 実行結果 X = (-0.82,-0.939999999999999) (0.74,0.2) (0.479999999999999,0.209999999999999) Iter = 1, Res = 7.37282538234167E-16, Info = 0