Merge pull request #47 from JuliaLinearAlgebra/bounds

Fix bug on small linear systems

src/aggregation.jl

@@ -35,9 +35,9 @@ function smoothed_aggregation(A::TA,
     end
     bsr_flag = false
     w = MultiLevelWorkspace(Val{bs}, eltype(A))
+    residual!(w, size(A, 1))
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
-        residual!(w, size(A, 1))
         A, B, bsr_flag = extend_hierarchy!(levels, strength, aggregate, smooth,
                                 improve_candidates, diagonal_dominance,
                                 keep, A, B, symmetry, bsr_flag)
@@ -46,6 +46,7 @@ function smoothed_aggregation(A::TA,
         #=if size(A, 1) <= max_coarse
             break
         end=#
+        residual!(w, size(A, 1))
     end
     #=A, B = extend_hierarchy!(levels, strength, aggregate, smooth,
                             improve_candidates, diagonal_dominance,

src/classical.jl

@@ -40,12 +40,13 @@ function ruge_stuben(_A::Union{TA, Symmetric{Ti, TA}, Hermitian{Ti, TA}},
         end
     end
     w = MultiLevelWorkspace(Val{bs}, eltype(A))
+    residual!(w, size(A, 1))
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
-        residual!(w, size(A, 1))
         A = extend_heirarchy!(levels, strength, CF, A, symmetric)
         coarse_x!(w, size(A, 1))
         coarse_b!(w, size(A, 1))
+        residual!(w, size(A, 1))
     end
 
     MultiLevel(levels, A, coarse_solver(A), presmoother, postsmoother, w)

test/REQUIRE

@@ -0,0 +1,2 @@
+JLD2
+FileIO

test/runtests.jl

@@ -2,6 +2,8 @@ using Compat, Compat.Test, Compat.LinearAlgebra
 using Compat.SparseArrays, Compat.DelimitedFiles, Compat.Random
 using IterativeSolvers, AlgebraicMultigrid
 import AlgebraicMultigrid: Pinv, Classical
+using JLD2
+using FileIO
 
 if VERSION < v"0.7-"
     const seed! = srand
@@ -304,4 +306,18 @@ for sz in [10, 5, 2]
     @test AlgebraicMultigrid.grid_complexity(ml) == 1
 end
 
+# Issue #46
+for f in (smoothed_aggregation, ruge_stuben)
+
+    a = load("bug.jld2")["G"]
+    ml = f(a)
+    p = aspreconditioner(ml)
+    b = zeros(size(a,1))
+    b[1] = 1
+    b[2] = -1
+    @test sum(abs2, a * solve(ml, b) - b) < 1e-10
+    @test sum(abs2, a * cg(a, b, Pl = p, maxiter = 1000) - b) < 1e-10
+
+end
+
 end