v0.7 support, cleanup, comments and bug fixes

Author: mohamed82008

REQUIRE

@@ -1,2 +1,3 @@
-julia 0.6
-IterativeSolvers 0.4.1
+julia 0.7
+IterativeSolvers 0.7.0
+Compat 1.0.0

src/AMG.jl

@@ -1,6 +1,8 @@
 module AMG
 
 import IterativeSolvers: gauss_seidel!
+using Compat, Compat.LinearAlgebra
+using Compat.SparseArrays, Compat.Printf
 using Base.Threads
 
 const MT = false

src/aggregate.jl

@@ -97,7 +97,7 @@ function aggregation(::StandardAggregation, S::SparseMatrixCSC{T,R}) where {T,R}
     if minimum(x) == -1
         mask = x .!= -1
         I = collect(R, 1:n)[mask]
-        J = x[mask] + R(1)
+        J = x[mask] .+ R(1)
         #J = x[mask] + 1
         V = ones(eltype(S), length(J))
         AggOp = sparse(J,I,V,N,M)

src/aggregation.jl

@@ -77,11 +77,11 @@ function extend_hierarchy!(levels, strength, aggregate, smooth,
 
     A, B, bsr_flag
 end
-construct_R(::HermitianSymmetry, P) = P'
+construct_R(::HermitianSymmetry, P) = copy(P')
 
 function fit_candidates(AggOp, B, tol = 1e-10)
 
-    A = AggOp.'
+    A = copy(AggOp')
     n_fine, n_coarse = size(A)
     n_col = n_coarse
 

src/classical.jl

@@ -7,13 +7,13 @@ struct Solver{S,T,P,PS}
     max_coarse::Int64
 end
 
-function ruge_stuben{Ti,Tv}(A::SparseMatrixCSC{Ti,Tv};
+function ruge_stuben(A::SparseMatrixCSC{Ti,Tv};
                 strength = Classical(0.25),
                 CF = RS(),
                 presmoother = GaussSeidel(),
                 postsmoother = GaussSeidel(),
                 max_levels = 10,
-                max_coarse = 10)
+                max_coarse = 10) where {Ti,Tv}
 
     s = Solver(strength, CF, presmoother,
                 postsmoother, max_levels, max_levels)
@@ -29,7 +29,7 @@ function ruge_stuben{Ti,Tv}(A::SparseMatrixCSC{Ti,Tv};
     MultiLevel(levels, A, presmoother, postsmoother)
 end
 
-function extend_heirarchy!{Ti,Tv}(levels::Vector{Level{Ti,Tv}}, strength, CF, A::SparseMatrixCSC{Ti,Tv})
+function extend_heirarchy!(levels::Vector{Level{Ti,Tv}}, strength, CF, A::SparseMatrixCSC{Ti,Tv}) where {Ti,Tv}
     S, T = strength_of_connection(strength, A)
     splitting = split_nodes(CF, S)
     P, R = direct_interpolation(A, T, splitting)
@@ -38,50 +38,43 @@ function extend_heirarchy!{Ti,Tv}(levels::Vector{Level{Ti,Tv}}, strength, CF, A:
 end
 
 function direct_interpolation(A, T, splitting)
-
-    fill!(T.nzval, eltype(A)(1.))
+    fill!(T.nzval, eltype(A)(1))
     T .= A .* T
-    Pp = rs_direct_interpolation_pass1(T, A, splitting)
-    Pp .= Pp .+ 1
-
+    Pp = rs_direct_interpolation_pass1(T, splitting)
     Px, Pj, Pp = rs_direct_interpolation_pass2(A, T, splitting, Pp)
 
-    Pj .= Pj .+ 1
-
     R = SparseMatrixCSC(maximum(Pj), size(A, 1), Pp, Pj, Px)
-    P = R'
+    P = copy(R')
 
     P, R
 end
 
-
-function rs_direct_interpolation_pass1(T, A, splitting)
-
-     Bp = zeros(Int, size(A.colptr))
-     n = size(A, 1)
-     nnz = 0
-     for i = 1:n
-         if splitting[i] == C_NODE
-             nnz += 1
-         else
-             for j in nzrange(T, i)
-                 row = T.rowval[j]
-                 if splitting[row] == C_NODE && row != i
-                     nnz += 1
+# calculates the number of nonzeros in each column of the interpolation matrix
+function rs_direct_interpolation_pass1(T, splitting)
+    n = size(T, 2)
+    Bp = ones(Int, n+1)
+    nnzplus1 = 1
+    for i = 1:n
+        if splitting[i] == C_NODE
+            nnzplus1 += 1
+        else
+            for j in nzrange(T, i)
+                row = T.rowval[j]
+                if splitting[row] == C_NODE && row != i
+                    nnzplus1 += 1
                 end
-             end
+            end
         end
-        Bp[i+1] = nnz
-     end
-     Bp
+        Bp[i+1] = nnzplus1
+    end
+    Bp
  end
 
 
- function rs_direct_interpolation_pass2(A::SparseMatrixCSC{Tv,Ti},
+function rs_direct_interpolation_pass2(A::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)
@@ -157,7 +150,7 @@ function rs_direct_interpolation_pass1(T, A, splitting)
         m[i] = sum
         sum += splitting[i]
     end
-    Bj .= m[Bj]
+    Bj .= m[Bj] .+ 1
 
-   Bx, Bj, Bp
+    Bx, Bj, Bp
 end

src/multilevel.jl

@@ -16,7 +16,7 @@ abstract type CoarseSolver end
 struct Pinv <: CoarseSolver
 end
 
-MultiLevel{Ti,Tv}(l::Vector{Level{Ti,Tv}}, A::SparseMatrixCSC{Ti,Tv}, presmoother, postsmoother) =
+MultiLevel(l::Vector{Level{Ti,Tv}}, A::SparseMatrixCSC{Ti,Tv}, presmoother, postsmoother) where {Ti,Tv} =
     MultiLevel(l, A, Pinv(), presmoother, postsmoother)
 Base.length(ml) = length(ml.levels) + 1
 
@@ -92,12 +92,12 @@ Keyword Arguments
 * log::Bool - return vector of residuals along with solution
 
 """
-function solve{T}(ml::MultiLevel, b::Vector{T},
+function solve(ml::MultiLevel, b::Vector{T},
                                     cycle::Cycle = V();
                                     maxiter::Int = 100,
                                     tol::Float64 = 1e-5,
                                     verbose::Bool = false,
-                                    log::Bool = false)
+                                    log::Bool = false) where {T}
 
     A = length(ml) == 1 ? ml.final_A : ml.levels[1].A                                   
     V = promote_type(eltype(A), eltype(b))
@@ -149,4 +149,4 @@ function __solve(v::V, ml, x, b, lvl)
     x
 end
 
-coarse_solver(::Pinv, A, b) = pinv(full(A)) * b
+coarse_solver(::Pinv, A, b) = pinv(Matrix(A)) * b

src/preconditioner.jl

@@ -1,4 +1,4 @@
-import Base: \, *, A_ldiv_B!, A_mul_B!
+import LinearAlgebra: \, *, ldiv!, mul!
 
 struct Preconditioner
     ml::MultiLevel
@@ -9,5 +9,5 @@ aspreconditioner(ml::MultiLevel) = Preconditioner(ml)
 \(p::Preconditioner, b) = solve(p.ml, b, maxiter = 1, tol = 1e-12)
 *(p::Preconditioner, b) = p.ml.levels[1].A * x
 
-A_ldiv_B!(x, p::Preconditioner, b) = copy!(x, p \ b)
-A_mul_B!(b, p::Preconditioner, x) = A_mul_B!(b, p.ml.levels[1].A, x)
+ldiv!(x, p::Preconditioner, b) = copyto!(x, p \ b)
+mul!(b, p::Preconditioner, x) = mul!(b, p.ml.levels[1].A, x)

src/smoother.jl

@@ -143,7 +143,7 @@ function weight(::LocalWeighting, S, ω)
     end
     D_inv_S = scale_rows(S, D)
     # eltype(S)(ω) * D_inv_S
-    scale!(D_inv_S, eltype(S)(ω))
+    rmul!(D_inv_S, eltype(S)(ω))
 end
 
 #approximate_spectral_radius(A) =