use callable structs

mohamed82008 · mohamed82008 · commit e89a333f6774 · 2018-08-12T11:22:15.000+02:00
diff --git a/src/AlgebraicMultigrid.jl b/src/AlgebraicMultigrid.jl
@@ -14,17 +14,17 @@ include("utils.jl")
 export approximate_spectral_radius
 
 include("strength.jl")
-export strength_of_connection, Classical, SymmetricStrength
+export Classical, SymmetricStrength
 
 include("splitting.jl")
-export split_nodes, RS
+export RS
 
 include("gallery.jl")
 export poisson
 
 include("smoother.jl")
 export GaussSeidel, SymmetricSweep, ForwardSweep, BackwardSweep,
-        smooth_prolongator, JacobiProlongation
+        JacobiProlongation
 
 include("multilevel.jl")
 export solve
@@ -33,7 +33,7 @@ include("classical.jl")
 export ruge_stuben
 
 include("aggregate.jl")
-export aggregation, StandardAggregation
+export StandardAggregation
 
 include("aggregation.jl")
 export fit_candidates, smoothed_aggregation
diff --git a/src/aggregate.jl b/src/aggregate.jl
@@ -1,7 +1,7 @@
 struct StandardAggregation
 end
 
-function aggregation(::StandardAggregation, S::SparseMatrixCSC{T,R}) where {T,R}
+function (::StandardAggregation)(S::SparseMatrixCSC{T,R}) where {T,R}
 
     n = size(S, 1)
     x = zeros(R, n)
diff --git a/src/aggregation.jl b/src/aggregation.jl
@@ -54,18 +54,18 @@ function extend_hierarchy!(levels, strength, aggregate, smooth,
                             symmetry, bsr_flag)
 
     # Calculate strength of connection matrix
-    S = strength_of_connection(strength, A, bsr_flag)
+    S = strength(A, bsr_flag)
 
     # Aggregation operator
-    AggOp = aggregation(aggregate, S)
+    AggOp = aggregate(S)
     # b = zeros(eltype(A), size(A, 1))
 
     # Improve candidates
     b = zeros(size(A,1))
-    relax!(improve_candidates, A, B, b)
+    improve_candidates(A, B, b)
     T, B = fit_candidates(AggOp, B)
 
-    P = smooth_prolongator(smooth, A, T, S, B)
+    P = smooth(A, T, S, B)
     R = construct_R(symmetry, P)
     push!(levels, Level(A, P, R))
 
diff --git a/src/classical.jl b/src/classical.jl
@@ -28,8 +28,8 @@ end
 
 function extend_heirarchy!(levels::Vector{Level{Ti,Tv}}, strength, CF, A::SparseMatrixCSC{Ti,Tv}) where {Ti,Tv}
     At = copy(A')
-    S, T = strength_of_connection(strength, At)
-    splitting = split_nodes(CF, S)
+    S, T = strength(At)
+    splitting = CF(S)
     P, R = direct_interpolation(At, T, splitting)
     push!(levels, Level(A, P, R))
     A = R * A * P
diff --git a/src/multilevel.jl b/src/multilevel.jl
@@ -130,7 +130,7 @@ end
 function __solve(v::V, ml, x, b, lvl)
 
     A = ml.levels[lvl].A
-    presmoother!(ml.presmoother, A, x, b)
+    ml.presmoother(A, x, b)
 
     res = b - A * x
     coarse_b = ml.levels[lvl].R * res
@@ -144,7 +144,7 @@ function __solve(v::V, ml, x, b, lvl)
 
     x .+= ml.levels[lvl].P * coarse_x
 
-    postsmoother!(ml.postsmoother, A, x, b)
+    ml.postsmoother(A, x, b)
 
     x
 end
diff --git a/src/preconditioner.jl b/src/preconditioner.jl
@@ -1,5 +1,3 @@
-import Compat.LinearAlgebra: \, *, ldiv!, mul!
-
 struct Preconditioner
     ml::MultiLevel
 end
diff --git a/src/smoother.jl b/src/smoother.jl
@@ -15,30 +15,17 @@ GaussSeidel(f::ForwardSweep) = GaussSeidel(f, 1)
 GaussSeidel(b::BackwardSweep) = GaussSeidel(b, 1)
 GaussSeidel(s::SymmetricSweep) = GaussSeidel(s, 1)
 
-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)
-
-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
-        gs!(A, b, x, 1, 1, size(A, 1))
-        gs!(A, b, x, size(A,1), -1, 1)
-    end
-end
-
-function smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b)
+function (s::GaussSeidel{S})(A, x, b) where {S<:Sweep}
     for i in 1:s.iter
-        gs!(A, b, x, size(A,1), -1, 1)
+        if S === ForwardSweep || S === SymmetricSweep
+            gs!(A, b, x, 1, 1, size(A, 1))
+        end
+        if S === BackwardSweep || S === SymmetricSweep
+            gs!(A, b, x, size(A, 1), -1, 1)
+        end
     end
 end
 
-
 function gs!(A, b, x, start, step, stop)
     n = size(A, 1)
     z = zero(eltype(A))
@@ -104,10 +91,7 @@ end
 struct LocalWeighting
 end
 
-function smooth_prolongator(j::JacobiProlongation,
-                                        A, T, S, B,
-                                        degree = 1,
-                                        weighting = LocalWeighting())
+function (j::JacobiProlongation)(A, T, S, B, degree = 1, weighting = LocalWeighting())
     D_inv_S = weight(weighting, A, j.ω)
     P = T
     for i = 1:degree
diff --git a/src/splitting.jl b/src/splitting.jl
@@ -17,7 +17,7 @@ function remove_diag!(a)
 	dropzeros!(a)
 end
 
-function split_nodes(::RS, S)
+function (::RS)(S)
 	remove_diag!(S)
 	RS_CF_splitting(S, copy(S'))
 end
diff --git a/src/strength.jl b/src/strength.jl
@@ -4,8 +4,7 @@ struct Classical{T} <: Strength
 end
 Classical(;θ = 0.25) = Classical(θ)
 
-function strength_of_connection(c::Classical, 
-                At::SparseMatrixCSC{Tv,Ti}) where {Ti,Tv}
+function (c::Classical)(At::SparseMatrixCSC{Tv,Ti}) where {Ti,Tv}
 
     θ = c.θ
 
@@ -75,7 +74,7 @@ struct SymmetricStrength{T} <: Strength
 end
 SymmetricStrength() = SymmetricStrength(0.)
 
-function strength_of_connection(s::SymmetricStrength{T}, A, bsr_flag = false) where {T}
+function (s::SymmetricStrength{T})(A, bsr_flag = false) where {T}
 
     θ = s.θ
 
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -16,13 +16,13 @@ ref_split = readdlm("ref_split_test.txt")
 # classical strength of connection
 A = poisson(5)
 A = float.(A)
-S, T = strength_of_connection(Classical(0.2), A)
+S, T = Classical(0.2)(A)
 @test Matrix(S) == [ 1.0  0.5  0.0  0.0  0.0
                    0.5  1.0  0.5  0.0  0.0
                    0.0  0.5  1.0  0.5  0.0
                    0.0  0.0  0.5  1.0  0.5
                    0.0  0.0  0.0  0.5  1.0 ]
-S, T = strength_of_connection(Classical(0.25), graph)
+S, T = Classical(0.25)(graph)
 diff = S - ref_S
 @test maximum(diff) < 1e-10
 
@@ -32,18 +32,18 @@ end
 
 # Ruge-Stuben splitting
 S = poisson(7)
-@test split_nodes(RS(), S) == [0, 1, 0, 1, 0, 1, 0]
+@test RS()(S) == [0, 1, 0, 1, 0, 1, 0]
 srand(0)
 S = sprand(10,10,0.1); S = S + S'
-@test split_nodes(RS(), S) ==  [0, 1, 1, 0, 0, 0, 0, 0, 1, 1]
+@test RS()(S) ==  [0, 1, 1, 0, 0, 0, 0, 0, 1, 1]
 
 a = load("thing.jld2")["G"]
-S, T = AlgebraicMultigrid.strength_of_connection(Classical(0.25), a)
-@test split_nodes(RS(), S) == [0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
+S, T = Classical(0.25)(a)
+@test RS()(S) == [0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
 0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0,
 1, 0]
 
-@test split_nodes(RS(), ref_S) == Int.(vec(ref_split))
+@test RS()(ref_S) == Int.(vec(ref_split))
 
 end
 
diff --git a/test/sa_tests.jl b/test/sa_tests.jl
@@ -1,4 +1,4 @@
-import AlgebraicMultigrid: scale_cols_by_largest_entry!, strength_of_connection, 
+import AlgebraicMultigrid: scale_cols_by_largest_entry!, 
             SymmetricStrength, poisson
 function symmetric_soc(A::SparseMatrixCSC{T,V}, θ) where {T,V}
     D = abs.(diag(A))
@@ -30,7 +30,7 @@ function test_symmetric_soc()
     for matrix in cases
         for θ in (0.0, 0.1, 0.5, 1., 10.)
             ref_matrix = symmetric_soc(matrix, θ)
-            calc_matrix = strength_of_connection(SymmetricStrength(θ), matrix)
+            calc_matrix = SymmetricStrength(θ)(matrix)
 
             @test sum(abs2, ref_matrix - calc_matrix) < 1e-6
         end
@@ -128,7 +128,7 @@ function test_standard_aggregation()
     for matrix in cases
         for θ in (0.0, 0.1, 0.5, 1., 10.)
             C = symmetric_soc(matrix, θ)
-            calc_matrix = aggregation(StandardAggregation(), matrix)
+            calc_matrix = StandardAggregation()(matrix)
             ref_matrix = stand_agg(matrix)
             @test sum(abs2, ref_matrix - calc_matrix) < 1e-6
         end
@@ -235,66 +235,66 @@ function test_approximate_spectral_radius()
 end
 
 # Test Gauss Seidel 
-import AlgebraicMultigrid: gs!, relax!
+import AlgebraicMultigrid: gs!
 function test_gauss_seidel()    
     N = 1
     A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
     b = zeros(N)
     s = GaussSeidel(ForwardSweep())
-    relax!(s, A, x, b)
+    s(A, x, b)
     @test sum(abs2, x - zeros(N)) < 1e-8
 
     N = 3 
     A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
     b = zeros(N)
     s = GaussSeidel(ForwardSweep())
-    relax!(s, A, x, b)
+    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(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
     b = zeros(N)
     s = GaussSeidel(BackwardSweep())
-    relax!(s, A, x, b)
+    s(A, x, b)
     @test sum(abs2, x - zeros(N)) < 1e-8
 
     N = 3 
     A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
     b = zeros(N)
     s = GaussSeidel(BackwardSweep())
-    relax!(s, A, x, b)
+    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(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
     b = eltype(A).([10.])
     s = GaussSeidel(ForwardSweep())
-    relax!(s, A, x, b)
+    s(A, x, b)
     @test sum(abs2, x - [5.]) < 1e-8
 
     N = 3 
     A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
     b = eltype(A).([10., 20., 30.])
     s = GaussSeidel(ForwardSweep())
-    relax!(s, A, x, b)
+    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(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = ones(eltype(A), N)
     b = zeros(eltype(A), N)
     s1 = GaussSeidel(ForwardSweep(), 200)
-    relax!(s1, A, x, b)
+    s1(A, x, b)
     resid1 = norm(A*x,2)
     x = ones(eltype(A), N)
     s2 = GaussSeidel(BackwardSweep(), 200)
-    relax!(s2, A, x, b)
+    s2(A, x, b)
     resid2 = norm(A*x,2)
     @test resid1 < 0.01 && resid2 < 0.01
     @test isapprox(resid1, resid2)
@@ -305,7 +305,7 @@ end
 function test_jacobi_prolongator()
     A = poisson(100)
     T = poisson(100)
-    x = smooth_prolongator(JacobiProlongation(4/3), A, T, 1, 1)
+    x = JacobiProlongation(4/3)(A, T, 1, 1)
     ref = load("ref_R.jld2")["G"]
     @test sum(abs2, x - ref) < 1e-6
 end
@@ -319,13 +319,12 @@ function nodes_not_agg()
 end
 
 # Issue 26
-import AlgebraicMultigrid: relax!
 function test_symmetric_sweep()
     A = poisson(10)
     s = GaussSeidel(SymmetricSweep(), 4)
     x = ones(size(A,1))
     b = zeros(size(A,1))
-    relax!(s, A, x, b)
+    s(A, x, b)
     @test sum(abs2, x - [0.176765; 0.353529; 0.497517; 0.598914; 
                             0.653311; 0.659104; 0.615597; 0.52275; 
                             0.382787; 0.203251]) < 1e-6        

Original file line number	Diff line number	Diff line change
`@@ -1,5 +1,3 @@`
`1`		`-import Compat.LinearAlgebra: \, *, ldiv!, mul!`
`2`		`-`
`3`	`1`	`struct Preconditioner`
`4`	`2`	`ml::MultiLevel`
`5`	`3`	`end`