Merge pull request #77 from Alexander-Barth/master

implement W and F cycles and n-dimensional Poisson matrix

README.md

@@ -5,7 +5,7 @@
 [![codecov.io](http://codecov.io/github/JuliaLinearAlgebra/AlgebraicMultigrid.jl/coverage.svg?branch=master)](http://codecov.io/github/JuliaLinearAlgebra/AlgebraicMultigrid.jl?branch=master)
 [![Build status](https://ci.appveyor.com/api/projects/status/0wnj4lhk1rvl5pjp?svg=true)](https://ci.appveyor.com/project/ranjanan/algebraicmultigrid-jl)
 
-This package lets you solve sparse linear systems using Algebraic Multigrid (AMG). This works especially well for symmetric positive definite matrices. 
+This package lets you solve sparse linear systems using Algebraic Multigrid (AMG). This works especially well for symmetric positive definite matrices.
 
 ## Usage
 
@@ -40,18 +40,18 @@ You can use AMG as a preconditioner for Krylov methods such as Conjugate Gradien
 ```julia
 import IterativeSolvers: cg
 p = aspreconditioner(ml)
-c = cg(A, A*ones(1000), Pl = p) 
+c = cg(A, A*ones(1000), Pl = p)
 ```
 
 ## Features and Roadmap
 
-This package currently supports: 
+This package currently supports:
 
 AMG Styles:
 * Ruge-Stuben Solver
 * Smoothed Aggregation (SA)
 
-Strength of Connection: 
+Strength of Connection:
 * Classical Strength of Connection
 * Symmetric Strength of Connection
 
@@ -60,12 +60,12 @@ Smoothers:
 * Damped Jacobi
 
 Cycling:
-* V cycle 
+* V, W and F cycles
 
 In the future, this package will support:
 1. Other splitting methods (like CLJP)
 2. SOR smoother
-3. W, F, AMLI cycles
+3. AMLI cycles
 
 #### Acknowledgements
-This package has been heavily inspired by the [`PyAMG`](http://github.com/pyamg/pyamg) project. 
+This package has been heavily inspired by the [`PyAMG`](http://github.com/pyamg/pyamg) project.

src/gallery.jl

@@ -1,3 +1,63 @@
-poisson(T, n) = sparse(Tridiagonal(fill(T(-1), n-1), 
+poisson(T, n) = sparse(Tridiagonal(fill(T(-1), n-1),
                         fill(T(2), n), fill(T(-1), n-1)))
 poisson(n) = poisson(Float64, n)
+
+function stencil_grid(T,stencil,sz)
+    # upper-bound for storage
+    n = prod(sz) * sum(.!iszero,stencil)
+
+    # indices and value of nonzero elements
+    Si = zeros(Int,n)
+    Sj = zeros(Int,n)
+    Ss = zeros(T,n)
+
+    linindices = LinearIndices(sz)
+    nnz = 0
+
+    stencil_sz = size(stencil)
+    offset = CartesianIndex((stencil_sz .+ 1) .÷ 2)
+
+    for i in CartesianIndices(sz)
+        for k in CartesianIndices(stencil_sz)
+            if stencil[k] != 0
+                j = i + k - offset
+                if checkbounds(Bool,linindices,j)
+                    nnz = nnz + 1
+                    Si[nnz] = linindices[i]
+                    Sj[nnz] = linindices[j]
+                    Ss[nnz] = stencil[k]
+                end
+            end
+        end
+    end
+
+    sparse((@view Si[1:nnz]),
+           (@view Sj[1:nnz]),
+           (@view Ss[1:nnz]),
+           prod(sz),prod(sz))
+end
+
+
+
+function poisson(T,sz::NTuple{N,Int}) where N
+    #=
+    In 2D stencil has the value:
+    stencil = [0  -1   0;
+               -1  4  -1;
+               0  -1   0]
+
+    =#
+    stencil = zeros(T,ntuple(i -> 3,Val(N)))
+    for i = 0:N-1
+        Ipre = CartesianIndex(ntuple(l -> 2,i))
+        Ipost = CartesianIndex(ntuple(l -> 2,N-i-1))
+        stencil[Ipre, 1, Ipost] = -1
+        stencil[Ipre, 3, Ipost] = -1
+    end
+    Icenter = CartesianIndex(ntuple(l -> 2,Val(N)))
+    stencil[Icenter] = 2*N
+
+    stencil_grid(T,stencil,sz)
+end
+
+poisson(sz::NTuple{N,Int}) where N = poisson(Float64,sz)

src/multilevel.jl

@@ -66,8 +66,8 @@ function Base.show(io::IO, ml::MultiLevel)
     op = operator_complexity(ml)
     g = grid_complexity(ml)
     c = ml.coarse_solver
-    total_nnz = nnz(ml.final_A) 
-    if !isempty(ml.levels) 
+    total_nnz = nnz(ml.final_A)
+    if !isempty(ml.levels)
         total_nnz += sum(nnz(level.A) for level in ml.levels)
     end
     lstr = ""
@@ -119,6 +119,12 @@ abstract type Cycle end
 struct V <: Cycle
 end
 
+struct W <: Cycle
+end
+
+struct F <: Cycle
+end
+
 """
     solve(ml::MultiLevel, b::AbstractArray, cycle, kwargs...)
 
@@ -140,7 +146,7 @@ Keyword Arguments
 
 """
 function solve(ml::MultiLevel, b::AbstractArray, args...; kwargs...)
-    n = length(ml) == 1 ? size(ml.final_A, 1) : size(ml.levels[1].A, 1) 
+    n = length(ml) == 1 ? size(ml.final_A, 1) : size(ml.levels[1].A, 1)
     V = promote_type(eltype(ml.workspace), eltype(b))
     x = zeros(V, size(b))
     return solve!(x, ml, b, args...; kwargs...)
@@ -183,8 +189,22 @@ function solve!(x, ml::MultiLevel, b::AbstractArray{T},
     # @show residuals
     log ? (x, residuals) : x
 end
-function __solve!(x, ml, v::V, b, lvl)
 
+function __solve_next!(x, ml, cycle::V, b, lvl)
+    __solve!(x, ml, cycle, b, lvl)
+end
+
+function __solve_next!(x, ml, cycle::W, b, lvl)
+    __solve!(x, ml, cycle, b, lvl)
+    __solve!(x, ml, cycle, b, lvl)
+end
+
+function __solve_next!(x, ml, cycle::F, b, lvl)
+    __solve!(x, ml, cycle, b, lvl)
+    __solve!(x, ml, V(), b, lvl)
+end
+
+function __solve!(x, ml, cycle::Cycle, b, lvl)
     A = ml.levels[lvl].A
     ml.presmoother(A, x, b)
 
@@ -200,7 +220,7 @@ function __solve!(x, ml, v::V, b, lvl)
     if lvl == length(ml.levels)
         ml.coarse_solver(coarse_x, coarse_b)
     else
-        coarse_x = __solve!(coarse_x, ml, v, coarse_b, lvl + 1)
+        coarse_x = __solve_next!(coarse_x, ml, cycle, coarse_b, lvl + 1)
     end
 
     mul!(res, ml.levels[lvl].P, coarse_x)

src/preconditioner.jl

@@ -1,10 +1,11 @@
-struct Preconditioner{ML<:MultiLevel}
+struct Preconditioner{ML<:MultiLevel,C<:Cycle}
     ml::ML
     init::Symbol
+    cycle::C
 end
-Preconditioner(ml) = Preconditioner(ml, :zero)
+Preconditioner(ml, cycle::Cycle) = Preconditioner(ml, :zero, cycle)
 
-aspreconditioner(ml::MultiLevel) = Preconditioner(ml)
+aspreconditioner(ml::MultiLevel, cycle::Cycle = V()) = Preconditioner(ml,cycle)
 
 import LinearAlgebra: \, *, ldiv!, mul!
 ldiv!(p::Preconditioner, b) = copyto!(b, p \ b)
@@ -14,7 +15,7 @@ function ldiv!(x, p::Preconditioner, b)
     else
         x .= b
     end
-    solve!(x, p.ml, b, maxiter = 1, calculate_residual = false)
+    solve!(x, p.ml, b, p.cycle, maxiter = 1, calculate_residual = false)
 end
 mul!(b, p::Preconditioner, x) = mul!(b, p.ml.levels[1].A, x)
 

test/cycle_tests.jl

@@ -0,0 +1,31 @@
+    import AlgebraicMultigrid: ruge_stuben, smoothed_aggregation,
+    poisson, aspreconditioner
+
+import IterativeSolvers: cg
+
+function test_cycles()
+
+    A = poisson((50,50))
+    b = A * ones(size(A,1))
+
+    reltol = 1e-8
+
+    for method in [ruge_stuben, smoothed_aggregation]
+        ml = method(A)
+
+        for cycle in [AlgebraicMultigrid.V(),AlgebraicMultigrid.W(),AlgebraicMultigrid.F()]
+            x,convhist = solve(ml, b, cycle; reltol = reltol, log = true)
+
+            @debug "number of iterations for $cycle using $method: $(length(convhist))"
+            @test norm(b - A*x) < reltol * norm(b)
+        end
+
+        for cycle in [AlgebraicMultigrid.V(),AlgebraicMultigrid.W(),AlgebraicMultigrid.F()]
+            p = aspreconditioner(ml,cycle)
+
+            x,convhist = cg(A, b, Pl = p; reltol = reltol, log = true)
+            @debug "CG: number of iterations for $cycle using $method: $(convhist.iters)"
+            @test norm(b - A*x) <= reltol * norm(b)
+        end
+    end
+end

test/runtests.jl

@@ -8,6 +8,7 @@ using FileIO
 using Random: seed!
 
 include("sa_tests.jl")
+include("cycle_tests.jl")
 
 @testset "AlgebraicMultigrid Tests" begin
 
@@ -139,6 +140,9 @@ ml = ruge_stuben(A)
 x = solve(ml, A * ones(100))
 @test sum(abs2, x - zeros(100)) < 1e-6
 
+
+
+
 end
 
 @testset "Preconditioning" begin
@@ -237,7 +241,7 @@ for (T,V) in ((Float64, Float64), (Float32,Float32),
         b = V.(b)
         c = cg(a, b, maxiter = 10)
         @test eltype(solve(ml, b)) == eltype(c)
-end    
+end
 
 end
 
@@ -271,6 +275,10 @@ end
     test_jacobi_prolongator()
 end
 
+@testset "Cycles" begin
+    test_cycles()
+end
+
 @testset "Int32 support" begin
     a = sparse(Int32.(1:10), Int32.(1:10), rand(10))
     @inferred smoothed_aggregation(a)