Add symmetric smoothing to improve CG convergence

Author: ranjanan

src/AMG.jl

@@ -12,6 +12,7 @@ include("gallery.jl")
 export poisson
 
 include("smoother.jl")
+export GaussSeidel, SymmetricSweep, ForwardSweep, BackwardSweep
 
 include("multilevel.jl")
 export solve

src/classical.jl

@@ -20,9 +20,9 @@ function ruge_stuben{Ti,Tv}(A::SparseMatrixCSC{Ti,Tv};
 
     levels = Vector{Level{Ti,Tv}}()
 
-    while length(levels) < max_levels
+    while length(levels) < max_levels #&& size(A, 1) < max_coarse
         A = extend_heirarchy!(levels, strength, CF, A)
-        if size(A, 1) < max_coarse
+        if size(A, 1) <= max_coarse
             break
         end
     end

src/multilevel.jl

@@ -61,19 +61,20 @@ end
 function solve{T}(ml::MultiLevel, b::Vector{T},
                                     maxiter = 100,
                                     cycle = V(),
-                                    tol = 1e-5)
+                                    tol = 1e-5;
+                                    verbose = false)
     x = zeros(T, size(b))
     residuals = Vector{T}()
     A = ml.levels[1].A
     normb = norm(b)
     if normb != 0
         tol *= normb
     end
-    push!(residuals, T(norm(b - A*x)))
+    push!(residuals, normb)
 
     lvl = 1
     while length(residuals) <= maxiter && residuals[end] > tol
-        if length(ml.levels) == 1
+        if length(ml) == 1
             x = coarse_solver(ml.coarse_solver, A, b)
         else
             x = __solve(cycle, ml, x, b, lvl)

src/smoother.jl

@@ -4,36 +4,46 @@ struct SymmetricSweep <: Sweep
 end
 struct ForwardSweep <: Sweep
 end
+struct BackwardSweep <: Sweep
+end
 struct GaussSeidel{S} <: Smoother
     sweep::S
 end
-GaussSeidel() = GaussSeidel(ForwardSweep())
+GaussSeidel() = GaussSeidel(SymmetricSweep())
+
+presmoother!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
+postsmoother!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
 
-presmoother!(s, A, x, b) = smoother(s, s.sweep, A, x, b)
-postsmoother!(s, A, x, b) = smoother(s, s.sweep, A, x, b)
+smoother!(s::GaussSeidel, ::ForwardSweep, A, x, b) =
+                    gs!(A, b, x, 1, 1, size(A, 1))
+                    # gauss_seidel!(x, A, b, maxiter = 1)
+function smoother!(s::GaussSeidel, ::SymmetricSweep, A, x, b)
+    smoother!(s, ForwardSweep(), A, x, b)
+    smoother!(s, BackwardSweep(), A, x, b)
+end
 
-smoother(s::GaussSeidel, ::ForwardSweep, A, x, b) =
-                    gs!(A, b, x)
-                        #gauss_seidel!(x, A, b, maxiter = 1)
+smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b) =
+    gs!(A, b, x, size(A,1), -1, 1)
 
-function gs!{T,Ti}(A::SparseMatrixCSC{T,Ti}, b::Vector{T}, x::Vector{T})
+
+function gs!{T,Ti}(A::SparseMatrixCSC{T,Ti}, b::Vector{T}, x::Vector{T}, start, step, stop)
     n = size(A, 1)
     z = zero(eltype(A))
-    for i = 1:n
-        # rsum = calc_weighted_sum(A, x)
+    for i = start:step:stop
         rsum = z
-        diag = z
+        d = z
         for j in nzrange(A, i)
             row = A.rowval[j]
             val = A.nzval[j]
             if i == row
-                diag = val
+                d = val
             else
                 rsum += val * x[row]
             end
         end
-        if diag != 0
-            x[i] = (b[i] - rsum) / diag
+
+        if d != 0
+            x[i] = (b[i] - rsum) / d
         end
     end
 end