From f8a53a3f51fe452bc9e89c96ee2356aa4ba9635e Mon Sep 17 00:00:00 2001 From: Leo Date: Sat, 28 Jul 2018 02:53:17 +0800 Subject: [PATCH 1/3] test linear operator --- src/expmv.jl | 14 ++++++++------ test/runtests.jl | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 6 deletions(-) diff --git a/src/expmv.jl b/src/expmv.jl index b7d9875..1ca6705 100644 --- a/src/expmv.jl +++ b/src/expmv.jl @@ -1,7 +1,7 @@ export expmv, expmv! """ - 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 +13,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::Real=0) 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,14 +23,17 @@ 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::Real=0) = 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) + tol::Real=1e-7, m::Int=min(30,size(A,1)), norm=Base.norm, anorm::Real=0) if size(vec,1) != size(A,2) error("dimension mismatch") end + if anorm == 0 # set default anorm + anorm = hasmethod(norm, Tuple{typeof(A), typeof(Inf)}) ? norm(A, Inf) : 1.0 + end # safety factors gamma = 0.9 @@ -39,7 +42,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/test/runtests.jl b/test/runtests.jl index 3b377f0..c7ba513 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,39 @@ function test_expmv3() return e1, e2 end +struct LinearOp + m +end + +Base.A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x) +Base.size(lo::LinearOp) = size(lo.m) +Base.size(lo::LinearOp, i::Int) = size(lo.m, i) +Base.eltype(lo::LinearOp) = eltype(lo.m) + +function test_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_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") From 9d6f284f9b633f9ebeb7751b254c3fb0ee88bc8b Mon Sep 17 00:00:00 2001 From: Leo Date: Wed, 8 Aug 2018 00:21:28 +0800 Subject: [PATCH 2/3] update anorm, throw warning. --- src/expmv.jl | 29 +++++++++++++++++++----- src/phimv.jl | 12 +++++----- test/runtests.jl | 57 ++++++++++++++++++++++++++++++++++++++---------- 3 files changed, 76 insertions(+), 22 deletions(-) diff --git a/src/expmv.jl b/src/expmv.jl index 1ca6705..b42f68c 100644 --- a/src/expmv.jl +++ b/src/expmv.jl @@ -1,5 +1,20 @@ export expmv, expmv! +function default_anorm(A) + try + norm(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], [anorm]) @@ -13,7 +28,7 @@ function expmv{T}( t::Number, A, vec::Vector{T}; tol::Real=1e-7, m::Int=min(30, size(A, 1)), - norm=Base.norm, anorm::Real=0) + norm=Base.norm, anorm=nothing) result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec)) expmv!(t, A, result; tol=tol, m=m, norm=norm, anorm=anorm) return result @@ -23,16 +38,18 @@ expmv!{T}( t::Number, A, vec::Vector{T}; tol::Real=1e-7, m::Int=min(30,size(A,1)), - norm=Base.norm, anorm::Real=0) = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm, anorm=anorm) + norm=Base.norm, anorm=nothing) = 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, anorm::Real=0) +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=nothing) if size(vec,1) != size(A,2) error("dimension mismatch") end - if anorm == 0 # set default anorm - anorm = hasmethod(norm, Tuple{typeof(A), typeof(Inf)}) ? norm(A, Inf) : 1.0 + + # set default anorm + if anorm === nothing + anorm = default_anorm(A) end # safety factors diff --git a/src/phimv.jl b/src/phimv.jl index ed931d7..4dee4d7 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=nothing) 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=nothing) = 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=nothing) if size(vec, 1) != size(A, 2) error("dimension mismatch") @@ -61,7 +61,9 @@ 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) + if anorm === nothing + anorm = default_anorm(A) + end rndoff = anorm*eps() # estimate first time-step and round to two significant digits diff --git a/test/runtests.jl b/test/runtests.jl index c7ba513..227b453 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,16 @@ using Expokit using Base.Test +struct LinearOp + m +end + +Base.A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x) +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,16 +50,7 @@ function test_expmv3() return e1, e2 end -struct LinearOp - m -end - -Base.A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x) -Base.size(lo::LinearOp) = size(lo.m) -Base.size(lo::LinearOp, i::Int) = size(lo.m, i) -Base.eltype(lo::LinearOp) = eltype(lo.m) - -function test_linop(n::Int) +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)[:] @@ -69,7 +70,7 @@ function test_linop(n::Int) end println("testing linear operator n=1000 (first expmv, then expm)") -res, t1, t2, t3 = test_linop(1000) +res, t1, t2, t3 = test_expmv_linop(1000) println("residuum: $res\n") @test res < 1e-6 @@ -280,3 +281,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 From 8254cb47b9163a034c0b6c35e785ec58ca366396 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 9 Aug 2018 23:26:54 +0800 Subject: [PATCH 3/3] minor change --- src/expmv.jl | 13 ++++--------- src/phimv.jl | 11 ++++------- test/runtests.jl | 8 +++++++- 3 files changed, 15 insertions(+), 17 deletions(-) diff --git a/src/expmv.jl b/src/expmv.jl index b42f68c..b48a42e 100644 --- a/src/expmv.jl +++ b/src/expmv.jl @@ -2,7 +2,7 @@ export expmv, expmv! function default_anorm(A) try - norm(A, Inf) + 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`. @@ -28,7 +28,7 @@ function expmv{T}( t::Number, A, vec::Vector{T}; tol::Real=1e-7, m::Int=min(30, size(A, 1)), - norm=Base.norm, anorm=nothing) + 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, anorm=anorm) return result @@ -38,20 +38,15 @@ expmv!{T}( t::Number, A, vec::Vector{T}; tol::Real=1e-7, m::Int=min(30,size(A,1)), - norm=Base.norm, anorm=nothing) = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm, anorm=anorm) + 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, anorm=nothing) + 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") end - # set default anorm - if anorm === nothing - anorm = default_anorm(A) - end - # safety factors gamma = 0.9 delta = 1.2 diff --git a/src/phimv.jl b/src/phimv.jl index 4dee4d7..6b68cf2 100644 --- a/src/phimv.jl +++ b/src/phimv.jl @@ -35,7 +35,7 @@ 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, anorm=nothing) + 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, anorm=anorm) return result @@ -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, anorm=nothing) = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm, anorm=anorm) + 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, anorm=nothing) + 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,9 +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 - if anorm === nothing - anorm = default_anorm(A) - end rndoff = anorm*eps() # estimate first time-step and round to two significant digits @@ -96,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 227b453..4478186 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,7 +5,13 @@ struct LinearOp m end -Base.A_mul_B!(y, lo::LinearOp, x) = A_mul_B!(y, lo.m, x) +@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: *