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 = b - A * x
	rnorm0 = norm(r)
	value = 0.0
	flag   = -1
	iter   = 0
	for iter = 1:maxIter
		println(iter)
		#第一次迭代
		for i = 1:n
			value = b[i]
			for j = 1:n
				for gidx = A.colptr[j] : A.colptr[j+1]-1
					if A.rowval[gidx] == i
						value -= A.nzval[gidx] * x[j]
					end
				end
			end
			x[i]  = x[i] + value*OmInvD[i]
		end

		#第二次迭代
		for i = n:-1:1
			value = b[i]
			for j = 1:n
				for gidx = A.colptr[j] : A.colptr[j+1]-1
					if A.rowval[gidx] == i
						value -= A.nzval[gidx] * x[j]
					end
				end
			end
			x[i]  = x[i] + value*OmInvD[i]
		end
		r = b - A*x
		println(norm(r) / rnorm0 )

		if norm(r) / rnorm0 < tol
			flag = 0; break
		end
	end

	println(r)
	println(x)
	return x,flag, norm(r) / rnorm0,iter

end