Merge pull request #22 from ranjanan/satest

Tests for Smooth Aggregation

Author: ranjanan

src/aggregation.jl

@@ -57,8 +57,7 @@ function extend_hierarchy!(levels, strength, aggregate, smooth,
     S = strength_of_connection(strength, A, bsr_flag)
 
     # Aggregation operator
-    AggOp = aggregation(aggregate, S')
-
+    AggOp = aggregation(aggregate, S)
     # b = zeros(eltype(A), size(A, 1))
 
     # Improve candidates
@@ -81,51 +80,24 @@ 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_fine, n_coarse = size(A)
     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)
-
+    # copy!(A.nzval, B)
+    for i = 1:n_col
+        for j in nzrange(A,i)
+            row = A.rowval[j]
+            A.nzval[j] = B[row]
+        end
+    end
+    k = 1
     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
@@ -136,82 +108,22 @@ function fit_candidates(AggOp, B, tol = 1e-10)
         end
         for j in nzrange(A, i)
             row = A.rowval[j]
-            Qx[row] *= scale
+            # Qx[row] *= scale
+            #@show k
+            # Qx[k] *= scale
+            # k += 1
+            A.nzval[j] *= 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
+    # SparseMatrixCSC(size(A)..., A.colptr, A.rowval, Qx), R
+    A, R
 end
 function norm_col(A, Qx, i)
     s = zero(eltype(A))
     for j in nzrange(A, i)
         val = Qx[A.rowval[j]]
+        # val = A.nzval[A.rowval[j]]
         s += val*val
     end
     sqrt(s)

src/gallery.jl

@@ -1 +1,3 @@
-poisson(n) = sparse(Tridiagonal(fill(-1, n-1), fill(2, n), fill(-1, n-1)))
+poisson(T, n) = sparse(Tridiagonal(fill(T(-1), n-1), 
+                        fill(T(2), n), fill(T(-1), n-1)))
+poisson(n) = poisson(Float64, n)

src/multilevel.jl

@@ -88,6 +88,7 @@ function solve{T}(ml::MultiLevel, b::Vector{T},
         push!(residuals, T(norm(b - A * x)))
     end
 
+    # @show residuals
     if log
         return x, residuals
     else

src/smoother.jl

@@ -18,8 +18,11 @@ presmoother!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
 postsmoother!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
 relax!(s, A, x, b) = smoother!(s, s.sweep, A, x, b)
 
-smoother!(s::GaussSeidel, ::ForwardSweep, A, x, b) =
-                    gs!(A, b, x, 1, 1, size(A, 1))
+function smoother!(s::GaussSeidel, ::ForwardSweep, A, x, b)
+    for i in 1:s.iter
+        gs!(A, b, x, 1, 1, size(A, 1))
+    end
+end
 
 function smoother!(s::GaussSeidel, ::SymmetricSweep, A, x, b)
     for i in 1:s.iter
@@ -28,8 +31,11 @@ function smoother!(s::GaussSeidel, ::SymmetricSweep, A, x, b)
     end
 end
 
-smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b) =
-    gs!(A, b, x, size(A,1), -1, 1)
+function smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b)
+    for i in 1:s.iter
+        gs!(A, b, x, size(A,1), -1, 1)
+    end
+end
 
 
 function gs!(A, b, x, start, step, stop)
@@ -94,11 +100,13 @@ end
 
 struct DiagonalWeighting
 end
+struct LocalWeighting
+end
 
 function smooth_prolongator(j::JacobiProlongation,
                                         A, T, S, B,
                                         degree = 1,
-                                        weighting = DiagonalWeighting())
+                                        weighting = LocalWeighting())
     D_inv_S = weight(weighting, A, j.ω)
     P = T
     for i = 1:degree
@@ -111,7 +119,14 @@ function weight(::DiagonalWeighting, S, ω)
     D_inv = 1 ./ diag(S)
     D_inv_S = scale_rows(S, D_inv)
     (eltype(S)(ω) / approximate_spectral_radius(D_inv_S)) * D_inv_S
-    #(ω) * D_inv_S
+    # (ω) * D_inv_S
+end
+
+function weight(::LocalWeighting, S, ω)
+    D = abs.(S) * ones(size(S, 1))
+    D_inv = 1 ./ D[find(D)]
+    D_inv_S = scale_rows(S, D_inv)
+    eltype(S)(ω) * D_inv_S
 end
 
 #approximate_spectral_radius(A) =
@@ -127,4 +142,4 @@ function scale_rows!(ret, S, v)
     end
     ret
 end
-scale_rows(S, v) = scale_rows!(deepcopy(S), S,  v)
+scale_rows(S, v) = scale_rows!(deepcopy(S), S,  v)

src/strength.jl

@@ -107,20 +107,26 @@ function strength_of_connection{T}(s::SymmetricStrength{T}, A, bsr_flag = false)
         for j in nzrange(A, i)
             row = A.rowval[j]
             val = A.nzval[j]
-            if val*val < eps_Aii * diags[row]
-                S.nzval[j] = 0
+            if row != i
+                if val*val < eps_Aii * diags[row]
+                    S.nzval[j] = 0
+                end
             end
         end
     end
 
     dropzeros!(S)
 
-    S.nzval .= abs.(S.nzval)
-    # for i = 1:size(S.nzval, 1)
+    # S.nzval .= abs.(S.nzval)
+    #for i = 1:size(S.nzval, 1)
     #     S.nzval[i] = abs(S.nzval[i])
-    # end
+    #end
 
     scale_cols_by_largest_entry!(S)
+    
+    for i = 1:size(S.nzval, 1)
+         S.nzval[i] = abs(S.nzval[i])
+    end
 
     S
 end

src/utils.jl

@@ -5,7 +5,7 @@ function approximate_spectral_radius(A, tol = 0.01,
     symmetric = false
 
     # Initial guess
-    v0 = rand(eltype(A), size(A,1))
+    v0 = rand(eltype(A), size(A,2))
     maxiter = min(size(A, 1), maxiter)
     ev = zeros(eltype(A), maxiter)
     max_index = 0
@@ -21,7 +21,7 @@ function approximate_spectral_radius(A, tol = 0.01,
         # m, max_index = findmax(abs.(ev))
         m, max_index = findmaxabs(ev)
         error = H[nvecs, nvecs-1] * evect[end, max_index]
-        if abs(error) / abs(ev[max_index]) < tol || flag
+        if (abs(error) / abs(ev[max_index]) < tol) || flag
             # v0 = X * evect[:, max_index]
             A_mul_B!(v0, X, evect[:, max_index])
             break
@@ -63,6 +63,7 @@ function approximate_eigenvalues(A, tol, maxiter, symmetric, v0)
     # V = Vector{Vector{eltype(A)}}(maxiter + 1)
     # V = zeros(size(A,1), maxiter)
     # V[1] = v0
+    breakdown = find_breakdown(eltype(A))
     flag = false
 
     for j = 1:maxiter
@@ -73,11 +74,11 @@ function approximate_eigenvalues(A, tol, maxiter, symmetric, v0)
         for (i,v) in enumerate(V)
         # for i = 1:j
             v = V[i]
-            H[i,j] = dot(v, w)
+            H[i,j] = dot(conj.(v), w)
             BLAS.axpy!(-H[i,j], v, w)
         end
         H[j+1,j] = norm(w)
-        if H[j+1, j] < eps()
+        if H[j+1, j] < breakdown
             flag = true
             if H[j+1,j] != 0
                 scale!(w, 1/H[j+1,j])
@@ -94,3 +95,6 @@ function approximate_eigenvalues(A, tol, maxiter, symmetric, v0)
 
     Vects, Eigs, H, V, flag
 end
+
+find_breakdown(::Type{Float64}) = eps(Float64) * 10^6
+find_breakdown(::Type{Float32}) = eps(Float64) * 10^3

test/ref_R.jld

@@ -4,6 +4,8 @@ using JLD
 using IterativeSolvers
 import AMG: V, coarse_solver, Pinv, Classical
 
+include("sa_tests.jl")
+
 @testset "AMG Tests" begin
 
 graph = load("test.jld")["G"]
@@ -234,4 +236,33 @@ end
 
 end
 
+# Smoothed Aggregation
+@testset "Smoothed Aggregation" begin
+
+@testset "Symmetric Strength of Connection" begin
+
+        test_symmetric_soc()
+end
+
+@testset "Standard Aggregation" begin
+
+        test_standard_aggregation()
+end
+
+@testset "Fit Candidates" begin
+        test_fit_candidates()
+end
+
+@testset "Approximate Spectral Radius" begin
+        test_approximate_spectral_radius()
+end
+end
+
+@testset "Gauss Seidel" begin
+        test_gauss_seidel()
+end
+
+@testset "Jacobi Prolongation" begin
+        test_jacobi_prolongator()
+end
 end