◆ ZSor_r()

Sub ZSor_r	(	N As	Long,
		B() As	Complex,
		X() As	Complex,
		Info As	Long,
		XX() As	Complex,
		YY() As	Complex,
		IRev 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]	B()	配列 B(LB - 1) (LB >= N) 右辺ベクトル b.
[in,out]	X()	配列 X(LX - 1) (LX >= N) [in] 解の初期推定値. [out] 求められた近似解.
[out]	Info	= 0: 正常終了. = i < 0: (-i)番目の入力パラメータの誤り. = 11: 最大反復回数を超えた.
[in,out]	XX()	配列 XX(LXX - 1) (LXX >= N) Matvec および Matsol 演算のための配列.
[in,out]	YY()	配列 YY(LYY - 1) (LYY >= N) Matvec および Matsol 演算のための配列.
[in,out]	IRev	リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に IRev = 0 に設定しておくこと. それ以降の呼び出し時にはIRevの値を変更してはならない(収束時を除く). [out] 0 以外のときには下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかは Info をチェックすること. = 1: Matvec 演算. AXX を求め YY に設定すること. 他の変数を変更してはならない. = 6: Matsol 演算. (D/ω + L)^(-1)XX を求め YY に設定すること. 他の変数を変更してはならない. ただし, ωは SOR 法の緩和パラメータ (0 < ω < 2), L は A の下三角部分, D は A の対角部分を表す. = 10: 収束判定を行うために反復ごとに戻る. 収束時には IRev = 11 として, それ以外のときは IRev を変更せずに再度呼び出すこと. X(), Iter および Res にはその反復時の最新の値が入っているので収束判定のために使用してよい. また, 中間結果出力のために使うこともできる.
[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 )

とする.
Sub Ex_ZSor_r()

Const N = 5, Nnz = 14, Omega = 1.05, 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(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)

IRev = 0

Do

Call ZSor_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res)

If IRev = 1 Then '- Matvec

Call CsrZusmv("N", N, N, Cmplx(1), A(), Ia(), Ja(), XX(), Cmplx(0), YY())

ElseIf IRev = 6 Then '- Matsol

Call Zcopy(N, XX(0), YY(0))

Call CsrZussv("L", "N", "N", N, A(), Ia(), Ja(), YY(), , , , Omega)

ElseIf IRev = 10 Then '- Check convergence

If Res < Tol Then IRev = 11

End If

Loop While IRev <> 0

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

Cmplx
Function Cmplx(R As Double, Optional I As Double=0) As Complex
複素数の作成

Cimag
Function Cimag(A As Complex) As Double
複素数の虚数部

Creal
Function Creal(A As Complex) As Double
複素数の実数部

CsrZussv
Sub CsrZussv(Uplo As String, Trans As String, Diag As String, N As Long, Val() As Complex, Rowptr() As Long, Colind() As Long, X() As Complex, Optional Info As Long, Optional Base As Long=-1, Optional IncX As Long=1, Optional Omega As Double=1)
Ax = b, ATx = b または AHx = b の解 (複素三角行列) (CSR)

CsrZusmv
Sub CsrZusmv(Trans As String, M As Long, N As Long, Alpha As Complex, Val() As Complex, Rowptr() As Long, Colind() As Long, X() As Complex, Beta As Complex, Y() As Complex, Optional Info As Long, Optional Base As Long=-1, Optional IncX As Long=1, Optional IncY As Long=1)
y <- αAx + βy, y <- αATx + βy または y <- αAHx + βy (複素行列) (CSR)

ZSor_r
Sub ZSor_r(N As Long, B() As Complex, X() As Complex, Info As Long, XX() As Complex, YY() As Complex, IRev As Long, Optional Iter As Long, Optional Res As Double, Optional MaxIter As Long=500)
逐次的過剰緩和(SOR)法による連立一次方程式 Ax = b の解 (複素行列) (リバースコミュニケーション版)

Zcopy
Sub Zcopy(N As Long, ZX_I As Complex, ZY_I As Complex, Optional IncX As Long=1, Optional IncY As Long=1)
y <- x (複素ベクトル) (BLAS 1)

実行結果: 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