Merge pull request #15 from ranjanan/standard_strength

Smoothed Aggregation

Author: ranjanan

README.md

@@ -34,7 +34,15 @@ ml = ruge_stuben(A) # Construct a Ruge-Stuben solver
 solve(ml, A * ones(1000)) # should return ones(1000)
 ```
 
-## Roadmap
+### As a Preconditioner
+You can use AMG as a preconditioner for Krylov methods such as Conjugate Gradients.
+```julia
+import IterativeSolvers: cg
+p = aspreconditioner(ml)
+c = cg(A, A*ones(1000), Pl = p) 
+```
+
+## Features and Roadmap
 
 This package currently supports: 
 1. Ruge-Stuben Solver
@@ -43,9 +51,14 @@ This package currently supports:
 4. Gauss-Siedel smoothers
 5. V cycle multigrid
 
-In the future, this package will support
+The following have experimental support:
 1. SmoothedAggregation Solver
 2. Standard Strength of Conneciton
-3. Other splitting methods (like CLJP)
-4. SOR, Jacobi smoothers
-5. W, F, AMLI cycles
+
+In the future, this package will support:
+1. Other splitting methods (like CLJP)
+2. SOR, Jacobi smoothers
+3. W, F, AMLI cycles
+
+#### Acknowledgements
+This package has been heavily inspired by the [`PyAMG`](http://github.com/pyamg/pyamg) project. 

src/AMG.jl

@@ -1,9 +1,13 @@
 module AMG
 
 import IterativeSolvers: gauss_seidel!
+using Base.Threads
+
+include("utils.jl")
+export approximate_spectral_radius
 
 include("strength.jl")
-export strength_of_connection, Classical
+export strength_of_connection, Classical, SymmetricStrength
 
 include("splitting.jl")
 export split_nodes, RS
@@ -12,14 +16,21 @@ include("gallery.jl")
 export poisson
 
 include("smoother.jl")
-export GaussSeidel, SymmetricSweep, ForwardSweep, BackwardSweep
+export GaussSeidel, SymmetricSweep, ForwardSweep, BackwardSweep,
+        smooth_prolongator, JacobiProlongation
 
 include("multilevel.jl")
 export solve
 
 include("classical.jl")
 export ruge_stuben
 
+include("aggregate.jl")
+export aggregation, StandardAggregation
+
+include("aggregation.jl")
+export fit_candidates, smoothed_aggregation
+
 include("preconditioner.jl")
 export aspreconditioner
 

src/aggregate.jl

@@ -0,0 +1,101 @@
+struct StandardAggregation
+end
+
+function aggregation(::StandardAggregation, S)
+
+    n = size(S, 1)
+    x = zeros(Int, n)
+    y = zeros(Int, n)
+
+    next_aggregate = 1
+
+    for i = 1:n
+        if x[i] != 0
+            continue
+        end
+
+        has_agg_neighbors = false
+        has_neighbors = false
+
+        for j in nzrange(S, i)
+            row = S.rowval[j]
+            if row != i
+                has_neighbors = true
+                if x[row] != 0
+                    has_agg_neighbors = true
+                    break
+                end
+            end
+        end
+
+        if !has_neighbors
+            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
+            end
+
+            next_aggregate += 1
+        end
+    end
+
+
+    # Pass 2
+    for i = 1:n
+        if x[i] != 0
+            continue
+        end
+
+        for j in nzrange(S, i)
+            row = S.rowval[j]
+            x_row = x[row]
+            if x_row > 0
+                x[i] = -x_row
+                break
+            end
+        end
+    end
+
+    next_aggregate -= 1
+
+    # Pass 3
+    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
+        end
+
+        x[i] = next_aggregate
+        y[next_aggregate + 1] = i
+
+        for j in nzrange(S, i)
+            row = S.rowval[j]
+
+            if x[row] == 0
+                x[row] = next_aggregate
+            end
+        end
+
+        next_aggregate += 1
+    end
+
+    y = y[1:next_aggregate]
+    M,N = (n, next_aggregate)
+
+    Tp = collect(1:n+1)
+    x .= x .+ 1
+    Tx = ones(length(x))
+
+    SparseMatrixCSC(N, M, Tp, x, Tx)
+end

src/aggregation.jl

@@ -0,0 +1,218 @@
+function smoothed_aggregation{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti},
+                        symmetry = HermitianSymmetry(),
+                        strength = SymmetricStrength(),
+                        aggregate = StandardAggregation(),
+                        smooth = JacobiProlongation(4.0/3.0),
+                        presmoother = GaussSeidel(),
+                        postsmoother = GaussSeidel(),
+                        improve_candidates = GaussSeidel(4),
+                        max_levels = 10,
+                        max_coarse = 10,
+                        diagonal_dominance = false,
+                        keep = false,
+                        coarse_solver = Pinv())
+
+
+    n = size(A, 1)
+    # B = kron(ones(n, 1), eye(1))
+    B = ones(n)
+
+    #=max_levels, max_coarse, strength =
+        levelize_strength_or_aggregation(max_levels, max_coarse, strength)
+    max_levels, max_coarse, aggregate =
+        levelize_strength_or_aggregation(max_levels, max_coarse, aggregate)
+
+    improve_candidates =
+        levelize_smooth_or_improve_candidates(improve_candidates, max_levels)=#
+    # str = [stength for _ in 1:max_levels - 1]
+    # agg = [aggregate for _ in 1:max_levels - 1]
+    # sm = [smooth for _ in 1:max_levels]
+
+    levels = Vector{Level{Tv,Ti}}()
+    bsr_flag = false
+
+    while length(levels) < max_levels # && size(A, 1) > max_coarse
+        A, B, bsr_flag = extend_hierarchy!(levels, strength, aggregate, smooth,
+                                improve_candidates, diagonal_dominance,
+                                keep, A, B, symmetry, bsr_flag)
+        if size(A, 1) <= max_coarse
+            break
+        end
+    end
+    #=A, B = extend_hierarchy!(levels, strength, aggregate, smooth,
+                            improve_candidates, diagonal_dominance,
+                            keep, A, B, symmetry)=#
+    MultiLevel(levels, A, presmoother, postsmoother)
+end
+
+struct HermitianSymmetry
+end
+
+function extend_hierarchy!(levels, strength, aggregate, smooth,
+                            improve_candidates, diagonal_dominance, keep,
+                            A, B,
+                            symmetry, bsr_flag)
+
+    # Calculate strength of connection matrix
+    S = strength_of_connection(strength, A, bsr_flag)
+
+    # Aggregation operator
+    AggOp = aggregation(aggregate, S')
+
+    # b = zeros(eltype(A), size(A, 1))
+
+    # Improve candidates
+    # relax!(improve_candidates, A, B, b)
+    T, B = fit_candidates(AggOp, B)
+
+    P = smooth_prolongator(smooth, A, T, S, B)
+    R = construct_R(symmetry, P)
+    push!(levels, Level(A, P, R))
+
+    A = R * A * P
+
+    dropzeros!(A)
+
+    bsr_flag = true
+
+    A, B, bsr_flag
+end
+construct_R(::HermitianSymmetry, P) = P'
+
+function fit_candidates(AggOp, B, tol = 1e-10)
+
+    # K1 = Int(size(B, 1) / N_fine)
+    # K2 = size(B, 2)
+
+    A = AggOp.'
+
+    n_coarse = size(A, 2)
+    n_fine = size(A, 1)
+    n_col = n_coarse
+
+    # R = zeros(eltype(B), N_coarse, K2, K2)
+    R = zeros(eltype(B), n_coarse)
+    # Qx = zeros(eltype(B), nnz(AggOp), K1, K2)
+    Qx = zeros(eltype(B), nnz(A))
+
+
+    # R = vec(R)
+    # Qx = vec(Qx)
+
+    # n_row = N_fine
+    # n_col = N_coarse
+
+    # BS = K1 * K2
+
+    #=for i = 1:n_col
+        Ax_start = 1 + BS * A.colptr[i]
+
+        for j in nzrange(A, i)
+            B_start = 1 + BS * A.rowval[j]
+            B_end = B_start + BS
+            @show B_start
+            @show B_end
+            for ind in B_start:B_end
+                A.nzval[ind + Ax_start] = B[ind]
+            end
+            Ax_start += BS
+        end
+    end=#
+    # copy!(A.nzval, B)
+    copy!(Qx, B)
+    # @show size(A.nzval)
+    # @show size(B)
+
+    for i = 1:n_col
+        norm_i = norm_col(A, Qx, i)
+        # norm_i = norm(A[:,i])
+        threshold_i = tol * norm_i
+        if norm_i > threshold_i
+            scale = 1 / norm_i
+            R[i] = norm_i
+        else
+            scale = 0
+            R[i] = 0
+        end
+        for j in nzrange(A, i)
+            row = A.rowval[j]
+            Qx[row] *= scale
+        end
+        #=col_start = A.colptr[i]
+        col_end = A.colptr[i+1]
+
+        Ax_start = 1 + BS * col_start
+        Ax_end = 1 + BS * col_end
+        R_start = 1 + i * K2 * K2
+
+        for bj = 1:K2
+            norm_j = zero(eltype(A))
+            Ax_col = Ax_start + bj
+            while Ax_col < Ax_end
+                norm_j += norm(A.nzval[Ax_col])
+                Ax_col += K2
+            end
+            norm_j = sqrt(norm_j)
+
+            threshold_j = tol * norm_j
+
+            for bi = 1:bj
+                dot_prod = zero(eltype(A))
+
+                Ax_bi = Ax_start + bj
+                Ax_bj = Ax_start + bj
+                while Ax_bi < Ax_end
+                    dot_prod += dot(A.nzval[Ax_bj], A.nzval[Ax_bi])
+                    Ax_bi    += K2
+                    Ax_bj    += K2
+                end
+
+                Ax_bi = Ax_start + bi;
+                Ax_bj = Ax_start + bj;
+                while Ax_bi < Ax_end
+                    A.nzval[Ax_bj] -= dot_prod * A.nzval[Ax_bi]
+                    Ax_bi  += K2
+                    Ax_bj  += K2
+                end
+
+                R[R_start + K2 * bi + bj] = dot_prod
+            end
+
+            norm_j = zero(eltype(A))
+            Ax_bj = Ax_start + bj
+            while Ax_bj < Ax_end
+                norm_j += norm(A.nzval[Ax_bj])
+                Ax_bj  += K2
+                norm_j = sqrt(norm_j)
+            end
+
+            if norm_j > threshold_j
+                scale = 1 / norm_j
+                R[R_start + K2 * bj + bj] = norm_j
+            else
+                scale = zero(eltype(A))
+                R[R_start + K2 * bj + bj] = 0
+            end
+
+            Ax_bj = Ax_start + bj
+            while Ax_bj < Ax_end
+                A.nzval[Ax_bj] *= scale
+                Ax_bj  += K2
+            end
+        end=#
+    end
+
+    #Q = SparseMatrixCSC(N_coarse, N_fine,
+                        #.colptr, A.rowval, Qx)
+
+    #R = reshape(R, N_coarse, K2)
+    SparseMatrixCSC(size(A)..., A.colptr, A.rowval, Qx), R
+end
+function norm_col(A, Qx, i)
+    s = zero(eltype(A))
+    for j in nzrange(A, i)
+        val = Qx[A.rowval[j]]
+        s += val*val
+    end
+    sqrt(s)
+end