Skip to content

Commit

Permalink
Add anorm keyword to support to linear maps (#30)
Browse files Browse the repository at this point in the history
* test linear operator

* update norm keyword, throw warning.
  • Loading branch information
GiggleLiu authored and acroy committed Aug 9, 2018
1 parent 3346967 commit b5e020e
Showing 3 changed files with 100 additions and 13 deletions.
28 changes: 21 additions & 7 deletions src/expmv.jl
Original file line number Diff line number Diff line change
@@ -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,20 +28,20 @@ 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

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
11 changes: 5 additions & 6 deletions src/phimv.jl
Original file line number Diff line number Diff line change
@@ -35,20 +35,20 @@ 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

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)
74 changes: 74 additions & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit b5e020e

Please sign in to comment.