|
|
◆ ZSsrSsorSolve()
| Sub ZSsrSsorSolve |
( |
Uplo As |
String, |
|
|
N As |
Long, |
|
|
Omega As |
Double, |
|
|
Val() As |
Complex, |
|
|
Rowptr() As |
Long, |
|
|
Colind() As |
Long, |
|
|
D() As |
Complex, |
|
|
B() As |
Complex, |
|
|
X() As |
Complex, |
|
|
Optional Info As |
Long, |
|
|
Optional Base As |
Long = -1 |
|
) |
| |
対称逐次的過剰緩和(SSOR)前処理 (複素対称行列) (CSR)
- 目的
- 連立一次方程式の疎で対称な係数行列 A に対する対称逐次的過剰緩和(SSOR)前処理を行う. すなわち, 連立一次方程式 M*x = b を解く. ここで, M は前処理行列である. A の上三角部分または下三角部分のみを使って計算を行う.
- 引数
-
| [in] | Uplo | = "U": 行列 A は上三角部分に格納されている.
= "L": 行列 A は下三角部分に格納されている. |
| [in] | N | 前処理行列の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る) |
| [in] | Omega | 緩和パラメータω. (0 < ω < 2) |
| [in] | Val() | 配列 Val(LVal - 1) (LVal >= Nnz)
行列 A の非ゼロ要素の値. (Nnz は非ゼロ要素数) |
| [in] | Rowptr() | 配列 Rowptr(LRowptr - 1) (LRowptr >= N + 1)
行列 A の行ポインタ. |
| [in] | Colind() | 配列 Colind(LColind - 1) (LColind >= Nnz)
行列 A の列インデクス. (Nnz は非ゼロ要素数) |
| [in] | D() | 配列 D(LD) (LD >= N)
CSR_SSOR()で求めた前処理行列 M の対角要素. |
| [in] | B() | 配列 B(LB - 1) (LB >= N)
右辺ベクトル b. |
| [out] | X() | 配列 X(LX - 1) (LX >= N)
解ベクトル x. |
| [out] | Info | (省略可)
= 0: 正常終了.
= i < 0: (-i)番目の入力パラメータの誤り.
= j > 0: 行列が特異である(j番目の対角要素が0). |
| [in] | Base | (省略可)
Rowptr() および Colind() のインデクス形式.
= 0: 0-ベース(C形式): 開始インデクス値が 0.
= 1: 1-ベース(Fortran形式): 開始インデクス値が 1.
(省略時: Rowptr(0) = 1 であれば 1, そうでなければ 0 とみなす) |
- 使用例
- 連立一次方程式 Ax = B を SSOR 前処理付き COCG 法で解く. ただし,
( 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 )
Sub Ex_ZCocg_Ssor_Csr()
Const N = 3, Nnz = 6, Omega = 0.9, 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(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) Dim D(N - 1) As Complex
Call ZCsxSsor(N, Omega, A(), Ia(), Ja(), D(), Info) If Info <> 0 Then Debug.Print "Ssor Info =" + Str(Info)
IRev = 0
Do
Call ZCocg_r(N, B(), X(), Info, XX(), YY(), IRev, Iter, Res) If IRev = 1 Then '- Matvec
ElseIf IRev = 3 Then '- Psolve
Call ZSsrSsorSolve("L", N, Omega, A(), Ia(), Ja(), D(), YY(), XX(), Info) If Info <> 0 Then Debug.Print "SsorSolve Info =" + Str(Info)
ElseIf IRev = 10 Then '- Check convergence
If Res < Tol Then IRev = 11
End If
Loop While IRev <> 0
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
Function Cmplx(R As Double, Optional I As Double=0) As Complex 複素数の作成 Function Cimag(A As Complex) As Double 複素数の虚数部 Function Creal(A As Complex) As Double 複素数の実数部 Sub SsrZusmv(Uplo As String, 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 (CSR) (複素対称行列) Sub ZCocg_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) COCG(Conjugate Orthogonal Conjugate Gradient)法による連立一次方程式 Ax = b の解 (複素対称行列) (リバースコミュニケーション版) Sub ZSsrSsorSolve(Uplo As String, N As Long, Omega As Double, Val() As Complex, Rowptr() As Long, Colind() As Long, D() As Complex, B() As Complex, X() As Complex, Optional Info As Long, Optional Base As Long=-1) 対称逐次的過剰緩和(SSOR)前処理 (複素対称行列) (CSR) Sub ZCsxSsor(N As Long, Omega As Double, Val() As Complex, Ptr() As Long, Ind() As Long, D() As Complex, Optional Info As Long, Optional Base As Long=-1) 対称逐次的過剰緩和(SSOR)前処理のための初期化 (複素行列) (CSC/CSR)
- 実行結果
X =
(-0.819999999999999,-0.939999999999999)
(0.74,0.2)
(0.48,0.21)
Iter = 3, Res = 7.87074880010101E-16, Info = 0
|
1.9.7