Fix GMRES memory usage (#131)

* Add improved versions of GMRES and CG

* Add preconditioners

* Add Hessenberg type

* Hessenberg solver

* Fix incorrect comment

* Benchmark Givens rotation + UpperTriangular vs \

* Rename to Hessenberg

* Use UnsafeVectorView when solving the small system

* Clean up some

* Add Identity preconditioner type

* Improve GMRES with specialized functions for preconditioners, unsafe views, and conditional logging

* Don't run by default

* Some default values for gmres bench

* Reduce the number of allocs slightly by solving the least-squares-problem with the Hessenberg matrix in-place

* Fix for complex numbers

* Move from Vector{Vector{T}} to Matrix{T} using views and make orthogonalization routine reusable

* Make the comparison fair and the matrix sparser to ~10 nonzeros per row

* Fix a bug for complex numbers

* Temporarily remove tests with matrices as preconditioner, since we don't have A_ldiv_B for them at the moment

* Add tests for various orthogonalization methods

* Fix a problem for complex numbers again

* Export Identity

* Add tests to runtests.jl

* Fix compatibility with Julia 0.5

* Override the current gmres method and move solve! to A_ldiv_B!

Author: haampie

benchmark/benchmark-hessenberg.jl

@@ -0,0 +1,53 @@
+module HessenbergBench
+
+using IterativeSolvers
+using BenchmarkTools
+
+function backslash_versus_givens(; n = 1_000, ms = 10 : 10 : 100)
+    results = BenchmarkGroup()
+    results["backslash"] = BenchmarkGroup()
+    results["givens_qr"] = BenchmarkGroup()
+
+    # Example matrix
+    A = Tridiagonal(fill(-1.5, n - 1), fill(3.0, n), fill(-1.0, n - 1))
+
+    # Different sizes of the Hessenberg matrix
+    for m = ms
+        println(m)
+
+        H = zeros(m + 1, m)
+        
+        # Create an orthonormal basis for the Krylov subspace
+        V = rand(n, m + 1)
+        V[:, 1] /= norm(V[:, 1])
+        
+        for i = 1 : m
+            # Expand
+            V[:, i + 1] = A * V[:, i]
+            
+            # Orthogonalize
+            H[1 : i, i] = view(V, :, 1 : i)' * V[:, i + 1]
+            V[:, i + 1] -= view(V, :, 1 : i) * H[1 : i, i]
+            
+            # Re-orthogonalize
+            update = view(V, :, 1 : i)' * V[:, i + 1]
+            H[1 : i, i] += update
+            V[:, i + 1] -= view(V, :, 1 : i) * update
+            
+            # Normalize
+            H[i + 1, i] = norm(V[:, i + 1])
+            V[:, i + 1] /= H[i + 1, i]
+        end
+
+        # Fist standard basis vector as rhs
+        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["backslash"][m] = @benchmark $H \ $rhs
+    end
+
+    return results
+end
+
+end

benchmark/benchmark-linear-systems.jl

@@ -41,4 +41,27 @@ function cg(; n = 1_000_000, tol = 1e-6, maxiter::Int = 200)
     @benchmark IterativeSolvers.cg!(x0, $A, $b, Pl = $P, maxiter = $maxiter, tol = $tol, log = false) setup=(x0 = copy($initial))
 end
 
-end
+function indefinite(n)
+    # Generate an indefinite "hard" matrix
+    srand(1)
+    A = speye(n) + sprand(n, n, 5.0 / n)
+    A = (A + A') / 2
+    x = ones(n)
+    b = A * x
+
+    A, b
+end
+
+function gmres(; n = 100_000, tol = 1e-5, restart::Int = 15, maxiter::Int = 210)
+    A, b = indefinite(n)
+
+    println("Matrix of size ", n, " with ~", nnz(A) / n, " nonzeros per row")
+    println("Tolerance = ", tol, "; restart = ", restart, "; max #iterations = ", maxiter)
+    
+    impr = @benchmark IterativeSolvers.improved_gmres($A, $b, tol = $tol, restart = $restart, maxiter = $maxiter, log = false)
+    old = @benchmark IterativeSolvers.gmres($A, $b, tol = $tol, restart = $restart, maxiter = $maxiter, log = false)
+
+    impr, old
+end
+
+end

src/IterativeSolvers.jl

@@ -10,11 +10,13 @@ using SugarBLAS
 using Compat
 
 include("common.jl")
+include("orthogonalize.jl")
 include("krylov.jl")
 include("history.jl")
 
 #Specialized factorizations
 include("factorization.jl")
+include("hessenberg.jl")
 
 #Linear solvers
 include("stationary.jl")

src/common.jl

@@ -1,9 +1,9 @@
 import  Base: eltype, length, ndims, real, size, *,
         A_mul_B!, Ac_mul_B, Ac_mul_B!
 
-using   LinearMaps
+using   LinearMaps, Compat
 
-export  A_mul_B
+export  A_mul_B, Identity
 
 #### Type-handling
 """
@@ -74,3 +74,10 @@ type PosSemidefException <: Exception
     msg :: AbstractString
     PosSemidefException(msg::AbstractString="Matrix was not positive semidefinite") = new(msg)
 end
+
+# Identity preconditioner
+immutable Identity end
+
+Base.:\(::Identity, x) = copy(x)
+Base.A_ldiv_B!(::Identity, x) = x
+Base.A_ldiv_B!(y, ::Identity, x) = copy!(y, x)