diff --git a/src/expmv.jl b/src/expmv.jl
index b7d9875..b48a42e 100644
--- a/src/expmv.jl
+++ b/src/expmv.jl
@@ -1,7 +1,22 @@
 export expmv, expmv!
 
+function default_anorm(A)
+    try
+        Compat.opnorm(A, Inf)
+    catch err
+        if err isa MethodError
+            warn("opnorm($(typeof(A)), Inf) is not defined, fall back to using `anorm = 1.0`.
+To suppress this warning, please specify the anorm parameter manually.")
+            1.0
+        else
+            throw(err)
+        end
+    end
+end
+
+
 """
-    expmv{T}(t, A, vec; [tol], [m], [norm])
+    expmv{T}(t, A, vec; [tol], [m], [norm], [anorm])
 
 Calculate matrix exponential acting on some vector, ``w = e^{tA}v``,
 using the Krylov subspace approximation.
@@ -13,9 +28,9 @@ function expmv{T}( t::Number,
                    A, vec::Vector{T};
                    tol::Real=1e-7,
                    m::Int=min(30, size(A, 1)),
-                   norm=Base.norm)
+                   norm=Base.norm, anorm=default_anorm(A))
     result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec))
-    expmv!(t, A, result; tol=tol, m=m, norm=norm)
+    expmv!(t, A, result; tol=tol, m=m, norm=norm, anorm=anorm)
     return result
 end
 
@@ -23,10 +38,10 @@ expmv!{T}( t::Number,
            A, vec::Vector{T};
            tol::Real=1e-7,
            m::Int=min(30,size(A,1)),
-           norm=Base.norm) = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm)
+           norm=Base.norm, anorm=default_anorm(A)) = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm, anorm=anorm)
 
-function expmv!{T}( w::Vector{T}, t::Number, A, vec::Vector{T};
-                    tol::Real=1e-7, m::Int=min(30,size(A,1)), norm=Base.norm)
+function expmv!{T}(w::Vector{T}, t::Number, A, vec::Vector{T};
+                   tol::Real=1e-7, m::Int=min(30,size(A,1)), norm=Base.norm, anorm=default_anorm(A))
 
     if size(vec,1) != size(A,2)
         error("dimension mismatch")
@@ -39,7 +54,6 @@ function expmv!{T}( w::Vector{T}, t::Number, A, vec::Vector{T};
     btol = 1e-7     # tolerance for "happy-breakdown"
     maxiter = 10    # max number of time-step refinements
 
-    anorm = norm(A, Inf)
     rndoff= anorm*eps()
 
     # estimate first time-step and round to two significant digits
diff --git a/src/phimv.jl b/src/phimv.jl
index ed931d7..6b68cf2 100644
--- a/src/phimv.jl
+++ b/src/phimv.jl
@@ -35,9 +35,9 @@ function phimv{T}( t::Number,
                    A, u::Vector{T}, vec::Vector{T};
                    tol::Real=1e-7,
                    m::Int=min(30, size(A, 1)),
-                   norm=Base.norm)
+                   norm=Base.norm, anorm=default_anorm(A))
     result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec))
-    phimv!(t, A, u, result; tol=tol, m=m, norm=norm)
+    phimv!(t, A, u, result; tol=tol, m=m, norm=norm, anorm=anorm)
     return result
 end
 
@@ -45,10 +45,10 @@ phimv!{T}( t::Number,
            A, u::Vector{T}, vec::Vector{T};
            tol::Real=1e-7,
            m::Int=min(30, size(A, 1)),
-           norm=Base.norm) = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm)
+           norm=Base.norm, anorm=default_anorm(A)) = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm, anorm=anorm)
 
 function phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
-                    tol::Real=1e-7, m::Int=min(30, size(A, 1)), norm=Base.norm)
+                   tol::Real=1e-7, m::Int=min(30, size(A, 1)), norm=Base.norm, anorm=default_anorm(A))
 
     if size(vec, 1) != size(A, 2)
         error("dimension mismatch")
@@ -61,7 +61,6 @@ function phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
     btol = 1e-7     # tolerance for "happy-breakdown"
     maxiter = 10    # max number of time-step refinements
 
-    anorm = norm(A, Inf)
     rndoff = anorm*eps()
 
     # estimate first time-step and round to two significant digits
@@ -94,7 +93,7 @@ function phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
         scale!(copy!(vm[1], A*w+u), 1/beta)
 
         for j = 1:m
-            A_mul_B!(p, A, vm[j])
+            Base.A_mul_B!(p, A, vm[j])
 
             for i = 1:j
                 hm[i, j] = dot(vm[i], p)
diff --git a/test/runtests.jl b/test/runtests.jl
index 3b377f0..4478186 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,22 @@
 using Expokit
 using Base.Test
 
+struct LinearOp
+    m
+end
+
+@static if VERSION < v"0.7-"
+    Base.A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x)
+else
+    import LinearAlgebra: mul!, A_mul_B!
+    mul!(y, lo::LinearOp, x) = mul!(y, lo.m, x)
+    A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x)
+end
+Base.size(lo::LinearOp, i::Int) = size(lo.m, i)
+Base.eltype(lo::LinearOp) = eltype(lo.m)
+import Base: *
+*(lo::LinearOp, v::Vector) = Base.A_mul_B!(similar(v), lo, v)
+
 function test_expmv(n::Int)
 
     A = sprand(n,n,0.4)
@@ -40,6 +56,30 @@ function test_expmv3()
     return e1, e2
 end
 
+function test_expmv_linop(n::Int)
+    A = LinearOp(sprand(n,n,0.2) + 1im*sprand(n,n,0.2))
+    v = eye(n,1)[:]+0im*eye(n,1)[:]
+
+    tic()
+    w1 = expmv(1.0, A, v, anorm=norm(A.m, Inf))
+    t1 = toc()
+
+    tic()
+    w2 = expmv(1.0, A, v)
+    t2 = toc()
+
+    tic()
+    w3 = expm(full(A.m))*v
+    t3 = toc()
+
+    return max(norm(w1-w2)/norm(w2), norm(w1-w3)/norm(w3)), t1, t2, t3
+end
+
+println("testing linear operator n=1000 (first expmv, then expm)")
+res, t1, t2, t3 = test_expmv_linop(1000)
+println("residuum: $res\n")
+@test res < 1e-6
+
 println("testing real n=100 (first expmv, then expm)")
 res, t1, t2 = test_expmv(100)
 println("residuum: $res\n")
@@ -247,3 +287,37 @@ println("testing complex n=1000 (first phimv, then expm)")
 res, t1, t2 = test_phimv2(1000)
 println("residuum: $res\n")
 @test res < 1e-6
+
+function test_phimv_linop(n::Int)
+    p = 0.1
+    found = false
+    A = LinearOp(sprand(n,n,p))
+    u = rand(n)
+    x = similar(u)
+    while !found
+        try
+            copy!(x, A.m\u)
+            found = true
+        catch
+            A = LinearOp(sprand(n,n,p))
+        end
+    end
+    vec = eye(n, 1)[:]
+
+    w1 = phimv(1.0, A, u, vec) # warmup
+    tic()
+    w1 = phimv(1.0, A, u, vec)
+    t1 = toc()
+
+    w2 = expm(full(A.m))*(vec+x)-x # warmup
+    tic()
+    w2 = expm(full(A.m))*(vec+x)-x
+    t2 = toc()
+
+    return norm(w1-w2)/norm(w2), t1, t2
+end
+
+println("testing real n=1000 (first phimv, then expm), the linear operator version.")
+res, t1, t2 = test_phimv_linop(1000)
+println("residuum: $res\n")
+@test res < 1e-6