Update, not reset the solution at restart (#119)

* Update, not reset the solution at restart

* Add test for optimality condition GMRES

Author: haampie

src/gmres.jl

@@ -114,7 +114,7 @@ function gmres_method!(log::ConvergenceHistory, x, A, b;
 
         @eval a = $(VERSION < v"0.4-" ? Triangular(H[1:N, 1:N], :U) \ s[1:N] : UpperTriangular(H[1:N, 1:N]) \ s[1:N])
         w = a[1:N] * K
-        @blas! x = solve(Pr,w) #Right preconditioner
+        @blas! x += solve(Pr,w) #Right preconditioner
 
         if (rho<tol) | ((macroiter-1)*restart+N >= maxiter)
             setconv(log, rho<tol)

test/runtests.jl

@@ -74,6 +74,10 @@ for T in (Float32, Float64, Complex64, Complex128)
     F = lufact(A)
     b = b/norm(b)
 
+    # Test optimality condition: residual should be non-increasing
+    x_gmres, c_gmres = gmres(A, b, log = true, restart = 3, maxiter = 10);
+    @fact all(diff(c_gmres[:resnorm]) .<= 0.0) --> true
+
     x_gmres, c_gmres = gmres(A, b, Pl=L, Pl=R, log=true)
     @fact c_gmres.isconverged --> true
     @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))
@@ -103,6 +107,10 @@ for T in (Float64, Complex128)
     F = lufact(A)
     b = b / norm(b)
 
+    # Test optimality condition: residual should be non-increasing
+    x_gmres, c_gmres = gmres(A, b, log = true, restart = 3, maxiter = 10);
+    @fact all(diff(c_gmres[:resnorm]) .<= 0.0) --> true
+
     x_gmres, c_gmres= gmres(A, b, Pl=L, Pr=R, log=true)
     @fact c_gmres.isconverged --> true
     @fact norm(A*x_gmres - b) --> less_than(√eps(real(one(T))))