diff --git a/src/aggregate.jl b/src/aggregate.jl index 5252e86..1eeae6f 100644 --- a/src/aggregate.jl +++ b/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) diff --git a/src/aggregation.jl b/src/aggregation.jl index 29c28cb..a05dcb9 100644 --- a/src/aggregation.jl +++ b/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,11 +165,7 @@ 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) @@ -172,7 +173,6 @@ function fit_candidates(AggOp, B::AbstractVector; tol=1e-10) 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 diff --git a/test/nns_test.jl b/test/nns_test.jl index 24b6f56..19d32e4 100644 --- a/test/nns_test.jl +++ b/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 diff --git a/test/sa_tests.jl b/test/sa_tests.jl index e0c5f37..af93638 100644 --- a/test/sa_tests.jl +++ b/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 diff --git a/test/test_complex.jl b/test/test_complex.jl index 4656e0f..c440ad0 100644 --- a/test/test_complex.jl +++ b/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