From b6472ec5e60f66aeb1236400546e15aa0c672499 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreasnoackjensen@gmail.com>
Date: Tue, 15 Sep 2015 11:19:26 +0200
Subject: [PATCH] Fix KrylovSubspace constructors to work for custom operators.
 Fixes #48

---
 src/krylov.jl    |  9 +++++----
 src/lanczos.jl   |  2 +-
 test/runtests.jl | 17 +++++++++++++++++
 3 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/src/krylov.jl b/src/krylov.jl
index f421b853..bf58db8a 100644
--- a/src/krylov.jl
+++ b/src/krylov.jl
@@ -17,12 +17,13 @@ type KrylovSubspace{T, OpT}
     mvps::Int  #Count of matrix-vector products
 end
 
-KrylovSubspace(A, n::Int, order::Int, T::Type)=KrylovSubspace{T,typeof(A)}(A, n, order, Vector{T}[], 0)
+KrylovSubspace(A, n::Int, order::Int, T::Type) = KrylovSubspace{T,typeof(A)}(A, n, order, Array{T,1}[], 0)
 
-KrylovSubspace{T}(A, n::Int, order::Int, v::Vector{Vector{T}}=Vector{T}[])=KrylovSubspace{T,typeof(A)}(A, n, order, v, 0)
+KrylovSubspace{T}(A, n::Int, order::Int, v::Vector{Vector{T}}) = KrylovSubspace{T,typeof(A)}(A, n, order, v, 0)
 
-KrylovSubspace{T}(A::AbstractMatrix{T}, order::Int, v::Vector{Vector{T}}=Vector{T}[])=
-    KrylovSubspace{T,typeof(A)}(A, size(A,2), order, v, 0)
+KrylovSubspace(A, n::Int, order::Int) = KrylovSubspace{eltype(A),typeof(A)}(A, n, order, Array{eltype(A),1}[], 0)
+
+KrylovSubspace(A, order::Int) = KrylovSubspace{eltype(A),typeof(A)}(A, size(A,2), order, Array{eltype(A),1}[], 0)
 
 #Reset an existing KrylovSubspace
 function KrylovSubspace{T}(A::KrylovSubspace{T}, order::Int=size(A,2), v::Vector{Vector{T}}=Vector{T}[])
diff --git a/src/lanczos.jl b/src/lanczos.jl
index afa9d4be..0dc9b1a2 100644
--- a/src/lanczos.jl
+++ b/src/lanczos.jl
@@ -19,7 +19,7 @@ function lanczos!{T}(K::KrylovSubspace{T})
 end
 
 function eigvals_lanczos(A, neigs::Int=size(A,1); tol::Real=size(A,1)^3*eps(), maxiter::Int=size(A,1))
-    K = KrylovSubspace(A, 2) #In Lanczos, only remember the last two vectors
+    K = KrylovSubspace(A, size(A, 1), 2) #In Lanczos, only remember the last two vectors
     initrand!(K)
     resnorms = zeros(maxiter)
     e1 = eigvals(lanczos!(K), 1:neigs)
diff --git a/test/runtests.jl b/test/runtests.jl
index 2499711a..e5d859bf 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -183,6 +183,16 @@ end
 
 #Lanczos methods
 
+type MyOp{T}
+    buf::Matrix{T}
+end
+
+import Base: *, size, eltype
+
+*(A::MyOp, x::AbstractVector) = A.buf*x
+size(A::MyOp, i) = size(A.buf, i)
+eltype(A::MyOp) = eltype(A.buf)
+
 facts("eigvals_lanczos") do
 for T in (Float32, Float64)
     context("Matrix{$T}") do
@@ -194,6 +204,13 @@ for T in (Float32, Float64)
     T==Float64 && @fact c_lanczos.isconverged --> true #XXX Lanczos needs to be made more robust for Float32
     @fact norm(v - eval_lanczos) --> less_than(√eps(T))
     end
+
+    context("Op{$T}") do
+        A = MyOp(convert(Matrix{T}, randn(5,5)) |> t -> t + t')
+        v = eigvals(Symmetric(A.buf))
+        @fact c_lanczos.isconverged --> true
+        @fact norm(v - eval_lanczos) --> less_than(√eps(T))
+    end
 end
 end