Fix KrylovSubspace constructors to work for custom operators. Fixes #48

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}[])

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)

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