◆ 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`
	)

Solution of linear system Ax = b using Successive over-relaxation (SOR) method (Complex matrices) (Reverse communication version)

Purpose: This routine solves the linear system Ax = b using successive over-relaxation (SOR) iterative method.

Parameters

[in]	N	Dimension of the matrix. (N >= 0) (if N = 0, returns without computation)
[in]	B()	Array B(LB - 1) (LB >= N) Right hand side vector b.
[in,out]	X()	Array X(LX - 1) (LX >= N) [in] Initial guess of solution. [out] Obtained approximate solution.
[out]	Info	= 0: Successful exit = i < 0: The (-i)-th argument is invalid. = 11: Maximum number of iterations exceeded.
[in,out]	XX()	Array XX(LXX - 1) (LXX >= N) Array for Matvec and Matsol operations.
[in,out]	YY()	Array YY(LYY - 1) (LYY >= N) Array for Matvec and Matsol operations.
[in,out]	IRev	Control variable for reverse communication. [in] Before first call, IRev should be initialized to zero. On succeeding calls, IRev should not be altered (except if converged). [out] If IRev is not zero, complete the following process and call this routine again. = 0: Computation finished. See return code in Info. = 1: Matvec operation. User should set AXX in YY. Do not alter any other variables. = 6: Matsol operation. User must compute (D/ω + L)^(-1)XX and set in YY. Do not alter any other variables. ω is the relaxation parameter of SOR method (0 < ω < 2). L and D are the lower triangular part and the diagonal part of the matrix A, respectively. = 10: To be returned for the convergence test on every iteration . Set IRev = 11 if converged. Do not alter IRev otherwise. The latest values in X(), Iter and Res can be used to decide the convergence. Further, these values may be used to output the intermediate results.
[out]	Iter	(Optional) Actual number of iterations performed for convergence.
[out]	Res	(Optional) Final residual norm value norm(b - A*x).
[in]	MaxIter	(Optional) Maximum number of iterations. (MaxIter > 0) (default = 500)

Note: To calculate (D/ω + L)^(-1)*b, CscZussv or CsrZussv can be used.

Example Program: Solve the system of linear equations Ax = B, where
( 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
Building complex number

Cimag
Function Cimag(A As Complex) As Double
Imaginary part of complex number

Creal
Function Creal(A As Complex) As Double
Real part of complex number

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)
Solution of Ax = b, ATx = b or AHx = b (Complex triangular matrices) (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 or y <- αAHx + βy (Complex matrices) (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)
Solution of linear system Ax = b using Successive over-relaxation (SOR) method (Complex matrices) (Re...

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 (complex vector) (BLAS 1)

Example Results: 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