Merge pull request #7 from ranjanan/morefix

Fix Interpolation bugs

Author: ranjanan

src/classical.jl

@@ -22,11 +22,11 @@ function ruge_stuben(A::SparseMatrixCSC;
 
     while length(levels) < max_levels
         A = extend_heirarchy!(levels, strength, CF, A)
-        if size(levels[end].A, 1) < max_coarse
+        if size(A, 1) < max_coarse
             break
         end
     end
-    MultiLevel(levels, presmoother, postsmoother)
+    MultiLevel(levels, A, presmoother, postsmoother)
 end
 
 function extend_heirarchy!(levels::Vector{Level}, strength, CF, A)
@@ -44,7 +44,7 @@ function direct_interpolation{T,V}(A::T, S::T, splitting::Vector{V})
     Pp = rs_direct_interpolation_pass1(S, A, splitting)
     Pp = Pp .+ 1
 
-    Px, Pj = rs_direct_interpolation_pass2(A, S, splitting, Pp)
+    Px, Pj, Pp = rs_direct_interpolation_pass2(A, S, splitting, Pp)
 
     # Px .= abs.(Px)
     Pj = Pj .+ 1
@@ -59,6 +59,7 @@ end
 function rs_direct_interpolation_pass1(S, A, splitting)
 
      Bp = zeros(Int, size(A.colptr))
+     T = S'
      #=Sp = S.colptr
      Sj = S.rowval
      n_nodes = size(A, 1)
@@ -82,8 +83,8 @@ function rs_direct_interpolation_pass1(S, A, splitting)
          if splitting[i] == C_NODE
              nnz += 1
          else
-             for j in nzrange(S, i)
-                 row = S.rowval[j]
+             for j in nzrange(T, i)
+                 row = T.rowval[j]
                  if splitting[row] == C_NODE && row != i
                      nnz += 1
                 end
@@ -101,8 +102,9 @@ function rs_direct_interpolation_pass1(S, A, splitting)
                                                 Bp::Vector{Ti})
 
 
-    Bx = zeros(Float64, Bp[end])
-    Bj = zeros(Ti, Bp[end])
+    T = S'
+    Bx = zeros(Float64, Bp[end] - 1)
+    Bj = zeros(Ti, Bp[end] - 1)
 
     n = size(A, 1)
 
@@ -113,9 +115,9 @@ function rs_direct_interpolation_pass1(S, A, splitting)
         else
             sum_strong_pos = zero(Tv)
             sum_strong_neg = zero(Tv)
-            for j in nzrange(S, i)
-                row = S.rowval[j]
-                sval = S.nzval[j]
+            for j in nzrange(T, i)
+                row = T.rowval[j]
+                sval = T.nzval[j]
                 if splitting[row] == C_NODE && row != i
                     if sval < 0
                         sum_strong_neg += sval
@@ -153,9 +155,9 @@ function rs_direct_interpolation_pass1(S, A, splitting)
 
             nnz = Bp[i]
 
-            for j in nzrange(S, i)
-                row = S.rowval[j]
-                sval = S.nzval[j]
+            for j in nzrange(T, i)
+                row = T.rowval[j]
+                sval = T.nzval[j]
                 if splitting[row] == C_NODE && row != i
                     Bj[nnz] = row
                     if sval < 0
@@ -175,10 +177,15 @@ function rs_direct_interpolation_pass1(S, A, splitting)
         m[i] = sum
         sum += splitting[i]
     end
-    for i = 1:Bp[n]
-        Bj[i] == 0 && continue
-        Bj[i] = m[Bj[i]]
-    end
+    #@show m
+    #@show Bj
+    #l = issymmetric(S)? Bp[n]: Bp[n] + 1
+    #@show l
+    #for i = 1:l
+        #Bj[i] == 0 && continue
+    #    Bj[i] = m[Bj[i]]
+    #end
+    Bj .= m[Bj]
 
     #=Ap = A.colptr
     Aj = A.rowval

src/multilevel.jl

@@ -4,8 +4,9 @@ struct Level{Ti,Tv}
     R::SparseMatrixCSC{Ti,Tv}
 end
 
-struct MultiLevel{L, S, Pre, Post}
+struct MultiLevel{L, S, Pre, Post, Ti, Tv}
     levels::Vector{L}
+    final_A::SparseMatrixCSC{Ti,Tv}
     coarse_solver::S
     presmoother::Pre
     postsmoother::Post
@@ -15,25 +16,28 @@ abstract type CoarseSolver end
 struct Pinv <: CoarseSolver
 end
 
-MultiLevel(l::Vector{Level}, presmoother, postsmoother; coarse_solver = Pinv()) =
-    MultiLevel(l, coarse_solver, presmoother, postsmoother)
+MultiLevel(l::Vector{Level}, A, presmoother, postsmoother; coarse_solver = Pinv()) =
+    MultiLevel(l, A, coarse_solver, presmoother, postsmoother)
+Base.length(ml) = length(ml.levels) + 1
 
 function Base.show(io::IO, ml::MultiLevel)
-    op = operator_complexity(ml.levels)
-    g = grid_complexity(ml.levels)
+    op = operator_complexity(ml)
+    g = grid_complexity(ml)
     c = ml.coarse_solver
-    total_nnz = sum(nnz(level.A) for level in ml.levels)
+    total_nnz = sum(nnz(level.A) for level in ml.levels) + nnz(ml.final_A)
     lstr = ""
     for (i, level) in enumerate(ml.levels)
         lstr = lstr *
             @sprintf "   %2d   %10d   %10d [%5.2f%%]\n" i size(level.A, 1) nnz(level.A) (100 * nnz(level.A) / total_nnz)
     end
+    lstr = lstr *
+        @sprintf "   %2d   %10d   %10d [%5.2f%%]" length(ml.levels) + 1 size(ml.final_A, 1) nnz(ml.final_A) (100 * nnz(ml.final_A) / total_nnz)
     str = """
     Multilevel Solver
     -----------------
-    Operator Complexity: $op
-    Grid Complexity: $g
-    No. of Levels: $(size(ml.levels, 1))
+    Operator Complexity: $(round(op, 3))
+    Grid Complexity: $(round(g, 3))
+    No. of Levels: $(length(ml))
     Coarse Solver: $c
     Level     Unknowns     NonZeros
     -----     --------     --------
@@ -42,12 +46,12 @@ function Base.show(io::IO, ml::MultiLevel)
     print(io, str)
 end
 
-function operator_complexity(ml::Vector{Level})
-    sum(nnz(level.A) for level in ml) / nnz(ml[1].A)
+function operator_complexity(ml::MultiLevel)
+    (sum(nnz(level.A) for level in ml.levels) + nnz(ml.final_A)) / nnz(ml.levels[1].A)
 end
 
-function grid_complexity(ml::Vector{Level})
-    sum(size(level.A, 1) for level in ml) / size(ml[1].A, 1)
+function grid_complexity(ml::MultiLevel)
+    (sum(size(level.A, 1) for level in ml.levels) + size(ml.final_A, 1)) / size(ml.levels[1].A, 1)
 end
 
 abstract type Cycle end
@@ -83,8 +87,8 @@ function __solve{T}(v::V, ml, x::Vector{T}, b::Vector{T}, lvl)
     coarse_b = ml.levels[lvl].R * res
     coarse_x = zeros(T, size(coarse_b))
 
-    if lvl == length(ml.levels) - 1
-        coarse_x = coarse_solver(ml.coarse_solver, ml.levels[end].A, coarse_b)
+    if lvl == length(ml.levels)
+        coarse_x = coarse_solver(ml.coarse_solver, ml.final_A, coarse_b)
     else
         coarse_x = __solve(v, ml, coarse_x, coarse_b, lvl + 1)
     end

test/randlap.jld

@@ -2,6 +2,8 @@ using AMG
 using Base.Test
 using JLD
 
+@testset "AMG Tests" begin
+
 graph = load("test.jld")["G"]
 ref_S = load("ref_S_test.jld")["G"]
 ref_split = readdlm("ref_split_test.txt")
@@ -27,7 +29,6 @@ end
 # Ruge-Stuben splitting
 S = poisson(7)
 @test split_nodes(RS(), S) == [0, 1, 0, 1, 0, 1, 0]
-@show "buzz"
 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]
@@ -54,7 +55,9 @@ P, R = AMG.direct_interpolation(A, copy(A), splitting)
               0.0  1.0  0.0
               0.0  0.5  0.5
               0.0  0.0  1.0 ]
-
+A = load("thing.jld")["G"]
+ml = ruge_stuben(A)
+@test size(ml.levels[2].A, 1) == 19
 end
 
 @testset "Coarse Solver" begin
@@ -67,13 +70,31 @@ end
 A = poisson(1000)
 A = float.(A) #FIXME
 ml = AMG.ruge_stuben(A)
-@test length(ml.levels) == 8
-s = [1000, 500, 250, 125, 62, 31, 15, 7]
-n = [2998, 1498, 748, 373, 184, 91, 43, 19]
-for i = 1:8
+@test length(ml) == 8
+s = [1000, 500, 250, 125, 62, 31, 15]
+n = [2998, 1498, 748, 373, 184, 91, 43]
+for i = 1:7
     @test size(ml.levels[i].A, 1) == s[i]
     @test nnz(ml.levels[i].A) == n[i]
 end
+@test size(ml.final_A, 1) == 7
+@test nnz(ml.final_A) == 19
+
+A = load("randlap.jld")["G"]
+ml = ruge_stuben(A)
+@test length(ml) == 3
+s = [100, 17]
+n = [2066, 289]
+for i = 1:2
+    @test size(ml.levels[i].A, 1) == s[i]
+    @test nnz(ml.levels[i].A) == n[i]
+end
+@test size(ml.final_A, 1) == 2
+@test nnz(ml.final_A) == 4
+@test round(AMG.operator_complexity(ml), 3) ≈ 1.142
+@test round(AMG.grid_complexity(ml), 3) ≈ 1.190
+
+
 end
 
 @testset "Solver" begin
@@ -82,4 +103,12 @@ A = float.(A)
 ml = ruge_stuben(A)
 x = solve(ml, A * ones(1000))
 @test sum(abs2, x - ones(1000)) < 1e-10
+
+A = load("randlap.jld")["G"]
+ml = ruge_stuben(A)
+x = solve(ml, A * ones(100))
+@test sum(abs2, x - zeros(100)) < 1e-10
+
+end
+
 end