MatrixSolver/ssor.jl

export ssor

"""
x,flag,err,iter,resvec = ssor(A::SparseMatrixCSC, b::Vector; omega::Real=2/3, tol::Real=1e-2, maxIter::Int=1,out::Int=0)

Symmetric successive over-relaxation applied to the linear system A*x = b.

!! Note that upper/lower triangular matrix are never build to save memory !!

Input:

	A       - sparse matrix CSC
	b       - right hand side vector
	omega   - relaxation parameter (omega=1.0: Gauss Seidel)
	tol     - error tolerance (default=1e-2)
	maxIter - maximum number of iterations (default=1)
	out     - flag for output (-1 : no output, 0 : only errors, 1 : final status, 2: error at each iteration)

Output:

	x       - solution
	flag    - exit flag (  0 : desired tolerance achieved,
	                      -1 : maxIter reached without converging)
	err     - error, i.e., norm(A*x-b)/norm(b)
	iter    - number of iterations
	resvec  - error at each iteration
"""
function ssor(A::SparseMatrixCSC,b::Vector;tol::Real=1e-2,maxIter::Int=1,omega::Real=1.0,out::Int=0,storeInterm::Bool=false;)
	n = length(b)
	OmInvD = omega./diag(A) #D矩阵 并乘以omega
	x = zeros(eltype(A),n)
	r = copy(b)
	rnorm0 = norm(b)
	res = 0
	value = 0.0
	flag   = -1
	iter   = 0
	dx = 0
	for iter = 1:maxIter
		println(iter)
		#第一次迭代 按照列来
		for j = 1:n
			dx =  r[j] * OmInvD[j]
			for gidx = A.colptr[j] : A.colptr[j+1]-1
				r[A.rowval[gidx]] -= A.nzval[gidx] * dx
			end
			x[j] += dx
		end
		for j = n:-1:1
			dx =  r[j] * OmInvD[j]
			for gidx = A.colptr[j] : A.colptr[j+1]-1
				r[A.rowval[gidx]] -= A.nzval[gidx] * dx
			end
			x[j] += dx
		end
		res = norm(A*x-b) / rnorm0

		println(res)

		if res < tol
			flag = 0; break
		end
	end

	return x,flag,iter

end