Generalize direct interpol for non-symmetric matrices and edge cases

mohamed82008 · mohamed82008 · commit 8bc5debee183 · 2018-08-02T22:41:11.000+02:00
diff --git a/src/classical.jl b/src/classical.jl
@@ -22,28 +22,26 @@ function ruge_stuben(A::SparseMatrixCSC{Ti,Tv};
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
         A = extend_heirarchy!(levels, strength, CF, A)
-        #if size(A, 1) <= max_coarse
-        #    break
-        #end
     end
     MultiLevel(levels, A, presmoother, postsmoother)
 end
 
 function extend_heirarchy!(levels::Vector{Level{Ti,Tv}}, strength, CF, A::SparseMatrixCSC{Ti,Tv}) where {Ti,Tv}
-    S, T = strength_of_connection(strength, A)
+    At = copy(A')
+    S, T = strength_of_connection(strength, At)
     splitting = split_nodes(CF, S)
-    P, R = direct_interpolation(A, T, splitting)
+    P, R = direct_interpolation(At, T, splitting)
     push!(levels, Level(A, P, R))
     A = R * A * P
 end
 
-function direct_interpolation(A, T, splitting)
-    fill!(T.nzval, eltype(A)(1))
-    T .= A .* T
+function direct_interpolation(At, T, splitting)
+    fill!(T.nzval, eltype(At)(1))
+    T .= At .* T
+    
     Pp = rs_direct_interpolation_pass1(T, splitting)
-    Px, Pj, Pp = rs_direct_interpolation_pass2(A, T, splitting, Pp)
-
-    R = SparseMatrixCSC(maximum(Pj), size(A, 1), Pp, Pj, Px)
+    Px, Pj, Pp = rs_direct_interpolation_pass2(At, T, splitting, Pp)
+    R = SparseMatrixCSC(maximum(Pj), size(At, 1), Pp, Pj, Px)
     P = copy(R')
 
     P, R
@@ -60,7 +58,7 @@ function rs_direct_interpolation_pass1(T, splitting)
         else
             for j in nzrange(T, i)
                 row = T.rowval[j]
-                if splitting[row] == C_NODE && row != i
+                if splitting[row] == C_NODE
                     nnzplus1 += 1
                 end
             end
@@ -71,15 +69,15 @@ function rs_direct_interpolation_pass1(T, splitting)
  end
 
 
-function rs_direct_interpolation_pass2(A::SparseMatrixCSC{Tv,Ti},
+function rs_direct_interpolation_pass2(At::SparseMatrixCSC{Tv,Ti},
                                                 T::SparseMatrixCSC{Tv, Ti},
                                                 splitting::Vector{Ti},
                                                 Bp::Vector{Ti}) where {Tv,Ti}
                                                 
     Bx = zeros(Tv, Bp[end] - 1)
     Bj = zeros(Ti, Bp[end] - 1)
 
-    n = size(A, 1)
+    n = size(At, 1)
 
     for i = 1:n
         if splitting[i] == C_NODE
@@ -91,7 +89,7 @@ function rs_direct_interpolation_pass2(A::SparseMatrixCSC{Tv,Ti},
             for j in nzrange(T, i)
                 row = T.rowval[j]
                 sval = T.nzval[j]
-                if splitting[row] == C_NODE && row != i
+                if splitting[row] == C_NODE
                     if sval < 0
                         sum_strong_neg += sval
                     else
@@ -102,9 +100,9 @@ function rs_direct_interpolation_pass2(A::SparseMatrixCSC{Tv,Ti},
             sum_all_pos = zero(Tv)
             sum_all_neg = zero(Tv)
             diag = zero(Tv)
-            for j in nzrange(A, i)
-                row = A.rowval[j]
-                aval = A.nzval[j]
+            for j in nzrange(At, i)
+                row = At.rowval[j]
+                aval = At.nzval[j]
                 if row == i
                     diag += aval
                 else
@@ -115,28 +113,43 @@ function rs_direct_interpolation_pass2(A::SparseMatrixCSC{Tv,Ti},
                     end
                 end
             end
-            alpha = sum_all_neg / sum_strong_neg
-            beta = sum_all_pos / sum_strong_pos
 
             if sum_strong_pos == 0
-                diag += sum_all_pos
-                beta = zero(beta)
+                beta = zero(diag)
+                if diag >= 0
+                    diag += sum_all_pos
+                end
+            else
+                beta = sum_all_pos / sum_strong_pos
             end
 
-            neg_coeff = -1 * alpha / diag
-            pos_coeff = -1 * beta / diag
+            if sum_strong_neg == 0
+                alpha = zero(diag)
+                if diag < 0
+                    diag += sum_all_neg
+                end
+            else
+                alpha = sum_all_neg / sum_strong_neg
+            end
 
-            nnz = Bp[i]
+            if isapprox(diag, 0, atol=eps(Tv))
+                neg_coeff = Tv(0)
+                pos_coeff = Tv(0)
+            else
+                neg_coeff = alpha / diag
+                pos_coeff = beta / diag
+            end
 
+            nnz = Bp[i]
             for j in nzrange(T, i)
                 row = T.rowval[j]
                 sval = T.nzval[j]
-                if splitting[row] == C_NODE && row != i
+                if splitting[row] == C_NODE
                     Bj[nnz] = row
                     if sval < 0
-                        Bx[nnz] = neg_coeff * sval
+                        Bx[nnz] = abs(neg_coeff * sval)
                     else
-                        Bx[nnz] = pos_coeff * sval
+                        Bx[nnz] = abs(pos_coeff * sval)
                     end
                     nnz += 1
                 end
diff --git a/src/multilevel.jl b/src/multilevel.jl
@@ -74,7 +74,7 @@ struct V <: Cycle
 end
 
 """
-    solve(ml::MultiLevel, b::Vector, cycle, kwargs...)
+    solve(ml::MultiLevel, b::AbstractVector, cycle, kwargs...)
 
 Execute multigrid cycling.
 
@@ -92,7 +92,7 @@ Keyword Arguments
 * log::Bool - return vector of residuals along with solution
 
 """
-function solve(ml::MultiLevel, b::Vector{T},
+function solve(ml::MultiLevel, b::AbstractVector{T},
                                     cycle::Cycle = V();
                                     maxiter::Int = 100,
                                     tol::Float64 = 1e-5,
diff --git a/src/preconditioner.jl b/src/preconditioner.jl
@@ -1,3 +1,5 @@
+import Compat.LinearAlgebra: \, *, ldiv!, mul!
+
 struct Preconditioner
     ml::MultiLevel
 end
@@ -10,6 +12,7 @@ aspreconditioner(ml::MultiLevel) = Preconditioner(ml)
     A_mul_B!(b, p::Preconditioner, x) = A_mul_B!(b, p.ml.levels[1].A, x)
 else
     import Compat.LinearAlgebra: \, *, ldiv!, mul!
+    ldiv!(p::Preconditioner, b) = copyto!(b, p \ b)
     ldiv!(x, p::Preconditioner, b) = copyto!(x, p \ b)
     mul!(b, p::Preconditioner, x) = mul!(b, p.ml.levels[1].A, x)
 end
diff --git a/src/strength.jl b/src/strength.jl
@@ -5,12 +5,12 @@ end
 Classical(;θ = 0.25) = Classical(θ)
 
 function strength_of_connection(c::Classical, 
-                A::SparseMatrixCSC{Tv,Ti}) where {Ti,Tv}
+                At::SparseMatrixCSC{Tv,Ti}) where {Ti,Tv}
 
     θ = c.θ
 
-    m, n = size(A)
-    T = copy(A')
+    m, n = size(At)
+    T = deepcopy(At)
 
     for i = 1:n
         _m = find_max_off_diag(T, i)
@@ -29,7 +29,7 @@ function strength_of_connection(c::Classical,
 
         end
     end
-
+    
     dropzeros!(T)
 
     scale_cols_by_largest_entry!(T)
diff --git a/test/sa_tests.jl b/test/sa_tests.jl
@@ -13,12 +13,12 @@ function symmetric_soc(A::SparseMatrixCSC{T,V}, θ) where {T,V}
 
     S = sparse(i,j,v, size(A)...) + spdiagm(0=>D)
 
-    scale_cols_by_largest_entry!(S)
-
     for i = 1:size(S.nzval,1)
         S.nzval[i] = abs(S.nzval[i])
     end
 
+    scale_cols_by_largest_entry!(S)
+
     S
 end
 
@@ -222,7 +222,7 @@ function test_approximate_spectral_radius()
     for A in cases
         A = A + A'
         @static if VERSION < v"0.7-"
-            E,V = eig(A)            
+            E,V = eig(A)
         else
             E,V = (eigen(A)...,)
         end
@@ -231,15 +231,12 @@ function test_approximate_spectral_radius()
         expected_eig = E[largest_eig]
 
         @test isapprox(approximate_spectral_radius(A), expected_eig)
-
     end
-
 end
 
 # Test Gauss Seidel 
 import AlgebraicMultigrid: gs!, relax!
-function test_gauss_seidel()
-    
+function test_gauss_seidel()    
     N = 1
     A = spdiagm(0 => 2 * ones(N), -1 => -ones(N-1), 1 => -ones(N-1))
     x = eltype(A).(collect(0:N-1))
@@ -331,7 +328,5 @@ function test_symmetric_sweep()
     relax!(s, A, x, b)
     @test sum(abs2, x - [0.176765; 0.353529; 0.497517; 0.598914; 
                             0.653311; 0.659104; 0.615597; 0.52275; 
-                            0.382787; 0.203251]) < 1e-6
-        
+                            0.382787; 0.203251]) < 1e-6        
 end
-
