Merge pull request #13 from ranjanan/moreopt

Add symmetric Gauss Seidel smoothing, improve smoother

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

@@ -16,8 +16,8 @@ abstract type CoarseSolver end
 struct Pinv <: CoarseSolver
 end
 
-MultiLevel{Ti,Tv}(l::Vector{Level{Ti,Tv}}, A::SparseMatrixCSC{Ti,Tv}, presmoother, postsmoother; coarse_solver = Pinv()) =
-    MultiLevel(l, A, coarse_solver, presmoother, postsmoother)
+MultiLevel{Ti,Tv}(l::Vector{Level{Ti,Tv}}, A::SparseMatrixCSC{Ti,Tv}, presmoother, postsmoother) =
+    MultiLevel(l, A, Pinv(), presmoother, postsmoother)
 Base.length(ml) = length(ml.levels) + 1
 
 function Base.show(io::IO, ml::MultiLevel)
@@ -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/preconditioner.jl

@@ -6,8 +6,8 @@ end
 
 aspreconditioner(ml::MultiLevel) = Preconditioner(ml)
 
-\(p::Preconditioner, b) = p * b
-*(p::Preconditioner, b) = solve(p.ml, b, 1, V(), 1e-12)
+\(p::Preconditioner, b) = solve(p.ml, b, 1, V(), 1e-12)
+*(p::Preconditioner, b) = p.ml.levels[1].A * x
 
 A_ldiv_B!(x, p::Preconditioner, b) = copy!(x, p \ b)
 A_mul_B!(b, p::Preconditioner, x) = A_mul_B!(b, p.ml.levels[1].A, x)

src/smoother.jl

@@ -4,13 +4,46 @@ struct SymmetricSweep <: Sweep
 end
 struct ForwardSweep <: Sweep
 end
+struct BackwardSweep <: Sweep
+end
 struct GaussSeidel{S} <: Smoother
     sweep::S
 end
-GaussSeidel(;sweep = ForwardSweep()) = GaussSeidel(sweep)
+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)
+
+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
 
-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, ::BackwardSweep, A, x, b) =
+    gs!(A, b, x, size(A,1), -1, 1)
 
-smoother(s::GaussSeidel, ::ForwardSweep, A, x, b) =
-                        gauss_seidel!(x, A, b, maxiter = 1)
+
+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 = start:step:stop
+        rsum = z
+        d = z
+        for j in nzrange(A, i)
+            row = A.rowval[j]
+            val = A.nzval[j]
+            if i == row
+                d = val
+            else
+                rsum += val * x[row]
+            end
+        end
+
+        if d != 0
+            x[i] = (b[i] - rsum) / d
+        end
+    end
+end

test/runtests.jl

@@ -100,23 +100,39 @@ end
 end
 
 @testset "Solver" begin
+fsmoother = GaussSeidel(ForwardSweep())
+
 A = poisson(1000)
 A = float.(A)
 ml = ruge_stuben(A)
 x = solve(ml, A * ones(1000))
-@test sum(abs2, x - ones(1000)) < 1e-10
+@test sum(abs2, x - ones(1000)) < 1e-8
+
+ml = ruge_stuben(A, presmoother = fsmoother,
+                     postsmoother = fsmoother)
+x = solve(ml, A * ones(1000))
+@test sum(abs2, x - ones(1000)) < 1e-8
+
 
 A = load("randlap.jld")["G"]
+
+ml = ruge_stuben(A, presmoother = fsmoother,
+                    postsmoother = fsmoother)
+x = solve(ml, A * ones(100))
+@test sum(abs2, x - zeros(100)) < 1e-8
+
 ml = ruge_stuben(A)
 x = solve(ml, A * ones(100))
-@test sum(abs2, x - zeros(100)) < 1e-10
+@test sum(abs2, x - zeros(100)) < 1e-6
 
 end
 
 @testset "Preconditioning" begin
 A = load("thing.jld")["G"]
 n = size(A, 1)
-ml = ruge_stuben(A)
+smoother = GaussSeidel(ForwardSweep())
+ml = ruge_stuben(A, presmoother = smoother,
+                    postsmoother = smoother)
 p = aspreconditioner(ml)
 b = zeros(n)
 b[1] = 1
@@ -163,5 +179,34 @@ diff = x - [ 0.82365077, -0.537589  , -0.30632349, -0.19370186, -0.14773294,
         0.04968258,  0.04968737,  0.05105749,  0.05009268,  0.04972329,
         0.04970173]
 @test sum(abs2, diff) < 1e-8
+
+# Symmetric GaussSeidel Smoothing
+
+ml = ruge_stuben(A)
+p = aspreconditioner(ml)
+
+x = cg(A, b, Pl = p, maxiter = 100_000, tol = 1e-6)
+diff = x - [0.823762, -0.537478, -0.306212, -0.19359, -0.147621, 0.685002,
+            -0.155389, -0.127703, -0.111867, 0.453735, -0.0856607, -0.0858715,
+            -0.0846678, 0.129962, 0.0281662, -0.0389642, -0.0593981, -0.0653311,
+            0.0545782, -0.0474255, -0.0519275, -0.0467483, -0.0448061, 0.056504,
+            0.0280386, -0.0227138, -0.0405172, -0.0431067, -0.0440778, 0.076042,
+            0.052232, 0.0447537, 0.05847, 0.0509098, 0.0484189, 0.0528356,
+            0.0503983, 0.0495933, 0.0497211, 0.0497731, 0.0497942, 0.049799,
+            0.0511691, 0.0502043, 0.0498349, 0.0498134]
+@test sum(abs2, diff) < 1e-8
+
+x = solve(ml, b, 1, V(), 1e-12)
+diff = x - [0.775725, -0.571202, -0.290989, -0.157001, -0.106981, 0.622652,
+            -0.122318, -0.0891874, -0.0709834, 0.392621, -0.055544, -0.0507485,
+            -0.0466376, 0.107175, 0.0267468, -0.0200843, -0.0282827, -0.0299929,
+            0.0420468, -0.0175585, -0.0181318, -0.0121591, -0.00902523, 0.0394795,
+            0.019981, -0.00270916, -0.0106855, -0.0093661, -0.00837619, 0.052532,
+            0.0301423, 0.0248904, 0.0333098, 0.0262179, 0.0246211, 0.026778,
+            0.0245746, 0.0238448, 0.0233892, 0.0231593, 0.0230526, 0.0229771,
+            0.0247913, 0.0238555, 0.0233681, 0.023096]
+@test sum(abs2, diff) < 1e-8
+
+
 end
 end