◆ ZCsxDsSolve()

Sub ZCsxDsSolve	(	N As	Long,
		D() As	Complex,
		B() As	Complex,
		X() As	Complex,
		Optional Info As	Long
	)

対角スケーリング前処理 (複素行列) (CSC/CSR)

目的: 連立一次方程式の疎な係数行列 A に対する対角スケーリング前処理を行う. すなわち, 連立一次方程式 M*x = b を解く. ここで, M (= A の対角行列) は前処理行列である.

引数

[in]	N	前処理行列の次数. (N >= 0) (N = 0 の場合, 処理を行わずに戻る)
[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).

使用例: 連立一次方程式 Ax = B を DS 前処理付き FGMRES 法で解く. ただし,
( 0.78+0.16i -0.9-1.46i 0.48-1.08i )

A = ( 0.73+0.63i 1.58-1.24 -0.41-0.91i )

( 0.23-1.37i 0.79+0.64i -0.73-1.5i )

( 0.2126-0.2904i )

B = ( -0.3028+0.3346i )

( -1.2905-1.0346i )

とする.
Sub Ex_ZFgmres_Ds()

Const N = 3, Nnz = N * N, 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, I As Long

A(0) = Cmplx(0.78, 0.16): A(1) = Cmplx(-0.9, -1.46): A(2) = Cmplx(0.48, -1.08): A(3) = Cmplx(0.73, 0.63): A(4) = Cmplx(1.58, -1.24): A(5) = Cmplx(-0.41, -0.91): A(6) = Cmplx(0.23, -1.37): A(7) = Cmplx(0.79, 0.64): A(8) = Cmplx(-0.73, -1.5)

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) = Cmplx(0.2126, -0.2904): B(1) = Cmplx(-0.3028, 0.3346): B(2) = Cmplx(-1.2905, -1.0346)

Dim D(N - 1) As Complex

Call ZCsxDs(N, A(), Ia(), Ja(), D(), Info)

If Info <> 0 Then Debug.Print "Ds Info =" + Str(Info)

IRev = 0

Do

Call ZFgmres_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 = 3 Then '- Psolve

Call ZCsxDsSolve(N, D(), YY(), XX())

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

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
複素数の実数部

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)

ZFgmres_r
Sub ZFgmres_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 M As Long=0, Optional MaxIter As Long=500)
最小残差(FGMRES)法による連立一次方程式 Ax = b の解 (複素行列) (リバースコミュニケーション版)

ZCsxDsSolve
Sub ZCsxDsSolve(N As Long, D() As Complex, B() As Complex, X() As Complex, Optional Info As Long)
対角スケーリング前処理 (複素行列) (CSC/CSR)

ZCsxDs
Sub ZCsxDs(N As Long, Val() As Complex, Ptr() As Long, Ind() As Long, D() As Complex, Optional Info As Long, Optional Base As Long=-1)
対角スケーリング前処理のための初期化 (複素行列) (CSC/CSR)

実行結果: X =

(0.59,-0.28)

(-0.2,-0.04)

(0.24,-0.49)

Iter = 3, Res = 3.24824839330284E-16, Info = 0