Merge pull request #29 from ranjanan/mem

Efficiency improvements, add multi-threaded matvec

Author: ranjanan

src/AMG.jl

@@ -3,6 +3,8 @@ module AMG
 import IterativeSolvers: gauss_seidel!
 using Base.Threads
 
+const MT = false
+
 include("utils.jl")
 export approximate_spectral_radius
 

src/smoother.jl

@@ -124,10 +124,26 @@ function weight(::DiagonalWeighting, S, ω)
 end
 
 function weight(::LocalWeighting, S, ω)
-    D = abs.(S) * ones(size(S, 1))
+    #=D = abs.(S) * ones(eltype(S), size(S, 1))
     D_inv = 1 ./ D[find(D)]
     D_inv_S = scale_rows(S, D_inv)
-    eltype(S)(ω) * D_inv_S
+    eltype(S)(ω) * D_inv_S=#
+    D = zeros(eltype(S), size(S,1))
+    for i = 1:size(S, 1)
+        for j in nzrange(S, i)
+            row = S.rowval[j]
+            val = S.nzval[j]
+            D[row] += abs(val)
+        end
+    end
+    for i = 1:size(D, 1)
+        if D[i] != 0
+            D[i] = 1/D[i]
+        end
+    end
+    D_inv_S = scale_rows(S, D)
+    # eltype(S)(ω) * D_inv_S
+    scale!(D_inv_S, eltype(S)(ω))
 end
 
 #approximate_spectral_radius(A) =

src/utils.jl

@@ -97,4 +97,45 @@ function approximate_eigenvalues(A, tol, maxiter, symmetric, v0)
 end
 
 find_breakdown(::Type{Float64}) = eps(Float64) * 10^6
-find_breakdown(::Type{Float32}) = eps(Float64) * 10^3
+find_breakdown(::Type{Float32}) = eps(Float64) * 10^3
+
+using Base.Threads
+#=function mul!(α::Number, A::SparseMatrixCSC, B::StridedVecOrMat, β::Number, C::StridedVecOrMat)
+    A.n == size(B, 1) || throw(DimensionMismatch())
+    A.m == size(C, 1) || throw(DimensionMismatch())
+    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
+    nzv = A.nzval
+    rv = A.rowval
+    if β != 1
+        β != 0 ? scale!(C, β) : fill!(C, zero(eltype(C)))
+    end
+    for k = 1:size(C, 2)
+        for col = 1:A.n
+            αxj = α*B[col,k]
+            for j = A.colptr[col]:(A.colptr[col + 1] - 1)
+                C[rv[j], k] += nzv[j]*αxj
+            end
+        end
+    end
+    C
+end
+spmatvec(A::SparseMatrixCSC{TA,S}, x::StridedVector{Tx}) where {TA,S,Tx} =
+    (T = promote_type(TA, Tx); mul!(one(T), A, x, zero(T), similar(x, T, A.m)))=#
+
+# export spmatvec
+
+@static if haskey(ENV, "JULIA_NUM_THREADS")
+    import Base: *
+    function *(A::SparseMatrixCSC{T,V}, b::Vector{T}) where {T,V}
+        m, n = size(A)
+        ret = zeros(T, m)
+        @threads for i = 1:n
+            for j in nzrange(A, i)
+                row = A.rowval[j]
+                val = A.nzval[j]
+                ret[row] += val * b[i]
+            end
+        end
+        ret
+    end
+end