Initial

getDivGrad.jl

@@ -0,0 +1,37 @@
+using LinearOperators
+using MatrixDepot
+function getDivGrad(n1,n2,n3)
+	# D = getDivGrad(n1,n2,n3)
+	# builds 3D divergence operator
+	D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
+	D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
+	D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))
+
+	Div = [D1 D2 D3]
+	return Div*Div'
+end
+
+function ddx(n)
+# generate 1D finite difference on staggered grid
+	return d = spdiags(ones(n)*[-1 1],[0,1],n,n+1)
+end
+
+function spdiags(B,d,m,n)
+# A = spdiags(B,d,m,n)
+# creates a sparse matrix from its diagonals
+	d = d[:]
+	p = length(d)
+	len = zeros(p+1,1)
+	for k = 1:p
+		len[k+1] = round.(Int,len[k]+length(max(1,1-d[k]):min(m,n-d[k])))
+	end
+	a = zeros(round.(Int,len[p+1]),3)
+	for k = 1:p
+		# Append new d[k]-th diagonal to compact form
+		i = max(1,1-d[k]):min(m,n-d[k])
+		a[(round.(Int,len[k])+1):round.(Int,len[k+1]),:] = [i i+d[k] B[i+(m>=n)*d[k],k]]
+	end
+
+	A = sparse(round.(Int,a[:,1]),round.(Int,a[:,2]),a[:,3],m,n)
+	return A
+end

ssor.jl

@@ -0,0 +1,72 @@
+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 gidx = A.colptr[1] : A.colptr[n]
+				if A.rowval[gidx] == i
+					value -= A.nzval[gidx] * x[gidxi]
+				end
+			end
+			x[i]  = x[i] + value*OmInvD[i]
+		end
+
+		#第二次迭代
+		for i = n:-1:1
+			value = b[i]
+			for gidx = A.colptr[1] : A.colptr[n]
+				if A.rowval[gidx] == i
+					value -= A.nzval[gidx] * x[i]
+				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

test.jl

@@ -0,0 +1,6 @@
+include("getDivGrad.jl")
+include("ssor.jl")
+A = getDivGrad(32,32,32)
+rhs = randn(size(A,1))
+tol = 1e-9
+x1,flag1,relres1,iter1    = ssor(A,rhs,tol=tol,maxIter=1000,out=2)