Add a number of gauss seidel smoothing tests

Author: ranjanan

src/aggregation.jl

@@ -64,7 +64,7 @@ function extend_hierarchy!(levels, strength, aggregate, smooth,
     # relax!(improve_candidates, A, B, b)
     T, B = fit_candidates(AggOp, B)
 
-    P = smooth_prolongator(smooth, A, T, S', B)
+    P = smooth_prolongator(smooth, A, T, S, B)
     R = construct_R(symmetry, P)
     push!(levels, Level(A, P, R))
 

src/multilevel.jl

@@ -88,6 +88,7 @@ function solve{T}(ml::MultiLevel, b::Vector{T},
         push!(residuals, T(norm(b - A * x)))
     end
 
+    @show residuals
     if log
         return x, residuals
     else

src/smoother.jl

@@ -18,8 +18,11 @@ presmoother!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
 postsmoother!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
 relax!(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))
+function smoother!(s::GaussSeidel, ::ForwardSweep, A, x, b)
+    for i in 1:s.iter
+        gs!(A, b, x, 1, 1, size(A, 1))
+    end
+end
 
 function smoother!(s::GaussSeidel, ::SymmetricSweep, A, x, b)
     for i in 1:s.iter
@@ -28,8 +31,11 @@ function smoother!(s::GaussSeidel, ::SymmetricSweep, A, x, b)
     end
 end
 
-smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b) =
-    gs!(A, b, x, size(A,1), -1, 1)
+function smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b)
+    for i in 1:s.iter
+        gs!(A, b, x, size(A,1), -1, 1)
+    end
+end
 
 
 function gs!(A, b, x, start, step, stop)

test/runtests.jl

@@ -258,4 +258,7 @@ end
 end
 end
 
+@testset "Gauss Seidel" begin
+        test_gauss_seidel()
+end
 end

test/sa_tests.jl

@@ -227,4 +227,79 @@ function test_approximate_spectral_radius()
 
     end
 
+end
+
+# Test Gauss Seidel 
+import AMG: gs!, relax!
+function test_gauss_seidel()
+    
+    N = 1
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                    (0, -1, 1), N, N) 
+    x = eltype(A).(collect(0:N-1))
+    b = zeros(N)
+    s = GaussSeidel(ForwardSweep())
+    relax!(s, A, x, b)
+    @test sum(abs2, x - zeros(N)) < 1e-8
+
+    N = 3 
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                        (0, -1, 1), N, N) 
+    x = eltype(A).(collect(0:N-1))
+    b = zeros(N)
+    s = GaussSeidel(ForwardSweep())
+    relax!(s, A, x, b)
+    @test sum(abs2, x - [1.0/2.0, 5.0/4.0, 5.0/8.0]) < 1e-8
+
+    N = 1
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                    (0, -1, 1), N, N) 
+    x = eltype(A).(collect(0:N-1))
+    b = zeros(N)
+    s = GaussSeidel(BackwardSweep())
+    relax!(s, A, x, b)
+    @test sum(abs2, x - zeros(N)) < 1e-8
+
+    N = 3 
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                        (0, -1, 1), N, N) 
+    x = eltype(A).(collect(0:N-1))
+    b = zeros(N)
+    s = GaussSeidel(BackwardSweep())
+    relax!(s, A, x, b)
+    @test sum(abs2, x - [1.0/8.0, 1.0/4.0, 1.0/2.0]) < 1e-8
+
+    N = 1
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                    (0, -1, 1), N, N) 
+    x = eltype(A).(collect(0:N-1))
+    b = eltype(A).([10.])
+    s = GaussSeidel(ForwardSweep())
+    relax!(s, A, x, b)
+    @test sum(abs2, x - [5.]) < 1e-8
+
+    N = 3 
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                        (0, -1, 1), N, N) 
+    x = eltype(A).(collect(0:N-1))
+    b = eltype(A).([10., 20., 30.])
+    s = GaussSeidel(ForwardSweep())
+    relax!(s, A, x, b)
+    @test sum(abs2, x - [11.0/2.0, 55.0/4, 175.0/8.0]) < 1e-8
+
+    N = 100
+    A = spdiagm((2 * ones(N), -ones(N-1), -ones(N-1)), 
+                    (0, -1, 1), N, N) 
+    x = ones(eltype(A), N)
+    b = zeros(eltype(A), N)
+    s1 = GaussSeidel(ForwardSweep(), 200)
+    relax!(s1, A, x, b)
+    resid1 = norm(A*x,2)
+    x = ones(eltype(A), N)
+    s2 = GaussSeidel(BackwardSweep(), 200)
+    relax!(s2, A, x, b)
+    resid2 = norm(A*x,2)
+    @test resid1 < 0.01 && resid2 < 0.01
+    @test isapprox(resid1, resid2)
+    
 end