Add user-defined near null space to aggregation-based AMG

.gitignore

@@ -2,3 +2,7 @@
 *.jl.*.cov
 *.jl.mem
 /Manifest.toml
+.vscode
+*.mtx
+*.csv
+*.geo

Project.toml

@@ -24,6 +24,9 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
+FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b"
 
 [targets]
-test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test"]
+test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test","Ferrite","FerriteGmsh","Downloads"]

src/AlgebraicMultigrid.jl

@@ -9,7 +9,7 @@ using Printf
 @reexport import CommonSolve: solve, solve!, init
 using Reexport
 
-using LinearAlgebra: rmul!
+using LinearAlgebra: rmul!,qr
 
 include("utils.jl")
 export approximate_spectral_radius

src/aggregation.jl

@@ -1,21 +1,21 @@
-function smoothed_aggregation(A::TA, 
-                        ::Type{Val{bs}}=Val{1};
-                        symmetry = HermitianSymmetry(),
-                        strength = SymmetricStrength(),
-                        aggregate = StandardAggregation(),
-                        smooth = JacobiProlongation(4.0/3.0),
-                        presmoother = GaussSeidel(),
-                        postsmoother = GaussSeidel(),
-                        improve_candidates = GaussSeidel(iter=4),
-                        max_levels = 10,
-                        max_coarse = 10,
-                        diagonal_dominance = false,
-                        keep = false,
-                        coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}}
+function smoothed_aggregation(A::TA, _B=nothing,
+    ::Type{Val{bs}}=Val{1};
+    symmetry=HermitianSymmetry(),
+    strength=SymmetricStrength(),
+    aggregate=StandardAggregation(),
+    smooth=JacobiProlongation(4.0 / 3.0),
+    presmoother=GaussSeidel(),
+    postsmoother=GaussSeidel(),
+    improve_candidates=GaussSeidel(iter=4),
+    max_levels=10,
+    max_coarse=10,
+    diagonal_dominance=false,
+    keep=false,
+    coarse_solver=Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}}
 
     n = size(A, 1)
-    # B = kron(ones(n, 1), eye(1))
-    B = ones(T,n)
+    B = isnothing(_B) ? ones(T, n, 1) : copy(_B)
+    @assert size(A, 1) == size(B, 1)
 
     #=max_levels, max_coarse, strength =
         levelize_strength_or_aggregation(max_levels, max_coarse, strength)
@@ -28,15 +28,15 @@ function smoothed_aggregation(A::TA,
     # agg = [aggregate for _ in 1:max_levels - 1]
     # sm = [smooth for _ in 1:max_levels]
 
-    levels = Vector{Level{TA, TA, Adjoint{T, TA}}}()
+    levels = Vector{Level{TA,TA,Adjoint{T,TA}}}()
     bsr_flag = false
     w = MultiLevelWorkspace(Val{bs}, eltype(A))
     residual!(w, size(A, 1))
 
     while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
         A, B, bsr_flag = extend_hierarchy!(levels, strength, aggregate, smooth,
-                                improve_candidates, diagonal_dominance,
-                                keep, A, B, symmetry, bsr_flag)
+            improve_candidates, diagonal_dominance,
+            keep, A, B, symmetry, bsr_flag)
         coarse_x!(w, size(A, 1))
         coarse_b!(w, size(A, 1))
         #=if size(A, 1) <= max_coarse
@@ -54,9 +54,9 @@ struct HermitianSymmetry
 end
 
 function extend_hierarchy!(levels, strength, aggregate, smooth,
-                            improve_candidates, diagonal_dominance, keep,
-                            A, B,
-                            symmetry, bsr_flag)
+    improve_candidates, diagonal_dominance, keep,
+    A, B,
+    symmetry, bsr_flag)
 
     # Calculate strength of connection matrix
     if symmetry isa HermitianSymmetry
@@ -70,7 +70,7 @@ function extend_hierarchy!(levels, strength, aggregate, smooth,
     # b = zeros(eltype(A), size(A, 1))
 
     # Improve candidates
-    b = zeros(size(A,1))
+    b = zeros(size(A, 1))
     improve_candidates(A, B, b)
     T, B = fit_candidates(AggOp, B)
 
@@ -88,59 +88,39 @@ function extend_hierarchy!(levels, strength, aggregate, smooth,
 end
 construct_R(::HermitianSymmetry, P) = P'
 
-function fit_candidates(AggOp, B, tol = 1e-10)
-
+function fit_candidates(AggOp, B; tol=1e-10)
     A = adjoint(AggOp)
-    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)
-    for i = 1:size(Qx, 1)
-        Qx[i] = B[i]
-    end
-    # 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)
-        threshold_i = tol * norm_i
-        if norm_i > threshold_i
-            scale = 1 / norm_i
-            R[i] = norm_i
-        else
-            scale = 0
-            R[i] = 0
-        end
-        for j in nzrange(A, i)
-            row = A.rowval[j]
-            # Qx[row] *= scale
-            #@show k
-            # Qx[k] *= scale
-            # k += 1
-            A.nzval[j] *= scale
+    n_fine, m = ndims(B) == 1 ? (length(B), 1) : size(B)
+    n_fine2, n_agg = size(A)
+    @assert n_fine2 == n_fine
+    n_coarse = m * n_agg
+    T = eltype(B)
+    Qs = spzeros(T, n_fine, n_coarse)
+    R = zeros(T, n_coarse, m)
+
+    for agg in 1:n_agg
+        rows = A.rowval[A.colptr[agg]:A.colptr[agg+1]-1]
+        M = @view B[rows, :]     # size(rows) × m
+
+
+        F = qr(M)
+        r = min(length(rows), m)
+        Qfull = Matrix(F.Q)
+        Qj = Qfull[:, 1:r]
+        Rj = F.R
+
+        offset = (agg - 1) * m
+
+        for local_i in 1:length(rows), local_j in 1:r
+            val = Qj[local_i, local_j]
+            if abs(val) >= tol
+                Qs[rows[local_i], offset+local_j] = val
+            end
         end
-    end
+        dropzeros!(Qs)
 
-    # 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)
-        if A.rowval[j] > length(Qx)
-            val = 1
-        else
-            val = Qx[A.rowval[j]]
-        end
-        # val = A.nzval[A.rowval[j]]
-        s += val*val
+        R[offset+1:offset+r, :] .= Rj[1:r, :]
     end
-    sqrt(s)
+
+    return Qs, R
 end

src/multilevel.jl

@@ -287,7 +287,12 @@ end
 abstract type AMGAlg end
 
 struct RugeStubenAMG  <: AMGAlg end
-struct SmoothedAggregationAMG  <: AMGAlg end
+struct SmoothedAggregationAMG  <: AMGAlg
+    B::Union{<:AbstractMatrix,Nothing}
+    function SmoothedAggregationAMG(B::Union{AbstractMatrix,Nothing}=nothing)
+        new(B)
+    end
+end
 
 function solve(A::AbstractMatrix, b::Vector, s::AMGAlg, args...; kwargs...)
     solt = init(s, A, b, args...; kwargs...)
@@ -296,9 +301,9 @@ end
 function init(::RugeStubenAMG, A, b, args...; kwargs...)
     AMGSolver(ruge_stuben(A; kwargs...), b)
 end
-function init(::SmoothedAggregationAMG, A, b; kwargs...)
-    AMGSolver(smoothed_aggregation(A; kwargs...), b)
+function init(sa::SmoothedAggregationAMG, A, b; kwargs...)
+    AMGSolver(smoothed_aggregation(A,sa.B; kwargs...), b)
 end
 function solve!(solt::AMGSolver, args...; kwargs...) 
     _solve(solt.ml, solt.b, args...; kwargs...)   
-end
+end