Merge pull request #9 from ranjanan/precond

Add option to use as preconditioner

Author: ranjanan

src/AMG.jl

@@ -19,4 +19,7 @@ export solve
 include("classical.jl")
 export ruge_stuben
 
+include("preconditioner.jl")
+export aspreconditioner
+
 end # module

src/multilevel.jl

@@ -65,6 +65,9 @@ function solve{T}(ml::MultiLevel, b::Vector{T}; maxiter = 100,
     residuals = Vector{T}()
     A = ml.levels[1].A
     normb = norm(b)
+    if normb != 0
+        tol *= normb
+    end
     push!(residuals, norm(b - A*x))
 
     lvl = 1
@@ -74,7 +77,7 @@ function solve{T}(ml::MultiLevel, b::Vector{T}; maxiter = 100,
         else
             x = __solve(cycle, ml, x, b, lvl)
         end
-        push!(residuals, norm(A*x - b))
+        push!(residuals, norm(b - A * x))
     end
     x
 end

src/preconditioner.jl

@@ -0,0 +1,13 @@
+import Base: \, *, A_ldiv_B!, A_mul_B!
+
+struct Preconditioner
+    ml::MultiLevel
+end
+
+aspreconditioner(ml::MultiLevel) = Preconditioner(ml)
+
+\(p::Preconditioner, b) = p * b
+*(p::Preconditioner, b) = solve(p.ml, b; cycle = V(), maxiter = 1, tol = 1e-12)
+
+A_ldiv_B!(x, p::Preconditioner, b) = copy!(x, p \ b)
+A_mul_B!(b, p::Preconditioner, x) = A_mul_B!(b, p.ml.levels[1].A, x)

test/runtests.jl

@@ -1,6 +1,7 @@
 using AMG
 using Base.Test
 using JLD
+using IterativeSolvers
 
 @testset "AMG Tests" begin
 
@@ -111,4 +112,55 @@ x = solve(ml, A * ones(100))
 
 end
 
+@testset "Preconditioning" begin
+A = load("thing.jld")["G"]
+n = size(A, 1)
+ml = ruge_stuben(A)
+p = aspreconditioner(ml)
+b = zeros(n)
+b[1] = 1
+b[2] = -1
+x = solve(p.ml, A * ones(n); maxiter = 1, tol = 1e-12)
+diff = x - [  1.88664780e-16,   2.34982727e-16,   2.33917697e-16,
+         8.77869044e-17,   7.16783490e-17,   1.43415460e-16,
+         3.69199021e-17,   9.70950385e-17,   4.77034895e-17,
+         3.77491328e-17,   5.07592420e-18,   6.32131628e-19,
+        -1.60361276e-18,  -1.36749626e-16,  -4.16651794e-17,
+        -3.48207590e-17,  -4.19334783e-17,  -4.60500098e-17,
+        -6.91113945e-17,  -1.35997904e-16,  -9.85940056e-17,
+         8.99433377e-17,   1.63924842e-16,  -2.43048120e-16,
+        -3.21888624e-16,  -1.56389252e-16,  -3.66495834e-17,
+         1.70606350e-16,   1.66788345e-16,  -3.26736922e-16,
+        -3.39591125e-16,  -3.67075849e-16,  -3.89891523e-16,
+        -4.34537627e-16,  -4.42863579e-16,  -6.62433333e-16,
+        -5.56056397e-16,  -5.70242981e-16,  -6.48075679e-16,
+        -6.58200572e-16,  -7.24886474e-16,  -7.55973538e-16,
+        -6.76965535e-16,  -7.00643227e-16,  -6.23581397e-16,
+        -7.03016682e-16]
+@test sum(abs2, diff) < 1e-8
+x = solve(p.ml, b; maxiter = 1, tol = 1e-12)
+diff = x - [ 0.76347046, -0.5498286 , -0.2705487 , -0.15047352, -0.10248021,
+        0.60292674, -0.11497073, -0.08460548, -0.06931461,  0.38230708,
+       -0.055664  , -0.04854558, -0.04577031,  0.09964325,  0.01825624,
+       -0.01990265, -0.02866185, -0.03049521,  0.03310897, -0.01709034,
+       -0.02038031, -0.01325201, -0.01051535,  0.02992818,  0.01493605,
+       -0.00633922, -0.01285614, -0.01155069, -0.01095907,  0.04415807,
+        0.02213755,  0.018686  ,  0.02625713,  0.02007781,  0.01898018,
+        0.02107552,  0.01909623,  0.01874986,  0.01852736,  0.01844719,
+        0.01841821,  0.01841695,  0.01953195,  0.01885713,  0.01864432,
+        0.0185079 ]
+@test sum(abs2, diff) < 1e-8
+x = cg(A, b, Pl = p)
+diff = x - [ 0.82365077, -0.537589  , -0.30632349, -0.19370186, -0.14773294,
+        0.68489145, -0.15550115, -0.1278148 , -0.11197922,  0.45362483,
+       -0.08577219, -0.08598307, -0.08477946,  0.12985118,  0.02805496,
+       -0.03907565, -0.05950957, -0.06544269,  0.05446686, -0.047537  ,
+       -0.05203899, -0.04685981, -0.04491762,  0.05639249,  0.02792704,
+       -0.02282528, -0.04062864, -0.04321821, -0.0441893 ,  0.07593055,
+        0.05212038,  0.04464215,  0.05835841,  0.05079815,  0.04830733,
+        0.05272397,  0.05028666,  0.0494817 ,  0.04960952,  0.0496615 ,
+        0.04968258,  0.04968737,  0.05105749,  0.05009268,  0.04972329,
+        0.04970173]
+@test sum(abs2, diff) < 1e-8
+end
 end