Adding preconditioner to idrs (#297)

* Adding preconditioner to idrs

* Update test/idrs.jl

Author: deltaeecs

src/idrs.jl

@@ -24,6 +24,7 @@ shadow space.
 ## Keywords
 
 - `s::Integer = 8`: dimension of the shadow space;
+- `Pl::precT`: left preconditioner,
 - `abstol::Real = zero(real(eltype(b)))`,
   `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
   tolerance for the stopping condition
@@ -46,6 +47,7 @@ shadow space.
 """
 function idrs!(x, A, b;
                s = 8,
+               Pl = Identity(),
                abstol::Real = zero(real(eltype(b))),
                reltol::Real = sqrt(eps(real(eltype(b)))),
                maxiter=size(A, 2),
@@ -55,7 +57,7 @@ function idrs!(x, A, b;
     history[:abstol] = abstol
     history[:reltol] = reltol
     log && reserve!(history, :resnorm, maxiter)
-    idrs_method!(history, x, A, b, s, abstol, reltol, maxiter; kwargs...)
+    idrs_method!(history, x, A, b, s, Pl, abstol, reltol, maxiter; kwargs...)
     log && shrink!(history)
     log ? (x, history) : x
 end
@@ -78,8 +80,8 @@ end
 end
 
 function idrs_method!(log::ConvergenceHistory, X, A, C::T,
-    s::Number, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false
-    ) where {T}
+    s::Number, Pl::precT, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false
+    ) where {T, precT}
 
     verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm")
     R = C - A*X
@@ -129,6 +131,9 @@ function idrs_method!(log::ConvergenceHistory, X, A, C::T,
                 Q .+= c[i-k+1] .* U[i]
             end
 
+            # Preconditioning
+            ldiv!(Pl, V)
+
             # Compute new U[:,k] and G[:,k], G[:,k] is in space G_j
             V .= R .- V
 
@@ -181,6 +186,10 @@ function idrs_method!(log::ConvergenceHistory, X, A, C::T,
         # Now we have sufficient vectors in G_j to compute residual in G_j+1
         # Note: r is already perpendicular to P so v = r
         copyto!(V, R)
+
+        # Preconditioning
+        ldiv!(Pl, V)
+
         mul!(Q, A, V)
         om = omega(Q, R)
         R .-= om .* Q

test/idrs.jl

@@ -42,6 +42,25 @@ end
     @test norm(A * x - b) / norm(b) ≤ reltol
 end
 
+@testset "SparseMatrixCSC{$T, $Ti} with preconditioner" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
+    A = sprand(T, 1000, 1000, 0.1) + 30 * I
+    b = rand(T, 1000)
+    reltol = √eps(real(T))
+
+    x, history = idrs(A, b, log=true)
+    @test history.isconverged
+    @test norm(A * x - b) / norm(b) ≤ reltol
+
+    Apre = lu(A)
+    xpre, historypre = idrs(A, b, Pl = Apre, log=true)
+    @test historypre.isconverged
+    @test norm(A * xpre - b) / norm(b) ≤ reltol
+
+    @test isapprox(x, xpre, rtol = 1e-3)
+    @test historypre.iters < history.iters
+
+end
+
 @testset "Maximum number of iterations" begin
     x, history = idrs(rand(5, 5), rand(5), log=true, maxiter=2)
     @test history.iters == 2