Fix smoothers and add a scratch

Project.toml

@@ -13,7 +13,9 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 CommonSolve = "0.2"
+JLD2 = "0.6.3"
 LinearSolve = "2, 3"
+MatrixDepot = "1.0.14"
 Reexport = "1"
 TimerOutputs = "0.5.29"
 julia = "1.6"
@@ -24,8 +26,9 @@ FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test"]
+test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "MatrixDepot", "Random", "Test"]

README.md

@@ -82,7 +82,7 @@ This package currently supports:
 
 AMG Styles:
 * Ruge-Stuben Solver
-* Smoothed Aggregation (SA)
+* Smoothed Aggregation (SA) with handling of user-provided near null spaces
 
 Strength of Connection:
 * Classical Strength of Connection
@@ -91,14 +91,16 @@ Strength of Connection:
 Smoothers:
 * Gauss Seidel (Symmetric, Forward, Backward)
 * Damped Jacobi
+* SOR and SSOR
 
 Cycling:
 * V, W and F cycles
 
 In the future, this package will support:
 1. Other splitting methods (like CLJP)
-2. SOR smoother
-3. AMLI cycles
+2. AMLI cycles
+3. Root-Node AMG
+4. Shared-memory parallelization of relevant internal functions
 
 #### Acknowledgements
 This package has been heavily inspired by the [`PyAMG`](http://github.com/pyamg/pyamg) project.

src/aggregate.jl

@@ -1,14 +1,17 @@
 struct StandardAggregation
+    min_aggregate_size::Int
 end
 
-function (::StandardAggregation)(S::SparseMatrixCSC{T,R}) where {T,R}
+StandardAggregation() = StandardAggregation(1)
 
+function (agg::StandardAggregation)(S::SparseMatrixCSC{T,R}) where {T,R}
     n = size(S, 1)
     x = zeros(R, n)
     y = zeros(R, n)
 
     next_aggregate = 1
 
+    # Pass 1: Tentative aggregation
     for i = 1:n
         if x[i] != 0
             continue
@@ -28,41 +31,56 @@ function (::StandardAggregation)(S::SparseMatrixCSC{T,R}) where {T,R}
             end
         end
 
-        if !has_neighbors
+        if !has_neighbors # Mark isolated node
             x[i] = -n
         elseif !has_agg_neighbors
             x[i] = next_aggregate
-            y[next_aggregate] = i
-
+            aggregate_size = 0
             for j in nzrange(S, i)
                 row = S.rowval[j]
                 x[row] = next_aggregate
+                aggregate_size += 1
             end
 
-            next_aggregate += 1
+            # Reject aggregate if it is too small
+            if aggregate_size < agg.min_aggregate_size
+                x[i] = 0
+                for j in nzrange(S, i)
+                    row = S.rowval[j]
+                    x[row] = 0
+                end
+            else
+                y[next_aggregate] = i
+                next_aggregate += 1
+            end
         end
     end
 
-
-    # Pass 2
+    # Pass 2: Enlarge tentative aggregates
     for i = 1:n
+        # Skip marked node
         if x[i] != 0
             continue
         end
 
+        s_best = zero(eltype(S))
+        x_best = 0
         for j in nzrange(S, i)
             row = S.rowval[j]
             x_row = x[row]
-            if x_row > 0
-                x[i] = -x_row
-                break
+            s_candidate = S.nzval[j]
+            if x_row > 0 && s_candidate > s_best # Assigned and stronger than previous best
+                s_best = s_candidate
+                x_best = x_row
             end
         end
+        if x_best > 0
+            x[i] = -x_best
+        end
     end
 
+    # Shift assignments by 1 and apply the temporary assignments
     next_aggregate -= 1
-
-    # Pass 3
     for i = 1:n
         xi = x[i]
         if xi != 0
@@ -91,9 +109,9 @@ function (::StandardAggregation)(S::SparseMatrixCSC{T,R}) where {T,R}
     end
 
     y = y[1:next_aggregate]
-    M,N = (n, next_aggregate)
+    M, N = (n, next_aggregate)
 
-    # Some nodes not aggregated
+    # Pass 3: Aggregation of leftovers
     if minimum(x) == -1
         mask = x .!= -1
         I = collect(R, 1:n)[mask]

src/aggregation.jl

@@ -12,6 +12,7 @@ function smoothed_aggregation(A::TA,
                         max_coarse = 10,
                         diagonal_dominance = false,
                         keep = false,
+                        verbose = false,
                         coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}}
 
     @timeit_debug "prologue" begin
@@ -30,25 +31,30 @@ function smoothed_aggregation(A::TA,
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
         @timeit_debug "extend_hierarchy!" A, B, bsr_flag = extend_hierarchy_sa!(levels, strength, aggregate, smooth,
                                 improve_candidates, diagonal_dominance,
-                                keep, A, B, symmetry, bsr_flag)
+                                keep, A, B, presmoother, postsmoother, symmetry, bsr_flag, verbose)
         size(A, 1) == 0 && break
         coarse_x!(w, size(A, 1))
         coarse_b!(w, size(A, 1))
         residual!(w, size(A, 1))
     end
 
     @timeit_debug "coarse solver setup" cs = coarse_solver(A)
-    @timeit_debug "ml setup" ml = MultiLevel(levels, A, cs, presmoother, postsmoother, w)
+    @timeit_debug "ml setup" ml = MultiLevel(levels, A, cs, w)
+
+    if verbose
+        @info ml
+    end
+
     return ml
 end
 
 struct HermitianSymmetry
 end
 
 function extend_hierarchy_sa!(levels, strength, aggregate, smooth,
-                            improve_candidates, diagonal_dominance, keep,
-                            A, B,
-                            symmetry, bsr_flag)
+                            improve_candidates, diagonal_dominance,
+                            keep, A, B, presmoother, postsmoother,
+                            symmetry, bsr_flag, verbose = false)
 
     # Calculate strength of connection matrix
     @timeit_debug "strength" if symmetry isa HermitianSymmetry
@@ -59,7 +65,6 @@ function extend_hierarchy_sa!(levels, strength, aggregate, smooth,
 
     # Aggregation operator
     @timeit_debug "aggregation" AggOp = aggregate(S)
-    # b = zeros(eltype(A), size(A, 1))
 
     # Improve candidates
     b = zeros(size(A,1),size(B,2))
@@ -73,16 +78,20 @@ function extend_hierarchy_sa!(levels, strength, aggregate, smooth,
 
     @timeit_debug "RAP" RAP = R * A * P
 
-    push!(levels, Level(A, P, R))
+    @timeit_debug "smoother setup" begin
+        pre = setup_smoother(presmoother, A)
+        post = setup_smoother(postsmoother, A)
+        push!(levels, Level(A, P, R, pre, post))
+    end
 
     bsr_flag = true
 
+    # RAP is the coarse matrix and B is the coarse null space
     RAP, B, bsr_flag
 end
 construct_R(::HermitianSymmetry, P) = P'
 
 function fit_candidates(AggOp, B::AbstractVector; tol=1e-10)
-
     A = adjoint(AggOp)
     n_fine, n_coarse = size(A)
     n_col = n_coarse
@@ -113,10 +122,6 @@ function fit_candidates(AggOp, B::AbstractVector; tol=1e-10)
         end
         for j in nzrange(A, i)
             row = A.rowval[j]
-            # Qx[row] *= scale
-            #@show k
-            # Qx[k] *= scale
-            # k += 1
             A.nzval[j] *= scale
         end
     end
@@ -127,12 +132,12 @@ end
 
 function fit_candidates(AggOp, B::AbstractMatrix; tol=1e-10)
     A = adjoint(AggOp)
-    n_fine, m = ndims(B) == 1 ? (length(B), 1) : size(B)
+    n_fine, m = size(B)
     n_fine2, n_agg = size(A)
     @assert n_fine2 == n_fine
     n_coarse = m * n_agg
     T = eltype(B)
-    Qs = spzeros(T, n_fine, n_coarse)
+    Qs = spzeros(T, n_fine, n_coarse) # TODO use CSR here as the algorithm becomes easier
     R = zeros(T, n_coarse, m)
 
     for agg in 1:n_agg
@@ -154,11 +159,11 @@ function fit_candidates(AggOp, B::AbstractMatrix; tol=1e-10)
                 Qs[rows[local_i], offset+local_j] = val
             end
         end
-        dropzeros!(Qs)
 
         R[offset+1:offset+r, :] .= Rj[1:r, :]
     end
 
+    dropzeros!(Qs)
     return Qs, R
 end
 

src/classical.jl

@@ -25,17 +25,17 @@ function ruge_stuben(_A::Union{TA, Symmetric{Ti, TA}, Hermitian{Ti, TA}},
     residual!(w, size(A, 1))
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
-        @timeit_debug "extend_hierarchy!" A = extend_hierarchy_rs!(levels, strength, CF, A, symmetric)
+        @timeit_debug "extend_hierarchy!" A = extend_hierarchy_rs!(levels, strength, CF, A, presmoother, postsmoother, symmetric)
         coarse_x!(w, size(A, 1))
         coarse_b!(w, size(A, 1))
         residual!(w, size(A, 1))
     end
 
     @timeit_debug "coarse solver setup" cs = coarse_solver(A)
-    return MultiLevel(levels, A, cs, presmoother, postsmoother, w)
+    return MultiLevel(levels, A, cs, w)
 end
 
-function extend_hierarchy_rs!(levels, strength, CF, A::SparseMatrixCSC{Ti,Tv}, symmetric) where {Ti,Tv}
+function extend_hierarchy_rs!(levels, strength, CF, A::SparseMatrixCSC{Ti,Tv}, presmoother, postsmoother, symmetric) where {Ti,Tv}
     if symmetric
         At = A
     else
@@ -45,7 +45,13 @@ function extend_hierarchy_rs!(levels, strength, CF, A::SparseMatrixCSC{Ti,Tv}, s
     @timeit_debug "splitting" splitting = CF(S)
     @timeit_debug "interpolation" P, R = direct_interpolation(At, T, splitting)
     @timeit_debug "RAP" RAP = R * A * P
-    push!(levels, Level(A, P, R))
+
+    @timeit_debug "smoother setup" begin
+        pre = setup_smoother(presmoother, A)
+        post = setup_smoother(postsmoother, A)
+        push!(levels, Level(A, P, R, pre, post))
+    end
+
     return RAP
 end
 

src/multilevel.jl

@@ -1,15 +1,19 @@
-struct Level{TA, TP, TR}
+struct Level{TA, TP, TR, TPre, TPost}
     A::TA
     P::TP
     R::TR
+    presmoother::TPre
+    postsmoother::TPost
 end
 
-struct MultiLevel{S, Pre, Post, TA, TP, TR, TW}
+function Base.show(io::IO, l::Level)
+    print(io, "Level with R $(size(l.R)) | A $(size(l.A)) | P $(size(l.P))")
+end
+
+struct MultiLevel{S, TA, TP, TR, TW}
     levels::Vector{Level{TA, TP, TR}}
     final_A::TA
     coarse_solver::S
-    presmoother::Pre
-    postsmoother::Post
     workspace::TW
 end
 
@@ -259,7 +263,7 @@ end
 
 function __solve!(x, ml, cycle::Cycle, b, lvl)
     A = ml.levels[lvl].A
-    @timeit_debug "Presmoother" ml.presmoother(A, x, b)
+    @timeit_debug "Presmoother" ldiv!(x, ml.levels[lvl].presmoother, b)
 
     res = ml.workspace.res_vecs[lvl]
     @timeit_debug "Residual eval" mul!(res, A, x)
@@ -279,7 +283,7 @@ function __solve!(x, ml, cycle::Cycle, b, lvl)
     @timeit_debug "Prolongation" mul!(res, ml.levels[lvl].P, coarse_x)
     x .+= res
 
-    @timeit_debug "Postsmoother" ml.postsmoother(A, x, b)
+    @timeit_debug "Postsmoother" ldiv!(x, ml.levels[lvl].postsmoother, b)
 
     x
 end