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

ZCocr_s() Sub ZCocr_s ( N As  Long, Matvec As  LongPtr, Psolve As  LongPtr, ChkConv As  LongPtr, B() As  Complex, X() As  Complex, Optional Info As  Long, Optional Iter As  Long, Optional Res As  Double, Optional Mode As  Long = 0, Optional MaxIter As  Long = 500  ) COCR(Conjugate Orthogonal Conjugate Residual)法による連立一次方程式 Ax = b の解 (複素対称行列) (サブルーチン形式) 目的 前処理付き反復法(共役勾配(COCR)法)により複素対称行列を係数とする連立一次方程式 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] Psolve 前処理行列の適用, すなわち方程式 M*x = b の解をを求めるユーザーサブルーチンで, 次のように定義すること. ここで, M は前処理行列である. Sub Psolve(N As Long, B() As Complex, X() As Complex) M*x = b の解を求め 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: 正常終了. < 0: (-Info)番目の入力パラメータの誤り. = 1: (警告) 行列 A が正定値でない(計算は続行する). = 2: (警告) 前処理行列 M が正定値でない(計算は続行する). = 11: 最大反復回数を超えた. = 12: 行列 A が特異である(対角要素が0). [out] Iter (省略可) 収束時の反復回数. [out] Res (省略可) 最終的な残差ノルム norm(b - A*x) の値. [in] Mode (省略可) 引数 ChkCnv および Res において返す残差ノルムを選択することができる. (省略時 = 0) = 0: norm(b - A*x) を返す. ただし, 反復ごとに１回のMatvec演算が追加で必要になる. = 1: M^(-1)*norm(b - A*x) を返す. [in] MaxIter (省略可) 最大反復回数. (MaxIter > 0) (省略時 = 500) 使用例 連立一次方程式 Ax = B を解く. ただし, ( 0.31+0.77i 0.25+0.23i -0.81-0.83i ) A = ( 0.25+0.23i 0.26-0.26i -0.58-0.08i ) ( -0.81-0.83i -0.56-0.08i 2.09+0.6i ) ( 0.3941-1.2711i ) B = ( 0.0036-0.72i ) ( 0.3628+1.9977i ) とする. Const N = 3, Nnz = N * N, 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 SsrZusmv("L", N, Cmplx(1), A(), Ia(), Ja(), X(), Cmplx(0), Y()) End Sub Sub Psolve(N As Long, B() As Complex, X() As Complex) Dim I As Long For I = 0 To N - 1 X(I) = B(I) Next 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_ZCocr_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(0.31, 0.77): A(1) = Cmplx(0.25, 0.23): A(2) = Cmplx(0.26, -0.26): A(3) = Cmplx(-0.81, -0.83): A(4) = Cmplx(-0.56, -0.08): A(5) = Cmplx(2.09, 0.6) 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(0.3941, -1.2711): B(1) = Cmplx(0.0036, -0.72): B(2) = Cmplx(0.3628, 1.9977) Call ZCocr_s(N, AddressOf Matvec, AddressOf Psolve, AddressOf ChkConv, B(), X(), Info, Iter, Res) 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.820000000000007,-0.940000000000006) (0.739999999999996,0.2) (0.480000000000016,0.210000000000004) Iter = 3, Res = 5.1647781609597E-14, Info = 0