reduced memory allocations to optimize cg.jl (#125)

Author: shivin9

.gitignore

@@ -1,3 +1,4 @@
 # mkdocs site
 /site
 /docs/build
+*.mem

src/cg.jl

@@ -51,7 +51,7 @@ function cg_method!(log::ConvergenceHistory, x, K, b;
         push!(log,:resnorm,resnorm)
         verbose && @printf("%3d\t%1.2e\n",iter,resnorm)
         resnorm < tol && break
-        z = solve(Pl,r)
+        solve!(z,Pl,r)
         oldγ = γ
         γ = dot(r, z)
         β = γ/oldγ

src/common.jl

@@ -44,8 +44,14 @@ end
 Solve `A\b` with a direct solver. When `A` is a function `A(b)` is dispatched instead.
 """
 solve(A::Function,b) = A(b)
+
 solve(A,b) = A\b
 
+solve!{T}(out::AbstractArray{T},A::Int,b::AbstractArray{T}) = scale!(out,b, 1/A)
+
+solve!{T}(out::AbstractArray{T},A,b::AbstractArray{T}) = A_ldiv_B!(out,A,b)
+solve!{T}(out::AbstractArray{T},A::Function,b::AbstractArray{T}) = copy!(out,A(b))
+
 """
     initrand!(v)
 Overwrite `v` with a random unitary vector of the same length.

src/gmres.jl

@@ -7,8 +7,8 @@ export gmres, gmres!
 gmres(A, b; kwargs...) = gmres!(zerox(A,b), A, b; kwargs...)
 
 function gmres!(x, A, b;
-    tol=sqrt(eps(typeof(real(b[1])))), restart::Int=min(20,length(b)),
-    maxiter::Int=restart, plot::Bool=false, log::Bool=false, kwargs...
+    tol = sqrt(eps(typeof(real(b[1])))), restart::Int=min(20,length(b)),
+    maxiter::Int = restart, plot::Bool=false, log::Bool=false, kwargs...
     )
     (plot & !log) && error("Can't plot when log keyword is false")
     history = ConvergenceHistory(partial=!log, restart=restart)
@@ -26,7 +26,7 @@ end
 
 #One Arnoldi iteration
 #Optionally takes a truncation parameter l
-function arnoldi!(K::KrylovSubspace, w; l=K.order)
+function arnoldi(K::KrylovSubspace, w; l=K.order)
     v = nextvec(K)
     w = copy(v)
     n = min(length(K.v), l)
@@ -86,7 +86,7 @@ function gmres_method!(log::ConvergenceHistory, x, A, b;
         for j = 1:restart
             nextiter!(log, mvps=1)
             #Calculate next orthonormal basis vector in the Krylov subspace
-            H[1:j+1, j] = arnoldi!(K, w)
+            H[1:j+1, j] = arnoldi(K, w)
 
             #Update QR factorization of H
             #The Q is stored as a series of Givens rotations in J