callable Pinv and some inplace

Author: mohamed82008

src/aggregation.jl

@@ -10,8 +10,7 @@ function smoothed_aggregation(A::SparseMatrixCSC{T,V},
                         max_coarse = 10,
                         diagonal_dominance = false,
                         keep = false,
-                        coarse_solver = Pinv()) where {T,V}
-
+                        coarse_solver = Pinv) where {T,V}
 
     n = size(A, 1)
     # B = kron(ones(n, 1), eye(1))
@@ -42,7 +41,7 @@ function smoothed_aggregation(A::SparseMatrixCSC{T,V},
     #=A, B = extend_hierarchy!(levels, strength, aggregate, smooth,
                             improve_candidates, diagonal_dominance,
                             keep, A, B, symmetry)=#
-    MultiLevel(levels, A, presmoother, postsmoother)
+    MultiLevel(levels, A, coarse_solver(A), presmoother, postsmoother)
 end
 
 struct HermitianSymmetry

src/multilevel.jl

@@ -13,11 +13,14 @@ struct MultiLevel{S, Pre, Post, Ti, Tv}
 end
 
 abstract type CoarseSolver end
-struct Pinv <: CoarseSolver
+struct Pinv{T} <: CoarseSolver
+    pinvA::Matrix{T}
+    Pinv(A) = new{eltype(A)}(pinv(Matrix(A)))
 end
+(p::Pinv)(x, b) = mul!(x, p.pinvA, b)
 
 MultiLevel(l::Vector{Level{Ti,Tv}}, A::SparseMatrixCSC{Ti,Tv}, presmoother, postsmoother) where {Ti,Tv} =
-    MultiLevel(l, A, Pinv(), presmoother, postsmoother)
+    MultiLevel(l, A, Pinv(A), presmoother, postsmoother)
 Base.length(ml) = length(ml.levels) + 1
 
 function Base.show(io::IO, ml::MultiLevel)
@@ -113,9 +116,9 @@ function solve(ml::MultiLevel, b::AbstractVector{T},
     lvl = 1
     while length(residuals) <= maxiter && residuals[end] > tol
         if length(ml) == 1
-            x = coarse_solver(ml.coarse_solver, A, b)
+            ml.coarse_solver(x, b)
         else
-            x = __solve(cycle, ml, x, b, lvl)
+            __solve!(x, ml, cycle, b, lvl)
         end
         push!(residuals, T(norm(b - A * x)))
     end
@@ -127,7 +130,7 @@ function solve(ml::MultiLevel, b::AbstractVector{T},
         return x
     end
 end
-function __solve(v::V, ml, x, b, lvl)
+function __solve!(x, ml, v::V, b, lvl)
 
     A = ml.levels[lvl].A
     ml.presmoother(A, x, b)
@@ -137,9 +140,9 @@ function __solve(v::V, ml, x, b, lvl)
     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)
+        ml.coarse_solver(coarse_x, coarse_b)
     else
-        coarse_x = __solve(v, ml, coarse_x, coarse_b, lvl + 1)
+        coarse_x = __solve!(coarse_x, ml, v, coarse_b, lvl + 1)
     end
 
     x .+= ml.levels[lvl].P * coarse_x
@@ -148,5 +151,3 @@ function __solve(v::V, ml, x, b, lvl)
 
     x
 end
-
-coarse_solver(::Pinv, A, b) = pinv(Matrix(A)) * b

test/runtests.jl

@@ -1,7 +1,7 @@
 using Compat, Compat.Test, Compat.LinearAlgebra
 using Compat.SparseArrays, Compat.DelimitedFiles, Compat.Random
 using IterativeSolvers, FileIO, AlgebraicMultigrid
-import AlgebraicMultigrid: V, coarse_solver, Pinv, Classical
+import AlgebraicMultigrid: Pinv, Classical
 
 include("sa_tests.jl")
 
@@ -68,7 +68,7 @@ end
 @testset "Coarse Solver" begin
 A = float.(poisson(10))
 b = A * ones(10)
-@test sum(abs2, coarse_solver(Pinv(), A, b) - ones(10)) < 1e-6
+@test sum(abs2, Pinv(A)(similar(b), b) - ones(10)) < 1e-6
 end
 
 @testset "Multilevel" begin