Get level construction working

Author: ranjanan

src/AMG.jl

@@ -9,6 +9,10 @@ export split_nodes, RS
 include("gallery.jl")
 export poisson
 
+include("smoother.jl")
+
+include("multilevel.jl")
+
 include("classical.jl")
 
 end # module

src/classical.jl

@@ -7,36 +7,38 @@ struct Solver{S,T,P,PS}
     max_coarse::Int64
 end
 
-struct Level{T}
-    A::T
-end
-
 function ruge_stuben(A::SparseMatrixCSC;
-                strength = Classical(),
+                strength = Classical(0.25),
                 CF = RS(),
                 presmoother = GaussSiedel(),
                 postsmoother = GaussSiedel(),
                 max_levels = 10,
                 max_coarse = 500)
 
-        s = Solver(strength, CF, presmoother,
-                    postsmoother, max_levels, max_levels)
+    s = Solver(strength, CF, presmoother,
+                postsmoother, max_levels, max_levels)
 
-        levels = [Level(A)]
+    levels = Vector{Level}()
 
-        while length(levels) < max_levels && size(levels[end].A, 1)
-            extend_heirarchy!(levels, strength, CF, A)
+    while length(levels) < max_levels
+        A = extend_heirarchy!(levels, strength, CF, A)
+        if size(levels[end].A, 1) < max_coarse
+            break
         end
+    end
+    levels
 end
 
 function extend_heirarchy!(levels::Vector{Level}, strength, CF, A)
     S = strength_of_connection(strength, A)
     splitting = split_nodes(CF, S)
     P, R = direct_interpolation(A, S, splitting)
+    push!(levels, Level(A, P, R))
+    A = R * A * P
 end
 
 function direct_interpolation{T,V}(A::T, S::T, splitting::Vector{V})
-    
+
     fill!(S.nzval, 1.)
     S = A .* S
     Pp = rs_direct_interpolation_pass1(S, A, splitting)
@@ -102,6 +104,7 @@ function rs_direct_interpolation_pass1(S, A, splitting)
             sum_strong_pos = 0
             sum_strong_neg = 0
             for jj = Sp[i]: Sp[i+1]-1
+                jj > length(Sp) && continue
                 if splitting[Sj[jj]] == C_NODE && Sj[jj] != i
                     if Sx[jj] < 0
                         sum_strong_neg += Sx[jj]
@@ -115,6 +118,7 @@ function rs_direct_interpolation_pass1(S, A, splitting)
             sum_all_neg = 0
             diag = 0;
             for jj = Ap[i]:Ap[i+1]
+                jj > length(Ap) && continue
                 if Aj[jj] == i
                     diag += Ax[jj]
                 else
@@ -139,6 +143,7 @@ function rs_direct_interpolation_pass1(S, A, splitting)
 
             nnz = Bp[i]
             for jj = Sp[i]:Sp[i+1]
+                jj > length(Sj) && continue
                 if splitting[Sj[jj]] == C_NODE && Sj[jj] != i
                     Bj[nnz] = Sj[jj]
                     if Sx[jj] < 0
@@ -159,6 +164,7 @@ function rs_direct_interpolation_pass1(S, A, splitting)
        sum += splitting[i]
    end
    for i = 1:Bp[n_nodes]
+       Bj[i] == 0 && continue
        Bj[i] = m[Bj[i]]
    end
 

src/multilevel.jl

@@ -1,3 +1,5 @@
 struct Level{Ti,Tv}
     A::SparseMatrixCSC{Ti,Tv}
+    P::SparseMatrixCSC{Ti,Tv}
+    R::SparseMatrixCSC{Ti,Tv}
 end

src/smoother.jl

@@ -0,0 +1,3 @@
+abstract type Smoother end
+struct GaussSiedel <: Smoother
+end

src/strength.jl

@@ -2,6 +2,7 @@ abstract type Strength end
 struct Classical{T} <: Strength
     θ::T
 end
+Classical(;θ = 0.25) = Classical(θ)
 
 function strength_of_connection{T}(c::Classical{T}, A::SparseMatrixCSC)
 
@@ -14,8 +15,8 @@ function strength_of_connection{T}(c::Classical{T}, A::SparseMatrixCSC)
 
     for i = 1:n
         neighbors = A[:,i]
-        m = find_max_off_diag(neighbors, i)
-        threshold = θ * m
+        _m = find_max_off_diag(neighbors, i)
+        threshold = θ * _m
         for j in nzrange(A, i)
             row = A.rowval[j]
             val = A.nzval[j]
@@ -26,7 +27,7 @@ function strength_of_connection{T}(c::Classical{T}, A::SparseMatrixCSC)
             end
         end
     end
-    S = sparse(I, J, V)
+    S = sparse(I, J, V, m, n)
 
     scale_cols_by_largest_entry(S)
 end
@@ -51,16 +52,16 @@ function scale_cols_by_largest_entry(A::SparseMatrixCSC)
 
     k = 1
     for i = 1:n
-        m = maximum(A[:,i])
+        _m = maximum(A[:,i])
         for j in nzrange(A, i)
             row = A.rowval[j]
             val = A.nzval[j]
             I[k] = row
             J[k] = i
-            V[k] = val / m
+            V[k] = val / _m
             k += 1
         end
     end
 
-    sparse(I,J,V)
+    sparse(I,J,V,m,n)
 end