ZBicg ZBicg_r ZBicg_s ZCgs ZCgs_r ZCgs_s ZCocg ZCocg_r ZCocg_s ZCocr ZCocr_r ZCocr_s ZDiom ZDiom_r ZDiom_s ZDqgmres ZDqgmres_r ZDqgmres_s ZFgmres ZFgmres_r ZFgmres_s ZFom ZFom_r ZFom_s ZGcr ZGcr_r ZGcr_s ZGpbicg ZGpbicg_r ZGpbicg_s ZOrthomin ZOrthomin_r ZOrthomin_s ZQmr ZQmr_r ZQmr_s ZSor ZSor_r ZSor_s ZTfqmr ZTfqmr_r ZTfqmr_s

ZSor_s() Sub ZSor_s ( N As  Long, Matvec As  LongPtr, Matsol As  LongPtr, ChkConv As  LongPtr, B() As  Complex, X() As  Complex, Optional Info As  Long, Optional Iter As  Long, Optional Res As  Double, Optional MaxIter As  Long = 500  ) 逐次的過剰緩和(SOR)法による連立一次方程式 Ax = b の解 (複素行列) (サブルーチン形式) 目的 反復法(逐次的過剰緩和(SOR)法)により連立一次方程式 Ax = b の解を求める. 引数 [in] N 行列 A の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) [in] Matvec 行列 A とベクトル x の積をを求めるユーザーサブルーチンで, 次のように定義すること. Sub Matvec(N As Long, X() As Complex, Y() As Complex) A*x を計算し Y() に入れる. End Sub [in] Matsol 行列 A の対角および下三角部分を使って (D/ω + L)^(-1)*b を求めるユーザーサブルーチンで, 次のように定義すること. Sub Matsol(N As Long, B() As Complex, X() As Complex) x = (D/ω + L)^(-1)*b を求め (あるいは (D/ω + L)*x = b を解き) 解 x を X() に返す, End Sub [in] ChkConv 反復ごとに呼び出され収束判定を行うユーザーサブルーチンで, 次のように定義すること. ここで, X() は現在の近似解, Res は現在の残差ノルム norm(b - A*x), Iterは現在の反復回数である. 本ルーチンは中間結果を出力するために使用することもできる. Sub ChkConv(N As Long, X() As Complex, Res As Double, Iter As Long, IChk As Long) 収束であれば IChk = 1, そうでなければ IChk = 0 に設定する. End Sub [in] B() 配列 B(LB - 1) (LB >= N) 右辺ベクトル b. [in,out] X() 配列 X(LX - 1) (LX >= N) [in] 解の初期推定値. [out] 求められた近似解. [out] Info (省略可) = 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = 11: 最大反復回数を超えた. = 12: 行列が特異である(対角要素が0). [out] Iter (省略可) 収束時の反復回数. [out] Res (省略可) 最終的な残差ノルム norm(b - A*x) の値. [in] MaxIter (省略可) 最大反復回数. (MaxIter > 0) (省略時 = 500) 参考 (D/ω + L)^(-1)*b の計算は CscZussv または CsrZussv を用いて行うことができる. 使用例 連立一次方程式 Ax = B を解く. ただし, ( 4 0 1 0.7 0 ) ( 2i 4 0 1 0.7 ) A = ( 0 2i 4 0 1 ) ( 0 0 2i 4 0 ) ( 0 0 0 2i 4 ) ( 5.7 ) ( 5.7 + 2i ) B = ( 5 + 2i ) ( 4 + 2i ) ( 4 + 2i ) とする. Const N = 5, Nnz = 14, Omega = 1.05, Tol = 0.0000000001 Dim A(Nnz - 1) As Complex, Ia(N) As Long, Ja(Nnz - 1) As Long Sub Matvec(N As Long, X() As Complex, Y() As Complex) Call CsrZusmv("N", N, N, Cmplx(1), A(), Ia(), Ja(), X(), Cmplx(0), Y()) End Sub Sub Matsol(N As Long, B() As Complex, X() As Complex) Dim I As Long For I = 0 To N - 1 X(I) = B(I) Next Call CsrZussv("L", "N", "N", N, A(), Ia(), Ja(), X(), , , , Omega) End Sub Sub ChkConv(N As Long, X() As Complex, Res As Double, Iter As Long, IChk As Long) If Res < Tol Then IChk = 1 Else IChk = 0 End If End Sub Sub Ex_ZSor_s() Dim B(N - 1) As Complex, X(N - 1) As Complex Dim Iter As Long, Res As Double, Info As Long A(0) = Cmplx(4): A(1) = Cmplx(1): A(2) = Cmplx(0.7): A(3) = Cmplx(0, 2): A(4) = Cmplx(4): A(5) = Cmplx(1): A(6) = Cmplx(0.7): A(7) = Cmplx(0, 2): A(8) = Cmplx(4): A(9) = Cmplx(1): A(10) = Cmplx(0, 2): A(11) = Cmplx(4): A(12) = Cmplx(0, 2): A(13) = Cmplx(4) Ia(0) = 0: Ia(1) = 3: Ia(2) = 7: Ia(3) = 10: Ia(4) = 12: Ia(5) = 14 Ja(0) = 0: Ja(1) = 2: Ja(2) = 3: Ja(3) = 0: Ja(4) = 1: Ja(5) = 3: Ja(6) = 4: Ja(7) = 1: Ja(8) = 2: Ja(9) = 4: Ja(10) = 2: Ja(11) = 3: Ja(12) = 3: Ja(13) = 4 B(0) = Cmplx(5.7): B(1) = Cmplx(5.7, 2): B(2) = Cmplx(5, 2): B(3) = Cmplx(4, 2): B(4) = Cmplx(4, 2) Call ZSor_s(N, AddressOf Matvec, AddressOf Matsol, AddressOf ChkConv, B(), X(), Info, Iter, Res) Debug.Print "X =" Debug.Print "(" + Str(Creal(X(0))) + "," + Str(Cimag(X(0))) + ")" Debug.Print "(" + Str(Creal(X(1))) + "," + Str(Cimag(X(1))) + ")" Debug.Print "(" + Str(Creal(X(2))) + "," + Str(Cimag(X(2))) + ")" Debug.Print "(" + Str(Creal(X(3))) + "," + Str(Cimag(X(3))) + ")" Debug.Print "(" + Str(Creal(X(4))) + "," + Str(Cimag(X(4))) + ")" Debug.Print "Iter =" + Str(Iter) + ", Res =" + Str(Res) + ", Info =" + Str(Info) End Sub 実行結果 X = ( .999999999990298, 5.65196842116663E-12) ( 1.00000000000765, 6.21908808449608E-12) ( 1.00000000000282,-5.93301631217148E-12) ( .999999999997143,-3.11025281348357E-13) ( 1.00000000000028, 1.21430012714769E-12) Iter = 12, Res = 4.58544841000432E-11, Info = 0