CscIlu CscIlu0 CscIluSolve CscSsorSolve CsrIlu CsrIlu0 CsrIluSolve CsrSsorSolve CsxDs CsxDsSolve CsxSsor

CsrIlu() Sub CsrIlu ( N As  Long, Val() As  Double, Rowptr() As  Long, Colind() As  Long, P As  Long, Val2() As  Double, Rowptr2() As  Long, Colind2() As  Long, D() As  Double, Optional Info As  Long, Optional Md As  Long = 0, Optional Base As  Long = -1, Optional Base2 As  Long = 0  ) 不完全LU分解(レベル指定)(ILU(p))前処理のための初期化 (CSR) 目的 連立一次方程式の疎な係数行列 A の不完全LU分解を求める. A = L * U + R ここで, Rは完全なLU分解との差分であるが, Rが小さいものとしてこの分解を使って連立一次方程式を解くことにすると次の前処理行列が得られる. M = L * U フィルインを抑えるために, 分解により生じる新たな非ゼロ要素のうち指定された値 p よりも大きなレベル値のものを廃棄する. これを ILU(p) と呼ぶ. レベル値の初期値を次のように設定する. lev(a(i,j)) = 0 (a(i,j) != 0 または i = j の場合), = ∞ (a(i,j) = 0 の場合) LU分解においては消去の過程で a(i,j) = a(i,j) - a(i,k)*a(k,j) という演算を行うが, このとき, a(i,j)のレベル値を次のように更新する. lev(a(i,j)) = min{ lev(a(i,j)), lev(a(i,k)) + lev(a(k,j)) + 1 } このようにすると, 分解前の A の非ゼロ要素の場所は最後までレベル0のままになる. 分解後にレベル値が p よりも大きい要素は廃棄される. Val2(), Rowptr2() および Colind2() に下三角行列 L と上三角行列 U を出力する. また, D() に U の対角要素をコピーする. Val2(), Rowptr2(), Colind2() および D() を CsrIluSolve() が使用する. 引数 [in] N 行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Val() 配列 Val(LVal - 1) (LVal >= Nnz) (Nnz は非ゼロ要素数) 行列 A の非ゼロ要素の値. [in] Rowptr() 配列 Rowptr(LRowptr - 1) (LRowptr >= N + 1) 行列 A の行ポインタ. [in] Colind() 配列 Colind(LColind - 1) (LColind >= Nnz) 行列 A の列インデクス. [in] P ILU(p) アルゴリズムのレベル p の値. (P >= 0) 注: P = 0 とすることができるが, CsrIlu0()を使用した方が高速である. [out] Val2() 配列 Val2(LVal2 - 1) (LVal2 > 0) 下三角行列 L と上三角行列 U (行列 L*U) の非ゼロ要素の値. Nnz2 = min(LVal2, LColind2) とすると, Nnz2 は行列 L*U の非ゼロ要素数の最大値である. 分解中に非ゼロ要素数がこれを超えるとエラー番号 -6 で終了する. [out] Rowptr2() 配列 Rowptr2(LRowptr2 - 1) (LRowptr2 >= N + 1) 行列 L*U の行ポインタ. [out] Colind2() 配列 Colind2(LColind2 - 1) (LColind2 > 0) 行列 L*U の列インデクス. [out] D() 配列 D(LD) (LD >= N) 上三角行列 U の対角要素. [out] Info (省略可) = 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = j > 0: 行列が特異である(j番目の対角要素が0). [in] Md (省略可) 廃棄した要素の分を対角要素に修正を加えるかどうか指定する. このような修正を行う方法を修正ILU分解(Modified ILU (MILU)分解)という. (省略時 = 0) = 0: ILU (修正なし). = 1: MILU (対角要素の修正あり). [in] Base (省略可) Rowptr() および Colind() のインデクス形式. = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. (省略時: Rowptr(0) = 1 であれば 1, そうでなければ 0 とみなす) [in] Base2 (省略可) Rowptr2() および Colind2() のインデクス形式. (省略時 = 0) = 0: 0-ベース(C形式): 開始インデクス値が 0. = 1: 1-ベース(Fortran形式): 開始インデクス値が 1. 使用例 連立一次方程式 Ax = B を ILU(p) 前処理付き FGMRES 法で解く. ただし, ( 0.2 -0.11 -0.93 ) ( -0.3727 ) A = ( -0.32 0.81 0.37 ), B = ( 0.4319 ) ( -0.8 -0.92 -0.29 ) ( -1.4247 ) とする. Sub Ex_Fgmres_Ilu_Csr() Const N = 3, Nnz = N * N, NnzM = Nnz, Tol = 0.0000000001 '1.0e-10 Dim A(Nnz - 1) As Double, Ia(N) As Long, Ja(Nnz - 1) As Long Dim B(N - 1) As Double, X(N - 1) As Double Dim XX(N - 1) As Double, YY(N - 1) As Double Dim Iter As Long, Res As Double, IRev As Long, Info As Long A(0) = 0.2: A(1) = -0.11: A(2) = -0.93: A(3) = -0.32: A(4) = 0.81: A(5) = 0.37: A(6) = -0.8: A(7) = -0.92: A(8) = -0.29 Ia(0) = 0: Ia(1) = 3: Ia(2) = 6: Ia(3) = 9 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 B(0) = -0.3727: B(1) = 0.4319: B(2) = -1.4247 Dim M(NnzM - 1) As Double, Im(N) As Long, Jm(NnzM - 1) As Long Dim P As Long, D(N - 1) As Double P = 3 Call CsrIlu(N, A(), Ia(), Ja(), P, M(), Im(), Jm(), D(), Info) If Info <> 0 Then Debug.Print "Ilu Info =" + Str(Info) IRev = 0 Do Call Fgmres_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res) If IRev = 1 Then '- Matvec Call CsrDusmv("N", N, N, 1, A(), Ia(), Ja(), XX(), 0, YY()) ElseIf IRev = 3 Then '- Psolve Call CsrIluSolve("N", N, M(), Im(), Jm(), D(), YY(), XX(), Info) If Info <> 0 Then Debug.Print "IluSolve Info =" + Str(Info) ElseIf IRev = 10 Then '- Check convergence If Res < Tol Then IRev = 11 End If Loop While IRev <> 0 Debug.Print "X =", X(0), X(1), X(2) Debug.Print "Iter = " + CStr(Iter) + ", Res = " + CStr(Res) + ", Info = " + CStr(Info) End Sub 実行結果 X = 0.86 0.64 0.51 Iter = 1, Res = 2.00597268669339E-16, Info = 0