Add a few fixes to fit candidates routine and tests

Author: ranjanan

src/aggregation.jl

@@ -57,15 +57,14 @@ 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
     # relax!(improve_candidates, A, B, b)
     T, B = fit_candidates(AggOp, B)
 
-    P = smooth_prolongator(smooth, A, T, S, B)
+    P = smooth_prolongator(smooth, A, T, S', B)
     R = construct_R(symmetry, P)
     push!(levels, Level(A, P, R))
 
@@ -82,16 +81,20 @@ construct_R(::HermitianSymmetry, P) = P'
 function fit_candidates(AggOp, B, tol = 1e-10)
 
     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)
     Qx = zeros(eltype(B), nnz(A))
 
     copy!(Qx, B)
-    copy!(A.nzval, 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)
@@ -107,19 +110,20 @@ function fit_candidates(AggOp, B, tol = 1e-10)
             row = A.rowval[j]
             # Qx[row] *= scale
             #@show k
-            Qx[k] *= scale
-            k += 1
+            # Qx[k] *= scale
+            # k += 1
+            A.nzval[j] *= scale
         end
     end
 
-    SparseMatrixCSC(size(A)..., A.colptr, A.rowval, Qx), R
-    # A, 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]]
+        val = Qx[A.rowval[j]]
+        # val = A.nzval[A.rowval[j]]
         s += val*val
     end
     sqrt(s)

src/smoother.jl

@@ -111,7 +111,7 @@ 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
 
 #approximate_spectral_radius(A) =
@@ -127,4 +127,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/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/sa_tests.jl

@@ -176,10 +176,10 @@ function generate_fit_candidates_cases()
         B = ones(T, 9)
         push!(cases, (AggOp, B))
 
-        #AggOp = SparseMatrixCSC(3, 9, collect(1:10), 
-                        #[3,2,1,1,2,3,2,1,3], ones(T,9))
-        #B = T.(collect(1:9))
-        #push!(cases, (AggOp, B))
+        AggOp = SparseMatrixCSC(3, 9, collect(1:10), 
+                        [3,2,1,1,2,3,2,1,3], ones(T,9))
+        B = T.(collect(1:9))
+        push!(cases, (AggOp, B))
     end
 
     cases