Fix Smoothed Aggregation corner cases and test

src/aggregate.jl

@@ -1,14 +1,21 @@
+"""
+    StandardAggregation()
+
+Implementation of Algorithm 5.1 from ,,Algebraic Multigrid by Smoothed
+Aggregation for Second and Fourth Order Elliptic Problems'' by Vanek et al. (1996).
+
+Note that isolated nodes are not aggregated.
+"""
 struct StandardAggregation
 end
 
 function (::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,73 +35,87 @@ 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
-
             for j in nzrange(S, i)
                 row = S.rowval[j]
-                x[row] = next_aggregate
+                if row != i
+                    x[row] = next_aggregate
+                end
             end
 
             next_aggregate += 1
         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
 
-    next_aggregate -= 1
+    # Record which nodes are still unaggregated after Pass 1 and Pass 2 *before*
+    # the 0-based shift below.  After the shift, accepted aggregate-0 nodes also
+    # have x[i] == 0, so we cannot use x[i] == 0 as an "unaggregated" sentinel
+    # inside the subsequent Pass 3 loop.
+    unagg = x .== 0
 
-    # Pass 3
+    # Shift all Pass 1 / Pass 2 assignments from 1-indexed to 0-indexed.
+    next_aggregate -= 1
     for i = 1:n
         xi = x[i]
-        if xi != 0
-            if xi > 0
-                x[i] = xi - 1
-            elseif xi == -n
-                x[i] = -1
-            else
-                x[i] = -xi - 1
-            end
-            continue
+        if xi > 0
+            x[i] = xi - 1
+        elseif xi == -n
+            x[i] = -1
+        elseif xi < 0  # Pass 2 tentative: x[i] = -x_best
+            x[i] = -xi - 1
         end
+        # unagg[i] == true nodes keep x[i] == 0; handled below.
+    end
+
+    # Pass 3: form new aggregates from nodes that were not reached in Pass 1 or 2.
+    # Each seed i pulls in every unaggregated neighbour (identified via `unagg`,
+    # not x[row] == 0, to avoid confusion with the 0-indexed aggregate-0).
+    for i = 1:n
+        unagg[i] || continue
 
         x[i] = next_aggregate
-        y[next_aggregate + 1] = i
 
         for j in nzrange(S, i)
             row = S.rowval[j]
-
-            if x[row] == 0
+            if unagg[row]
                 x[row] = next_aggregate
+                unagg[row] = false
             end
         end
-
+        unagg[i] = false
         next_aggregate += 1
     end
 
-    y = y[1:next_aggregate]
     M, N = (n, next_aggregate)
 
-    # Some nodes not aggregated
-    if minimum(x) == -1
+    # Aggregation of leftovers
+    if isempty(x) || minimum(x) == -1
         mask = x .!= -1
         I = collect(R, 1:n)[mask]
         J = x[mask] .+ R(1)

src/aggregation.jl

@@ -94,10 +94,12 @@ function smoothed_aggregation(A::TA,
     end
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
+        prev_n_levels = length(levels)
         @timeit_debug "extend_hierarchy!" A, B, bsr_flag = extend_hierarchy_sa!(levels, strength, aggregate, smooth,
                                 improve_candidates, diagonal_dominance,
                                 keep, A, B, presmoother, postsmoother, symmetry, bsr_flag, verbose)
-        size(A, 1) == 0 && break
+        # Break if no coarsening happened (all nodes isolated) or degenerate coarse matrix
+        length(levels) == prev_n_levels && break
         coarse_x!(w, size(A, 1))
         coarse_b!(w, size(A, 1))
         residual!(w, size(A, 1))
@@ -128,8 +130,11 @@ function extend_hierarchy_sa!(levels, strength, aggregate, smooth,
     # Aggregation operator
     @timeit_debug "aggregation" AggOp = aggregate(S)
 
+    # Cannot coarsen further if no aggregates were formed (e.g. all nodes isolated)
+    size(AggOp, 1) == 0 && return A, B, bsr_flag
+
     # Improve candidates
-    b = zeros(size(A,1),size(B,2))
+    b = zeros(eltype(A), size(A,1),size(B,2))
     @timeit_debug "improve candidates" improve_candidates(A, B, b)
     @timeit_debug "fit candidates" T, B = fit_candidates(AggOp, B)
 
@@ -160,19 +165,14 @@ function fit_candidates(AggOp, B::AbstractVector; tol=1e-10)
     n_col = n_coarse
 
     R = zeros(eltype(B), n_coarse)
-    Qx = zeros(eltype(B), nnz(A))
-    # copy!(Qx, B)
-    for i = 1:size(Qx, 1)
-        Qx[i] = B[i]
-    end
+    Qx = copy(B)
     # copy!(A.nzval, B)
     for i = 1:n_col
         for j in nzrange(A,i)
             row = A.rowval[j]
             A.nzval[j] = B[row]
         end
     end
-    k = 1
     for i = 1:n_col
         norm_i = norm_col(A, Qx, i)
         threshold_i = tol * norm_i
@@ -233,11 +233,7 @@ end
 function norm_col(A, Qx, i)
     s = zero(eltype(A))
     for j in nzrange(A, i)
-        if A.rowval[j] > length(Qx)
-            val = 1
-        else
-            val = Qx[A.rowval[j]]
-        end
+        val = Qx[A.rowval[j]]
         # val = A.nzval[A.rowval[j]]
         s += val*val
     end

test/nns_test.jl

@@ -73,6 +73,24 @@ end
         1e-20 .* hcat(ones(4), collect(0:3))
     ))
 
+    # 4. Isolated intermediate node: node 3 has no aggregate assignment.
+    # AggOp is (n_fine × n_agg): rows = fine nodes, cols = aggregates.
+    # The masking loop below zeros fine[3,:] before fit_candidates is called.
+    push!(cases, (
+        sparse([1, 2, 4, 5], [1, 1, 2, 2], ones(4), 5, 2),
+        hcat(ones(5), Float64[1, 2, 3, 4, 5])
+    ))
+    # Same with 3 near-null vectors so length(rows) < m for singleton aggregates
+    push!(cases, (
+        sparse([1, 2, 4, 5], [1, 1, 2, 2], ones(4), 5, 2),
+        hcat(ones(5), Float64[1,2,3,4,5], Float64[5,4,3,2,1])
+    ))
+    # Isolated nodes at both ends: nodes 1 and 7 isolated, 2 aggregates
+    push!(cases, (
+        sparse([2, 3, 4, 5, 6], [1, 1, 2, 2, 2], ones(5), 7, 2),
+        hcat(ones(7), Float64[1,2,3,4,5,6,7])
+    ))
+
     # Run tests
     for (AggOp, fine) in cases
         # mask dofs not in aggregation

test/sa_tests.jl

@@ -64,7 +64,7 @@ end
 function stand_agg(C, ϵ=0)
     n = size(C, 1)
 
-    # FIXME manouvering around the implementation begin CSC but pretending to be CSR and
+    # FIXME maneuvering around the implementation begin CSC but pretending to be CSR and
     #       the fact that the implementation can only handle symmetric matrices while
     #       this test covers has non-symmetric ones...
     NϵT(C, i, ϵ) = [j for j in 1:size(C, 2) if abs(C[i, j]) > ϵ * sqrt(C[i,i]*C[j,j])]
@@ -117,19 +117,14 @@ function stand_agg(C, ϵ=0)
             continue
         end
 
-        Ni = union(Set(Nϵ(C, i, ϵ)), R)
-        # Do not aggregate isolated nodes
-        if length(Ni) > 0
-            continue
-        end
-        j += 1
-        push!(Cpts, Ni)
-
+        Ni = intersect(Set(Nϵ(C, i, ϵ)), R)
+        push!(Ni, i)  # always include the seed node itself
+        push!(Cpts, i)
+        setdiff!(R, Ni)
         for x in Ni
-            if x in R
-                aggregates[x] = j
-            end
+            aggregates[x] = j
         end
+        j += 1
     end
 
     mask = aggregates .> -1
@@ -139,6 +134,55 @@ function stand_agg(C, ϵ=0)
     return sparse(I, J, V, maximum(I; init=0), n)
 end
 
+
+
+# Corner case tests for StandardAggregation
+function test_standard_aggregation_corner_cases()
+
+    # Off-by-one in aggregate_size: a 4-node chain whose strength matrix has no
+    # diagonal entries should form 2 aggregates of size 2.
+    S_chain = sparse([1,2,2,3,3,4],[2,1,3,2,4,3], ones(Float64,6), 4, 4)
+    AggOp_chain = StandardAggregation()(S_chain)
+    @test size(AggOp_chain, 1) == 2                          # exactly 2 aggregates
+    @test all(vec(sum(AggOp_chain, dims=1)) .== 1)           # every node in exactly one
+
+    # Disconnected graph: two independent 3-node chains.
+    # Both components must be fully aggregated without mixing.
+    rows_d = [1,2,2,3,4,5,5,6]; cols_d = [2,1,3,2,5,4,6,5]
+    S_disc = sparse(rows_d, cols_d, ones(Float64,8), 6, 6) + spdiagm(0 => ones(6))
+    ref_disc  = stand_agg(S_disc)
+    calc_disc = StandardAggregation()(S_disc)
+    @test sum(abs2, calc_disc - ref_disc) < 1e-6
+
+    # All nodes isolated (diagonal-only strength matrix) → nothing is aggregated.
+    S_iso = spdiagm(0 => ones(Float64, 5))
+    @test sum(abs2, StandardAggregation()(S_iso) - stand_agg(S_iso)) < 1e-6
+    @test nnz(StandardAggregation()(S_iso)) == 0
+
+    # Empty 0×0 strength matrix must not crash (guards minimum() on empty array).
+    S_empty = spzeros(Float64, 0, 0)
+    AggOp_empty = StandardAggregation()(S_empty)
+    @test size(AggOp_empty) == (0, 0)
+
+    # smoothed_aggregation on a large diagonal matrix (all nodes isolated at every level)
+    # must return a valid 1-level hierarchy rather than crashing at solve time.
+    A_diag = spdiagm(0 => 2.0 * ones(Float64, 20))
+    ml_diag = smoothed_aggregation(A_diag)
+    @test length(ml_diag) == 1
+    @test size(ml_diag.final_A) == (20, 20)
+
+    # Intermediate isolated node: a 5-node chain where node 3 has a large diagonal
+    # that severs the connection between aggregates {1,2} and {4,5}.
+    # Verifies AggOp has a zero column for node 3 and the full pipeline solves.
+    A_iso = spdiagm(-1 => fill(-0.5, 4), 0 => [1.0,1.0,100.0,1.0,1.0],
+                                          1 => fill(-0.5, 4))
+    S_iso5, _ = SymmetricStrength(0.25)(A_iso)
+    AggOp_iso5 = StandardAggregation()(S_iso5)
+    @test size(AggOp_iso5, 1) == 2              # exactly 2 aggregates
+    @test nnz(AggOp_iso5[:, 3]) == 0           # node 3 isolated (zero column)
+
+end
+
 # Standard aggregation tests
 function test_standard_aggregation()
 
@@ -200,6 +244,24 @@ function generate_fit_candidates_cases()
                         [3,2,1,1,2,3,2,1,3], ones(T,9))
         B = T.(collect(1:9))
         push!(cases, (AggOp, B))
+
+        # Isolated intermediate node: fine node 3 has no aggregate (zero column).
+        # AggOp is 2×5, colptr has zero-width at column 3.
+        # Regression test for Qx allocation bug: for non-unit B the norm of
+        # aggregate 2 (nodes {4,5}) was computed incorrectly when fine-node
+        # indices exceeded nnz(AggOp).
+        AggOp = SparseMatrixCSC(2, 5, Int[1,2,3,3,4,5],
+                        Int[1,1,2,2], ones(T,4))
+        B = T.([1, 1, 5, 2, 3])
+        push!(cases, (AggOp, B))
+
+        # Two isolated intermediate nodes: fine nodes 3 and 7 have no aggregate.
+        # Aggregates: agg 1 = {nodes 1,2}, agg 2 = {nodes 4,5,6}, agg 3 = {nodes 8,9}.
+        # AggOp is 3×9 with zero-width columns at positions 3 and 7.
+        AggOp = SparseMatrixCSC(3, 9, Int[1,2,3,3,4,5,6,6,7,8],
+                        Int[1,1,2,2,2,3,3], ones(T,7))
+        B = T.(collect(1:9))
+        push!(cases, (AggOp, B))
     end
 
     cases
@@ -325,6 +387,14 @@ function test_jacobi_prolongator()
     @test sum(abs2, x - ref) < 1e-6
 end
 
+# Issue #24
+function nodes_not_agg()
+    A = include("onetoall.jl")
+    ml = smoothed_aggregation(A)
+    @test size(ml.levels[2].A) == (11,11)
+    @test size(ml.final_A) == (2,2)
+end
+
 # Smoothed Aggregation
 @testset "Smoothed Aggregation" begin
     @testset "Symmetric Strength of Connection" begin
@@ -335,6 +405,10 @@ end
         test_standard_aggregation()
     end
 
+    @testset "Standard Aggregation Corner Cases" begin
+        test_standard_aggregation_corner_cases()
+    end
+
     @testset "Fit Candidates" begin
         test_fit_candidates()
     end
@@ -351,4 +425,9 @@ end
         a = sparse(Int32.(1:10), Int32.(1:10), rand(10))
         @inferred smoothed_aggregation(a)
     end
+
+    # Issue #24
+    @testset "Regression Test #24" begin
+        nodes_not_agg()
+    end
 end

test/test_complex.jl

@@ -2,6 +2,7 @@
     A = poisson((5, 5))
     Ac = A .* 1/√2 + A .* im/√2
 
+    Random.seed!(1337)
     u = rand(Complex{Float64}, 5*5)
     b = Ac*u