Merge pull request #21 from ranjanan/single

Support single precision

Author: ranjanan

src/aggregate.jl

@@ -95,7 +95,7 @@ function aggregation(::StandardAggregation, S)
 
     Tp = collect(1:n+1)
     x .= x .+ 1
-    Tx = ones(length(x))
+    Tx = ones(eltype(S), length(x))
 
     SparseMatrixCSC(N, M, Tp, x, Tx)
 end

src/aggregation.jl

@@ -1,4 +1,4 @@
-function smoothed_aggregation{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti},
+function smoothed_aggregation(A::SparseMatrixCSC{T,V},
                         symmetry = HermitianSymmetry(),
                         strength = SymmetricStrength(),
                         aggregate = StandardAggregation(),
@@ -10,12 +10,12 @@ function smoothed_aggregation{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti},
                         max_coarse = 10,
                         diagonal_dominance = false,
                         keep = false,
-                        coarse_solver = Pinv())
+                        coarse_solver = Pinv()) where {T,V}
 
 
     n = size(A, 1)
     # B = kron(ones(n, 1), eye(1))
-    B = ones(n)
+    B = ones(T,n)
 
     #=max_levels, max_coarse, strength =
         levelize_strength_or_aggregation(max_levels, max_coarse, strength)
@@ -28,7 +28,7 @@ function smoothed_aggregation{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti},
     # agg = [aggregate for _ in 1:max_levels - 1]
     # sm = [smooth for _ in 1:max_levels]
 
-    levels = Vector{Level{Tv,Ti}}()
+    levels = Vector{Level{T,V}}()
     bsr_flag = false
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse

src/classical.jl

@@ -39,7 +39,7 @@ end
 
 function direct_interpolation(A, T, splitting)
 
-    fill!(T.nzval, 1.)
+    fill!(T.nzval, eltype(A)(1.))
     T .= A .* T
     Pp = rs_direct_interpolation_pass1(T, A, splitting)
     Pp .= Pp .+ 1
@@ -77,13 +77,13 @@ function rs_direct_interpolation_pass1(T, A, splitting)
  end
 
 
- function rs_direct_interpolation_pass2{Tv, Ti}(A::SparseMatrixCSC{Tv,Ti},
+ function rs_direct_interpolation_pass2(A::SparseMatrixCSC{Tv,Ti},
                                                 T::SparseMatrixCSC{Tv, Ti},
                                                 splitting::Vector{Ti},
-                                                Bp::Vector{Ti})
+                                                Bp::Vector{Ti}) where {Tv,Ti}
 
-
-    Bx = zeros(Float64, Bp[end] - 1)
+                                                
+    Bx = zeros(Tv, Bp[end] - 1)
     Bj = zeros(Ti, Bp[end] - 1)
 
     n = size(A, 1)
@@ -152,7 +152,7 @@ function rs_direct_interpolation_pass1(T, A, splitting)
     end
 
     m = zeros(Ti, n)
-    sum = zero(eltype(m))
+    sum = zero(Ti)
     for i = 1:n
         m[i] = sum
         sum += splitting[i]

src/multilevel.jl

@@ -1,7 +1,7 @@
-struct Level{Ti,Tv}
-    A::SparseMatrixCSC{Ti,Tv}
-    P::SparseMatrixCSC{Ti,Tv}
-    R::SparseMatrixCSC{Ti,Tv}
+struct Level{T,V}
+    A::SparseMatrixCSC{T,V}
+    P::SparseMatrixCSC{T,V}
+    R::SparseMatrixCSC{T,V}
 end
 
 struct MultiLevel{S, Pre, Post, Ti, Tv}
@@ -66,8 +66,10 @@ function solve{T}(ml::MultiLevel, b::Vector{T},
                                     tol = 1e-5;
                                     verbose = false,
                                     log = false)
-    x = zeros(T, size(b))
-    residuals = Vector{T}()
+    V = promote_type(eltype(ml.levels[1].A), eltype(b))
+    x = zeros(V, size(b))
+    tol = eltype(b)(tol)
+    residuals = Vector{V}()
     A = length(ml) == 1 ? ml.final_A : ml.levels[1].A
     normb = norm(b)
     if normb != 0
@@ -91,14 +93,14 @@ function solve{T}(ml::MultiLevel, b::Vector{T},
         return x
     end
 end
-function __solve{T}(v::V, ml, x::Vector{T}, b::Vector{T}, lvl)
+function __solve(v::V, ml, x, b, lvl)
 
     A = ml.levels[lvl].A
     presmoother!(ml.presmoother, A, x, b)
 
     res = b - A * x
     coarse_b = ml.levels[lvl].R * res
-    coarse_x = zeros(T, size(coarse_b))
+    coarse_x = zeros(eltype(coarse_b), size(coarse_b))
 
     if lvl == length(ml.levels)
         coarse_x = coarse_solver(ml.coarse_solver, ml.final_A, coarse_b)
@@ -113,5 +115,4 @@ function __solve{T}(v::V, ml, x::Vector{T}, b::Vector{T}, lvl)
     x
 end
 
-coarse_solver{Tv,Ti}(::Pinv, A::SparseMatrixCSC{Tv,Ti}, b::Vector{Tv}) =
-                                                        pinv(full(A)) * b
+coarse_solver(::Pinv, A, b) = pinv(full(A)) * b

src/smoother.jl

@@ -32,7 +32,7 @@ smoother!(s::GaussSeidel, ::BackwardSweep, A, x, b) =
     gs!(A, b, x, size(A,1), -1, 1)
 
 
-function gs!{T,Ti}(A::SparseMatrixCSC{T,Ti}, b::Vector{T}, x::Vector{T}, start, step, stop)
+function gs!(A, b, x, start, step, stop)
     n = size(A, 1)
     z = zero(eltype(A))
     for i = start:step:stop
@@ -110,7 +110,7 @@ end
 function weight(::DiagonalWeighting, S, ω)
     D_inv = 1 ./ diag(S)
     D_inv_S = scale_rows(S, D_inv)
-    (ω / approximate_spectral_radius(D_inv_S)) * D_inv_S
+    (eltype(S)(ω) / approximate_spectral_radius(D_inv_S)) * D_inv_S
     #(ω) * D_inv_S
 end
 

src/strength.jl

@@ -4,7 +4,8 @@ struct Classical{T} <: Strength
 end
 Classical(;θ = 0.25) = Classical(θ)
 
-function strength_of_connection{T, Ti, Tv}(c::Classical{T}, A::SparseMatrixCSC{Tv, Ti})
+function strength_of_connection(c::Classical{T}, 
+                A::SparseMatrixCSC{Tv,Ti}) where {T,Ti,Tv}
 
     θ = c.θ
 
@@ -77,13 +78,11 @@ SymmetricStrength() = SymmetricStrength(0.)
 
 function strength_of_connection{T}(s::SymmetricStrength{T}, A, bsr_flag = false)
 
-
-
     θ = s.θ
 
     if bsr_flag && θ == 0
         S = SparseMatrixCSC(size(A)...,
-                    A.colptr, A.rowval, ones(size(A.rowval)))
+                    A.colptr, A.rowval, ones(eltype(A), size(A.rowval)))
         return S
     else
         S = deepcopy(A)

src/utils.jl

@@ -5,11 +5,11 @@ function approximate_spectral_radius(A, tol = 0.01,
     symmetric = false
 
     # Initial guess
-    v0 = rand(size(A,1))
+    v0 = rand(eltype(A), size(A,1))
     maxiter = min(size(A, 1), maxiter)
     ev = zeros(eltype(A), maxiter)
     max_index = 0
-    X = zeros(size(A,1), maxiter)
+    X = zeros(eltype(A), size(A,1), maxiter)
 
     for i in 1:restart+1
         evect, ev, H, V, flag =
@@ -90,8 +90,7 @@ function approximate_eigenvalues(A, tol, maxiter, symmetric, v0)
         scale!(w, 1/H[j+1,j])
         push!(V, w)
     end
-
-    Eigs, Vects = eig(H[1:maxiter, 1:maxiter], eye(maxiter))
+    Eigs, Vects = eig(H[1:maxiter, 1:maxiter], eye(eltype(A), maxiter))
 
     Vects, Eigs, H, V, flag
 end

test/runtests.jl

@@ -51,6 +51,7 @@ end
 # Direct Interpolation
 using AMG
 A = poisson(5)
+A = Float64.(A)
 splitting = [1,0,1,0,1]
 P, R = AMG.direct_interpolation(A, copy(A), splitting)
 @test P ==  [ 1.0  0.0  0.0
@@ -216,4 +217,21 @@ diff = x - [0.775725, -0.571202, -0.290989, -0.157001, -0.106981, 0.622652,
 
 end
 
+@testset "Precision" begin
+
+a = poisson(100)
+b = rand(size(a,1))
+
+# Iterate through all types
+for (T,V) in ((Float64, Float64), (Float32,Float32),
+        (Float64,Float32), (Float32,Float64))
+        a = T.(a)
+        ml = smoothed_aggregation(a)
+        b = V.(b)
+        c = cg(a, b, maxiter = 10)
+        @test eltype(solve(ml, b)) == eltype(c)
+end    
+
 end
+
+end