Merge pull request #177 from haampie/remove-solve-function

Remove solve function

Author: haampie

benchmark/benchmark-hessenberg.jl

@@ -43,7 +43,7 @@ function backslash_versus_givens(; n = 1_000, ms = 10 : 10 : 100)
         rhs = [i == 1 ? 1.0 : 0.0 for i = 1 : size(H, 1)]
 
         # Run the benchmark
-        results["givens_qr"][m] = @benchmark IterativeSolvers.solve!(IterativeSolvers.Hessenberg(myH), myRhs) setup = (myH = copy($H); myRhs = copy($rhs))
+        results["givens_qr"][m] = @benchmark A_ldiv_B!(IterativeSolvers.Hessenberg(myH), myRhs) setup = (myH = copy($H); myRhs = copy($rhs))
         results["backslash"][m] = @benchmark $H \ $rhs
     end
 

src/cg.jl

@@ -66,7 +66,7 @@ end
 #####################
 
 function next(it::PCGIterable, iteration::Int)
-    solve!(it.c, it.Pl, it.r)
+    A_ldiv_B!(it.c, it.Pl, it.r)
 
     ρ_prev = it.ρ
     it.ρ = dot(it.c, it.r)

src/chebyshev.jl

@@ -29,7 +29,7 @@ done(c::ChebyshevIterable, iteration::Int) = iteration ≥ c.maxiter || converge
 function next(cheb::ChebyshevIterable, iteration::Int)
     T = eltype(cheb.x)
 
-    solve!(cheb.c, cheb.Pl, cheb.r)
+    A_ldiv_B!(cheb.c, cheb.Pl, cheb.r)
 
     if iteration == 1
         cheb.α = T(2) / cheb.λ_avg

src/common.jl

@@ -16,18 +16,9 @@ Build a zeros vector `Vector{T}`, where `T` is `Adivtype(A,b)`.
 """
 zerox(A, b) = zeros(Adivtype(A, b), size(A, 2))
 
-#### Numerics
 """
-    solve(A,b)
-Solve `A\\b` with a direct solver. When `A` is a function `A(b)` is dispatched instead.
+No-op preconditioner
 """
-solve(A::Function,b) = A(b)
-solve(A,b) = A\b
-solve!(out::AbstractArray{T},A::Int,b::AbstractArray{T}) where {T} = scale!(out,b, 1/A)
-solve!(out::AbstractArray{T},A,b::AbstractArray{T}) where {T} = A_ldiv_B!(out,A,b)
-solve!(out::AbstractArray{T},A::Function,b::AbstractArray{T}) where {T} = copy!(out,A(b))
-
-# Identity preconditioner
 struct Identity end
 
 \(::Identity, x) = copy(x)

test/cg.jl

@@ -2,8 +2,17 @@ using IterativeSolvers
 using LinearMaps
 using Base.Test
 
+import Base.A_ldiv_B!
+
 include("laplace_matrix.jl")
 
+struct JacobiPrec
+    diagonal
+end
+
+
+A_ldiv_B!(y, P::JacobiPrec, x) = y .= x ./ P.diagonal
+
 @testset "Conjugate Gradients" begin
 
 srand(1234321)
@@ -40,48 +49,35 @@ srand(1234321)
 end
 
 @testset "Sparse Laplacian" begin
-    A = laplace_matrix(Float64, 10, 3)
-    L = tril(A)
-    D = diag(A)
-    U = triu(A)
-
-    JAC(x) = D .\ x
-    SGS(x) = L \ (D .* (U \ x))
+    A = laplace_matrix(Float64, 10, 2)
+    P = JacobiPrec(diag(A))
 
     rhs = randn(size(A, 2))
-    rhs /= norm(rhs)
+    scale!(rhs, inv(norm(rhs)))
     tol = 1e-5
 
     @testset "Matrix" begin
         xCG = cg(A, rhs; tol=tol, maxiter=100)
-        xJAC = cg(A, rhs; Pl=JAC, tol=tol, maxiter=100)
-        xSGS = cg(A, rhs; Pl=SGS, tol=tol, maxiter=100)
+        xJAC = cg(A, rhs; Pl=P, tol=tol, maxiter=100)
         @test norm(A * xCG - rhs) ≤ tol
-        @test norm(A * xSGS - rhs) ≤ tol
         @test norm(A * xJAC - rhs) ≤ tol
     end
 
     Af = LinearMap(A)
     @testset "Function" begin
         xCG = cg(Af, rhs; tol=tol, maxiter=100)
-        xJAC = cg(Af, rhs; Pl=JAC, tol=tol, maxiter=100)
-        xSGS = cg(Af, rhs; Pl=SGS, tol=tol, maxiter=100)
+        xJAC = cg(Af, rhs; Pl=P, tol=tol, maxiter=100)
         @test norm(A * xCG - rhs) ≤ tol
-        @test norm(A * xSGS - rhs) ≤ tol
         @test norm(A * xJAC - rhs) ≤ tol
     end
 
     @testset "Function with specified starting guess" begin
-        tol = 1e-4
         x0 = randn(size(rhs))
         xCG, hCG = cg!(copy(x0), Af, rhs; tol=tol, maxiter=100, log=true)
-        xJAC, hJAC = cg!(copy(x0), Af, rhs; Pl=JAC, tol=tol, maxiter=100, log=true)
-        xSGS, hSGS = cg!(copy(x0), Af, rhs; Pl=SGS, tol=tol, maxiter=100, log=true)
+        xJAC, hJAC = cg!(copy(x0), Af, rhs; Pl=P, tol=tol, maxiter=100, log=true)
         @test norm(A * xCG - rhs) ≤ tol
-        @test norm(A * xSGS - rhs) ≤ tol
         @test norm(A * xJAC - rhs) ≤ tol
         @test niters(hJAC) == niters(hCG)
-        @test niters(hSGS) ≤ niters(hJAC)
     end
 end
 

test/hessenberg.jl

@@ -1,7 +1,7 @@
 using IterativeSolvers
 using Base.Test
 
-import IterativeSolvers: Hessenberg, solve!
+import IterativeSolvers: Hessenberg
 
 @testset "Hessenberg" begin