Merge pull request #2 from ranjanan/levels

Construct multilevel solver

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
+    MultiLevel(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)
@@ -65,7 +67,8 @@ function rs_direct_interpolation_pass1(S, A, splitting)
          if splitting[i] == C_NODE
              nnz += 1
          else
-            for jj = Sp[i]:Sp[i+1]-1
+            for jj = Sp[i]:Sp[i+1]
+                jj > length(Sj) && continue
                 if splitting[Sj[jj]] == C_NODE && Sj[jj] != i
                     nnz += 1
                 end
@@ -92,7 +95,6 @@ function rs_direct_interpolation_pass1(S, A, splitting)
     Bj = zeros(Ti, Bp[end])
     Bx = zeros(Float64, Bp[end])
     n_nodes = size(A, 1)
-    #Bp += 1
 
     for i = 1:n_nodes
         if splitting[i] == C_NODE
@@ -101,7 +103,8 @@ function rs_direct_interpolation_pass1(S, A, splitting)
         else
             sum_strong_pos = 0
             sum_strong_neg = 0
-            for jj = Sp[i]: Sp[i+1]-1
+            for jj = Sp[i]: Sp[i+1]
+                jj > length(Sj) && 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(Aj) && continue
                 if Aj[jj] == i
                     diag += Ax[jj]
                 else
@@ -134,19 +138,20 @@ function rs_direct_interpolation_pass1(S, A, splitting)
                 beta = 0
             end
 
-            neg_coeff = -alpha/diag;
-            pos_coeff = -beta/diag;
+            neg_coeff = -alpha / diag
+            pos_coeff = -beta / diag
 
             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
                         Bx[nnz] = neg_coeff * Sx[jj]
                     else
                         Bx[nnz] = pos_coeff * Sx[jj]
-                        nnz += 1
                     end
+                    nnz += 1
                 end
             end
         end
@@ -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,48 @@
 struct Level{Ti,Tv}
     A::SparseMatrixCSC{Ti,Tv}
+    P::SparseMatrixCSC{Ti,Tv}
+    R::SparseMatrixCSC{Ti,Tv}
+end
+
+struct MultiLevel{L, S}
+    levels::Vector{L}
+    coarse_solver::S
+end
+
+abstract type CoarseSolver end
+struct Pinv <: CoarseSolver
+end
+MultiLevel(l::Vector{Level}; coarse_solver = Pinv()) =
+    MultiLevel(l, coarse_solver)
+
+function Base.show(io::IO, ml::MultiLevel)
+    op = operator_complexity(ml.levels)
+    g = grid_complexity(ml.levels)
+    c = ml.coarse_solver
+    total_nnz = sum(nnz(level.A) for level in ml.levels)
+    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
+    str = """
+    Multilevel Solver
+    -----------------
+    Operator Complexity: $op
+    Grid Complexity: $g
+    No. of Levels: $(size(ml.levels, 1))
+    Coarse Solver: $c
+    Level     Unknowns     NonZeros
+    -----     --------     --------
+    $lstr
+    """
+    print(io, str)
+end
+
+function operator_complexity(ml::Vector{Level})
+    sum(nnz(level.A) for level in ml) / nnz(ml[1].A)
+end
+
+function grid_complexity(ml::Vector{Level})
+    sum(size(level.A, 1) for level in ml) / size(ml[1].A, 1)
 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,19 +27,17 @@ 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
 
 function find_max_off_diag(neighbors, col)
-    max_offdiag = 0
-    for (i,v) in enumerate(neighbors)
-        if col != i
-            max_offdiag = max(max_offdiag, abs(v))
-        end
+    maxval = zero(eltype(neighbors))
+    for i in 1:length(neighbors.nzval)
+        maxval = max(maxval, ifelse(neighbors.nzind[i] == col, 0, abs(neighbors.nzval[i])))
     end
-    max_offdiag
+    return maxval
 end
 
 function scale_cols_by_largest_entry(A::SparseMatrixCSC)
@@ -51,16 +50,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