|
|
◆ Sor_r()
| Sub Sor_r |
( |
N As |
Long, |
|
|
B() As |
Double, |
|
|
X() As |
Double, |
|
|
Info As |
Long, |
|
|
XX() As |
Double, |
|
|
YY() As |
Double, |
|
|
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 演算. A*XX を求め 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 の計算は CscDussv または CsrDussv を用いて行うことができる.
- 使用例
- 連立一次方程式 Ax = B を解く. ただし,
( 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_Sor_r()
Const N = 3, Nnz = N * N, Omega = 0.4, 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, I 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
IRev = 0
Do
Call Sor_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 = 6 Then '- Matsol
For I = 0 To N - 1
YY(I) = XX(I)
Next
Call CsrDussv("L", "N", "N", N, A(), Ia(), Ja(), YY(), Omega:=Omega) 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
Sub CsrDussv(Uplo As String, Trans As String, Diag As String, N As Long, Val() As Double, Rowptr() As Long, Colind() As Long, X() As Double, Optional Info As Long, Optional Base As Long=-1, Optional IncX As Long=1, Optional Omega As Double=1) Ax = b または ATx = b の解 (三角行列) (CSR) Sub CsrDusmv(Trans As String, M As Long, N As Long, Alpha As Double, Val() As Double, Rowptr() As Long, Colind() As Long, X() As Double, Beta As Double, Y() As Double, 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 (CSR) Sub Sor_r(N As Long, B() As Double, X() As Double, Info As Long, XX() As Double, YY() As Double, IRev As Long, Optional Iter As Long, Optional Res As Double, Optional MaxIter As Long=500) 逐次的過剰緩和(SOR)法による連立一次方程式 Ax = b の解 (リバースコミュニケーション版)
- 実行結果
X = 0.859999999965658 0.640000000003731 0.509999999994652
Iter = 57, Res = 2.83730635879748E-11, Info = 0
|
1.9.7