diff --git a/README.md b/README.md
index 8706541..7e43d01 100644
--- a/README.md
+++ b/README.md
@@ -50,13 +50,22 @@ In both cases, `ts` is an array of time snapshots for u, with `U[:,j] ≈ u(ts[j
 
 Apart from keyword arguments that affect the computation of Krylov subspaces (see the Arnoldi iteration section), you can also adjust the timestepping behavior using the arguments. By setting `adaptive=true`, the time step and Krylov subsapce size adaptation scheme of Niesen & Wright is used and the relative tolerance of which can be set using the keyword parameter `tol`. The `delta` and `gamma` parameter of the adaptation scheme can also be adjusted. The `tau` parameter adjusts the size of the timestep (and for `adaptive=true`, the initial timestep). By default, it is calculated using a heuristic formula by Niesen & Wright.
 
+### Support for matrix-free operators
+
+You can use any object as the "matrix" `A` as long as it implements the following linear operator interface:
+
+* `LinearAlgebra.mul!(y, A, x)` (for computing `y = A * x` in place).
+* `Base.size(A, dim)`
+* `LinearAlgebra.opnorm(A, p=Inf)`. If this is not implemented or the default implementation can be slow, you can manually pass in the operator norm (a rough estimate is fine) using the keyword argument `opnorm`.
+* `LinearAlgebra.ishermitian(A)`. If this is not implemented or the default implementation can be slow, you can manually pass in the value using the keyword argument `ishermitian`.
+
 ## Matrix-phi function `phi`
 
 ```julia
 phi(z,k[;cache]) -> [ϕ_0(z),ϕ_1(z),...,ϕ_k(z)]
 ```
 
-Compute ϕ function directly. `z` can both be a scalar or a Matrix. This is used by the caching versions of the ExpRK integrators to precompute the operators.
+Compute ϕ function directly. `z` can both be a scalar or a `AbstractMatrix` (note that unlike the previous functions, you *need* to use a concrete matrix). This is used by the caching versions of the ExpRK integrators to precompute the operators.
 
 Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used [2].
 
diff --git a/src/ExponentialUtilities.jl b/src/ExponentialUtilities.jl
index 1098fbe..a5c2a48 100644
--- a/src/ExponentialUtilities.jl
+++ b/src/ExponentialUtilities.jl
@@ -1,12 +1,15 @@
 module ExponentialUtilities
-using LinearAlgebra, SparseArrays
-using LinearAlgebra: exp!
+using LinearAlgebra, SparseArrays, Printf
+using LinearAlgebra: exp!, BlasInt
+using LinearAlgebra.LAPACK: stegr!
 
+include("utils.jl")
 include("phi.jl")
 include("arnoldi.jl")
 include("krylov_phiv.jl")
 include("krylov_phiv_adaptive.jl")
-include("krylov_expv.jl")
+include("StegrWork.jl")
+include("krylov_phiv_error_estimate.jl")
 
 export phi, phi!, KrylovSubspace, arnoldi, arnoldi!, lanczos!, ExpvCache, PhivCache,
     expv, expv!, phiv, phiv!, expv_timestep, expv_timestep!, phiv_timestep, phiv_timestep!,
diff --git a/src/stegr_cache.jl b/src/StegrWork.jl
similarity index 74%
rename from src/stegr_cache.jl
rename to src/StegrWork.jl
index f708955..de6a7d6 100644
--- a/src/stegr_cache.jl
+++ b/src/StegrWork.jl
@@ -1,5 +1,7 @@
-import LinearAlgebra.BLAS: @blasfunc
-import LinearAlgebra: BlasInt
+module Stegr
+using LinearAlgebra
+using LinearAlgebra.BLAS: @blasfunc
+using LinearAlgebra: BlasInt
 import LinearAlgebra.LAPACK: stegr!
 const liblapack = Base.liblapack_name
 
@@ -101,31 +103,5 @@ function StegrWork(::Type{T}, n::BlasInt,
     sw
 end
 
-mutable struct StegrCache{T,R<:Real} <: HermitianSubspaceCache{T}
-    v::Vector{T} # Subspace-propagated vector
-    w::Vector{T}
-    sw::StegrWork{R}
-    StegrCache(::Type{T}, n::Integer) where T = new{T,real(T)}(
-        Vector{T}(undef, n), Vector{T}(undef, n),
-        StegrWork(real(T), BlasInt(n)))
-end
-
-StegrCache(Ks::KrylovSubspace{B,T,U}) where {B,T,U} =
-    StegrCache(T, Ks.maxiter)
-
-"""
-    expT!(α, β, t, cache)
-
-Calculate the subspace exponential `exp(t*T)` for a tridiagonal
-subspace matrix `T` with `α` on the diagonal and `β` on the
-super-/subdiagonal, diagonalizing via `stegr!`.
-"""
-function expT!(α::AbstractVector{R}, β::AbstractVector{R}, t::Number,
-               cache::StegrCache{T,R}) where {T,R<:Real}
-    stegr!(α, β, cache.sw)
-    sel = 1:length(α)
-    @inbounds for i = sel
-        cache.w[i] = exp(t*cache.sw.w[i])*cache.sw.Z[1,i]
-    end
-    mul!(@view(cache.v[sel]), @view(cache.sw.Z[sel,sel]), @view(cache.w[sel]))
+export StegrWork
 end
diff --git a/src/arnoldi.jl b/src/arnoldi.jl
index 870d23a..d6b2807 100644
--- a/src/arnoldi.jl
+++ b/src/arnoldi.jl
@@ -52,14 +52,9 @@ function Base.show(io::IO, Ks::KrylovSubspace)
     println(IOContext(io, :limit => true), getH(Ks))
 end
 
-# Helper functions that returns the real part if that is what is
-# required (for Hermitian matrices), otherwise returns the value
-# as-is.
-coeff(::Type{T},α::T) where {T} = α
-coeff(::Type{U},α::T) where {U<:Real,T<:Complex} = real(α)
-
 #######################################
 # Arnoldi/Lanczos with custom IOP
+## High-level interface
 """
     arnoldi(A,b[;m,tol,opnorm,iop]) -> Ks
 
@@ -78,7 +73,7 @@ The default value of 0 indicates full Arnoldi. For symmetric/Hermitian `A`,
 
 Refer to `KrylovSubspace` for more information regarding the output.
 
-Happy-breakdown occurs whenver `norm(v_j) < tol * opnorm(A, Inf)`, in this case
+Happy-breakdown occurs whenver `norm(v_j) < tol * opnorm`, in this case
 the dimension of `Ks` is smaller than `m`.
 
 [^1]: Koskela, A. (2015). Approximating the matrix exponential of an
@@ -86,11 +81,12 @@ advection-diffusion operator using the incomplete orthogonalization method. In
 Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353).
 Springer, Cham.
 """
-function arnoldi(A, b; m=min(30, size(A, 1)), tol=1e-7, opnorm=LinearAlgebra.opnorm, iop=0, cache=nothing )
+function arnoldi(A, b; m=min(30, size(A, 1)), kwargs...)
     Ks = KrylovSubspace{eltype(b)}(length(b), m)
-    arnoldi!(Ks, A, b; m=m, tol=tol, opnorm=opnorm, iop=iop)
+    arnoldi!(Ks, A, b; m=m, kwargs...)
 end
 
+## Low-level interface
 """
     arnoldi_step!(j, iop, n, A, V, H)
 
@@ -110,6 +106,7 @@ function arnoldi_step!(j::Integer, iop::Integer, A,
     @. y /= beta
     beta
 end
+
 """
     arnoldi!(Ks,A,b[;tol,m,opnorm,iop]) -> Ks
 
@@ -117,8 +114,9 @@ Non-allocating version of `arnoldi`.
 """
 function arnoldi!(Ks::KrylovSubspace{B, T, U}, A, b::AbstractVector{T};
                   tol::Real=1e-7, m::Int=min(Ks.maxiter, size(A, 1)),
-                  opnorm=LinearAlgebra.opnorm, iop::Int=0, cache=nothing) where {B, T <: Number, U <: Number}
-    if ishermitian(A)
+                  ishermitian::Bool=LinearAlgebra.ishermitian(A),
+                  opnorm=LinearAlgebra.opnorm(A,Inf), iop::Int=0) where {B, T <: Number, U <: Number}
+    if ishermitian
         return lanczos!(Ks, A, b; tol=tol, m=m, opnorm=opnorm)
     end
     if m > Ks.maxiter
@@ -127,7 +125,8 @@ function arnoldi!(Ks::KrylovSubspace{B, T, U}, A, b::AbstractVector{T};
         Ks.m = m # might change if happy-breakdown occurs
     end
     V, H = getV(Ks), getH(Ks)
-    vtol = tol * opnorm(A, Inf)
+    # vtol = tol * opnorm
+    vtol = tol * (opnorm isa Function ? opnorm(A,Inf) : opnorm) # backward compatibility
     if iop == 0
         iop = m
     end
@@ -167,28 +166,6 @@ function lanczos_step!(j::Integer, A,
     β[j]
 end
 
-macro diagview(A,d::Integer=0)
-    s = d<=0 ? 1+abs(d) : :(m+$d)
-    quote
-        m = size($(esc(A)),1)
-        @view($(esc(A))[($s):m+1:end])
-    end
-end
-
-"""
-    realview(R, V)
-
-Returns a view of the real components of the complex vector `V`.
-"""
-realview(::Type{R}, V::AbstractVector{C}) where {R,C<:Complex} =
-    @view(reinterpret(R, V)[1:2:end])
-"""
-    realview(R, V)
-
-Returns a view of the real components of the real vector `V`.
-"""
-realview(::Type{R}, V::AbstractVector{R}) where {R} = V
-
 """
     lanczos!(Ks,A,b[;tol,m,opnorm]) -> Ks
 
@@ -197,14 +174,15 @@ Hermitian matrices.
 """
 function lanczos!(Ks::KrylovSubspace{B, T, U}, A, b::AbstractVector{T};
                   tol=1e-7, m=min(Ks.maxiter, size(A, 1)),
-                  opnorm=LinearAlgebra.opnorm, cache=nothing) where {B, T <: Number, U <: Number}
+                  opnorm=LinearAlgebra.opnorm(A,Inf)) where {B, T <: Number, U <: Number}
     if m > Ks.maxiter
         resize!(Ks, m)
     else
         Ks.m = m # might change if happy-breakdown occurs
     end
     V, H = getV(Ks), getH(Ks)
-    vtol = tol * opnorm(A, Inf)
+    # vtol = tol * opnorm
+    vtol = tol * (opnorm isa Function ? opnorm(A,Inf) : opnorm) # backward compatibility
     # Safe checks
     n = size(V, 1)
     @assert length(b) == size(A,1) == size(A,2) == n "Dimension mismatch"
diff --git a/src/krylov_phiv.jl b/src/krylov_phiv.jl
index 198ee5f..a53de76 100644
--- a/src/krylov_phiv.jl
+++ b/src/krylov_phiv.jl
@@ -32,32 +32,36 @@ employed by setting the kwarg `mode=:error_estimate`.
 
 Compute the expv product using a pre-constructed Krylov subspace.
 """
-function expv(t, A, b; m=min(30, size(A, 1)), tol=1e-7, rtol=√(tol),
-              opnorm=LinearAlgebra.opnorm, cache=nothing, iop=0,
-              mode=:happy_breakdown)
-    Ks = if mode == :happy_breakdown
-        arnoldi(A, b; m=m, tol=tol, opnorm=opnorm, iop=iop)
+function expv(t, A, b; mode=:happy_breakdown, kwargs...)
+    # Master dispatch function
+    if mode == :happy_breakdown
+        _expv_hb(t, A, b; kwargs...)
     elseif mode == :error_estimate
-        n = size(A,1)
-        T = promote_type(typeof(t), eltype(A), eltype(b))
-        U = ishermitian(A) ? real(T) : T
-        KrylovSubspace{T,U}(n, m)
+        _expv_ee(t, A, b; kwargs...)
     else
-        error("Unknown Krylov iteration termination mode, $(mode)")
+        throw(ArgumentError("Unknown Krylov iteration termination mode, $(mode)"))
     end
+end
+function _expv_hb(t, A, b; cache=nothing, kwargs_arnoldi...)
+    # Happy-breakdown mode: first construct Krylov subspace then expv!
+    Ks = arnoldi(A, b; kwargs_arnoldi...)
     w = similar(b)
-    if mode == :happy_breakdown
-        expv!(w, t, Ks; cache=cache)
-    elseif mode == :error_estimate
-        expv!(w, t, A, b,
-              Ks, get_subspace_cache(Ks);
-              atol=tol, rtol=rtol)
-    end
+    expv!(w, t, Ks; cache=cache)
+end
+function _expv_ee(t, A, b; m=min(30, size(A, 1)), tol=1e-7, rtol=√(tol),
+    ishermitian::Bool=LinearAlgebra.ishermitian(A))
+    # Error-estimate mode: construction of Krylov subspace and expv! at the same time
+    n = size(A,1)
+    T = promote_type(typeof(t), eltype(A), eltype(b))
+    U = ishermitian ? real(T) : T
+    Ks = KrylovSubspace{T,U}(n, m)
+    w = similar(b)
+    expv!(w, t, A, b, Ks, get_subspace_cache(Ks); atol=tol, rtol=rtol)
 end
-function expv(t, Ks::KrylovSubspace{B, T, U}; cache=nothing) where {B, T, U}
+function expv(t, Ks::KrylovSubspace{B, T, U}; kwargs...) where {B, T, U}
     n = size(getV(Ks), 1)
     w = Vector{T}(undef, n)
-    expv!(w, t, Ks; cache=cache)
+    expv!(w, t, Ks; kwargs...)
 end
 """
     expv!(w,t,Ks[;cache]) -> w
@@ -138,17 +142,15 @@ Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.
 the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.
 Formula (10).
 """
-function phiv(t, A, b, k; m=min(30, size(A, 1)), tol=1e-7, opnorm=LinearAlgebra.opnorm, iop=0,
-              cache=nothing, correct=false, errest=false)
-    Ks = arnoldi(A, b; m=m, tol=tol, opnorm=opnorm, iop=iop)
+function phiv(t, A, b, k; cache=nothing, correct=false, errest=false, kwargs_arnoldi...)
+    Ks = arnoldi(A, b; kwargs_arnoldi...)
     w = Matrix{eltype(b)}(undef, length(b), k+1)
     phiv!(w, t, Ks, k; cache=cache, correct=correct, errest=errest)
 end
-function phiv(t, Ks::KrylovSubspace{B, T, U}, k; cache=nothing, correct=false,
-              errest=false) where {B, T, U}
+function phiv(t, Ks::KrylovSubspace{B, T, U}, k; kwargs...) where {B, T, U}
     n = size(getV(Ks), 1)
     w = Matrix{T}(undef, n, k+1)
-    phiv!(w, t, Ks, k; cache=cache, correct=correct, errest=errest)
+    phiv!(w, t, Ks, k; kwargs...)
 end
 """
     phiv!(w,t,Ks,k[;cache,correct,errest]) -> w[,errest]
diff --git a/src/krylov_phiv_adaptive.jl b/src/krylov_phiv_adaptive.jl
index a6fc769..c34e62f 100644
--- a/src/krylov_phiv_adaptive.jl
+++ b/src/krylov_phiv_adaptive.jl
@@ -100,17 +100,20 @@ function phiv_timestep!(u::AbstractVector{T}, t::tType, A, B::AbstractMatrix{T};
     return u
 end
 function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractMatrix{T}; tau::Real=0.0,
-                        m::Int=min(10, size(A, 1)), tol::Real=1e-7, opnorm=LinearAlgebra.opnorm, iop::Int=0,
+                        m::Int=min(10, size(A, 1)), tol::Real=1e-7, opnorm=LinearAlgebra.opnorm(A,Inf), iop::Int=0,
                         correct::Bool=false, caches=nothing, adaptive=false, delta::Real=1.2,
+                        ishermitian::Bool=LinearAlgebra.ishermitian(A),
                         gamma::Real=0.8, NA::Int=0, verbose=false) where {T <: Number, tType <: Real}
     # Choose initial timestep
-    abstol = tol * opnorm(A, Inf)
+    if opnorm isa Function
+        opnorm = opnorm(A,Inf) # backward compatibility
+    end
+    abstol = tol * opnorm
     verbose && println("Absolute tolerance: $abstol")
     if iszero(tau)
-        Anorm = opnorm(A, Inf)
         b0norm = norm(@view(B[:, 1]), Inf)
-        tau = 10/Anorm * (abstol * ((m+1)/ℯ)^(m+1) * sqrt(2*pi*(m+1)) /
-                          (4*Anorm*b0norm))^(1/m)
+        tau = 10/opnorm * (abstol * ((m+1)/ℯ)^(m+1) * sqrt(2*pi*(m+1)) /
+                          (4*opnorm*b0norm))^(1/m)
         verbose && println("Initial time step unspecified, chosen to be $tau")
     end
     # Initialization
@@ -135,7 +138,7 @@ function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractM
     copyto!(u, @view(B[:, 1])) # u(0) = b0
     coeffs = ones(tType, p);
     if adaptive # initialization step for the adaptive scheme
-        if ishermitian(A)
+        if ishermitian
             iop = 2 # does not have an effect on arnoldi!, just for flops estimation
         end
         if iszero(NA)
@@ -176,7 +179,7 @@ function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractM
             while omega > delta # inner loop of Algorithm 3
                 m_new, tau_new, q, kappa = _phiv_timestep_adapt(
                     m, tau, epsilon, m_old, tau_old, epsilon_old, q, kappa,
-                    gamma, omega, maxtau, n, p, NA, iop, opnorm(getH(Ks), 1), verbose)
+                    gamma, omega, maxtau, n, p, NA, iop, LinearAlgebra.opnorm(getH(Ks), 1), verbose)
                 m, m_old = m_new, m
                 tau, tau_old = tau_new, tau
                 # Compute ϕp(tau*A)wp using the new parameters
diff --git a/src/krylov_expv.jl b/src/krylov_phiv_error_estimate.jl
similarity index 63%
rename from src/krylov_expv.jl
rename to src/krylov_phiv_error_estimate.jl
index bba7822..4be63c6 100644
--- a/src/krylov_expv.jl
+++ b/src/krylov_phiv_error_estimate.jl
@@ -1,15 +1,45 @@
-using Printf
+# Alternative Krylov phiv methods using error estimates of Saad to automatically
+# terminate Arnoldi/Lanczos iterations.
+# Currently only expv for Lanczos is implemented.
 
+########################################
+# Cache types
 abstract type SubspaceCache{T} end
 abstract type HermitianSubspaceCache{T} <: SubspaceCache{T} end
 
-include("stegr_cache.jl")
+mutable struct StegrCache{T,R<:Real} <: HermitianSubspaceCache{T}
+    v::Vector{T} # Subspace-propagated vector
+    w::Vector{T}
+    sw::Stegr.StegrWork{R}
+    StegrCache(::Type{T}, n::Integer) where T = new{T,real(T)}(
+        Vector{T}(undef, n), Vector{T}(undef, n),
+        Stegr.StegrWork(real(T), BlasInt(n)))
+end
+
+"""
+    expT!(α, β, t, cache)
+
+Calculate the subspace exponential `exp(t*T)` for a tridiagonal
+subspace matrix `T` with `α` on the diagonal and `β` on the
+super-/subdiagonal, diagonalizing via `stegr!`.
+"""
+function expT!(α::AbstractVector{R}, β::AbstractVector{R}, t::Number,
+               cache::StegrCache{T,R}) where {T,R<:Real}
+    stegr!(α, β, cache.sw)
+    sel = 1:length(α)
+    @inbounds for i = sel
+        cache.w[i] = exp(t*cache.sw.w[i])*cache.sw.Z[1,i]
+    end
+    mul!(@view(cache.v[sel]), @view(cache.sw.Z[sel,sel]), @view(cache.w[sel]))
+end
 
 get_subspace_cache(Ks::KrylovSubspace{B,T,U}) where {B,T,U<:Complex} =
     error("Subspace exponential caches not yet available for non-Hermitian matrices.")
 get_subspace_cache(Ks::KrylovSubspace{B,T,U}) where {B,T,U<:Real} =
-    StegrCache(Ks)
+    StegrCache(T, Ks.maxiter)
 
+########################################
+# Phiv with error estimate as termination condition
 """
     expv!(w, t, A, b, Ks, cache)
 
diff --git a/src/phi.jl b/src/phi.jl
index 72147b1..a01af46 100644
--- a/src/phi.jl
+++ b/src/phi.jl
@@ -48,9 +48,9 @@ formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software
 package for computing matrix exponentials. ACM Transactions on Mathematical
 Software (TOMS), 24(1), 130-156. Theorem 1).
 """
-function phiv_dense(A, v, k; cache=nothing)
+function phiv_dense(A, v, k; kwargs...)
     w = Matrix{eltype(A)}(undef, length(v), k+1)
-    phiv_dense!(w, A, v, k; cache=cache)
+    phiv_dense!(w, A, v, k; kwargs...)
 end
 """
     phiv_dense!(w,A,v,k[;cache]) -> w
@@ -100,14 +100,14 @@ Calls `phiv_dense` on each of the basis vectors to obtain the answer. If A
 is `Diagonal`, instead calls the scalar `phi` on each diagonal element and the
 return values are also `Diagonal`s
 """
-function phi(A::AbstractMatrix{T}, k; caches=nothing) where {T <: Number}
+function phi(A::AbstractMatrix{T}, k; kwargs...) where {T <: Number}
     m = size(A, 1)
     if A isa Diagonal
         out = [similar(A) for i = 1:k+1]
     else
         out = [Matrix{T}(undef, m, m) for i = 1:k+1]
     end
-    phi!(out, A, k; caches=caches)
+    phi!(out, A, k; kwargs...)
 end
 """
     phi!(out,A,k[;caches]) -> out
diff --git a/src/utils.jl b/src/utils.jl
new file mode 100644
index 0000000..a6a723f
--- /dev/null
+++ b/src/utils.jl
@@ -0,0 +1,29 @@
+# Utility functions
+
+"""
+    coeff(::Type,α)
+
+Helper functions that returns the real part if that is what is
+required (for Hermitian matrices), otherwise returns the value
+as-is.
+"""
+coeff(::Type{T},α::T) where {T} = α
+coeff(::Type{U},α::T) where {U<:Real,T<:Complex} = real(α)
+
+"""
+    @diagview(A,d) -> view of the `d`th diagonal of `A`.
+"""
+macro diagview(A,d::Integer=0)
+    s = d<=0 ? 1+abs(d) : :(m+$d)
+    quote
+        m = size($(esc(A)),1)
+        @view($(esc(A))[($s):m+1:end])
+    end
+end
+
+"""
+    realview(::Type, V) -> real view of `V`
+"""
+realview(::Type{R}, V::AbstractVector{C}) where {R,C<:Complex} =
+    @view(reinterpret(R, V)[1:2:end])
+realview(::Type{R}, V::AbstractVector{R}) where {R} = V