Skip to content

Commit

Permalink
Upgrade to Julia 0.7/1.0 (#34)
Browse files Browse the repository at this point in the history
* Upgrade to julia 0.7/1.0 (eg use mul! instead of A_mul_B!)
* drop support for julia 0.6
  • Loading branch information
acroy authored Aug 13, 2018
1 parent 1e462ba commit 4f91a8a
Show file tree
Hide file tree
Showing 9 changed files with 148 additions and 172 deletions.
1 change: 0 additions & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@ os:
- linux
- osx
julia:
- 0.6
- 0.7
- nightly
matrix:
Expand Down
4 changes: 2 additions & 2 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ Pkg.add("Expokit")
w = expmv!{T}( w::Vector{T}, t::Number, A, v::Vector{T}; kwargs...)
```
The function `expmv!` calculates `w = exp(t*A)*v`, where `A` is a
matrix or any type that supports `size` and `A_mul_B!` and `v` a dense vector by using Krylov subspace projections. The result is
matrix or any type that supports `size`, `eltype` and `mul!` and `v` is a dense vector by using Krylov subspace projections. The result is
stored in `w`.

The following keywords are supported
Expand All @@ -43,7 +43,7 @@ w = expmv{T}( t::Number, A, v::Vector{T}; kwargs...)
w = phimv!{T}( w::Vector{T}, t::Number, A, u::Vector{T}, v::Vector{T}; kwargs...)
```
The function `phimv!` calculates `w = e^{tA}v + t φ(t A) u` with `φ(z) = (exp(z)-1)/z`, where `A` is a
matrix or any type that supports `size` and `A_mul_B!`, `u` and `v` are dense vectors by using Krylov subspace projections. The result is stored in `w`. The supported keywords are the same as for `expmv!`.
matrix or any type that supports `size`, `eltype` and `mul!`, `u` and `v` are dense vectors by using Krylov subspace projections. The result is stored in `w`. The supported keywords are the same as for `expmv!`.

## chbv

Expand Down
3 changes: 1 addition & 2 deletions REQUIRE
Original file line number Diff line number Diff line change
@@ -1,2 +1 @@
julia 0.6
Compat
julia 0.7
12 changes: 2 additions & 10 deletions src/Expokit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,18 +11,10 @@ R.B. Sidje, ACM Trans. Math. Softw., 24(1):130-156, 1998 (or its
"""
module Expokit

using Compat, Compat.LinearAlgebra
import Compat:view, String
const LinearAlgebra = Compat.LinearAlgebra

if VERSION < v"0.7-"
nothing
else
using LinearAlgebra, SparseArrays
end
using LinearAlgebra

const axpy! = LinearAlgebra.axpy!
const expm! = LinearAlgebra.expm!
const exp! = LinearAlgebra.exp!

include("expmv.jl")
include("padm.jl")
Expand Down
4 changes: 2 additions & 2 deletions src/chbv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ chbv!(A, vec::Vector{T}) where {T} = chbv!(vec, A, copy(vec))

function chbv!(w::Vector{T}, A, vec::Vector{T}) where {T<:Real}
p = min(length(θ), length(α))
scale!(copy!(w, vec), α0)
rmul!(copyto!(w, vec), α0)
@inbounds for i = 1:p
w .+= real((A - θ[i]*I) \ (α[i] * vec))
end
Expand All @@ -86,7 +86,7 @@ end

function chbv!(w::Vector{T}, A, vec::Vector{T}) where {T<:Complex}
p = min(length(θ), length(α))
scale!(copy!(w, vec), α0)
rmul!(copyto!(w, vec), α0)
t = [θ; θconj]
a = 0.5 * [α; αconj]
@inbounds for i = 1:2*p
Expand Down
34 changes: 17 additions & 17 deletions src/expmv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,11 @@ export expmv, expmv!

function default_anorm(A)
try
Compat.opnorm(A, Inf)
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.")
@warn "opnorm($(typeof(A)), Inf) is not defined, fall back to using `anorm = 1.0`.
To suppress this warning, please specify the `anorm` keyword manually."
1.0
else
throw(err)
Expand All @@ -28,7 +28,7 @@ function expmv( t::Number,
A, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30, size(A, 1)),
norm=Base.norm, anorm=default_anorm(A)) where {T}
norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T}

result = convert(Vector{promote_type(eltype(A), T, typeof(t))}, copy(vec))
expmv!(t, A, result; tol=tol, m=m, norm=norm, anorm=anorm)
Expand All @@ -39,10 +39,10 @@ expmv!( t::Number,
A, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30,size(A,1)),
norm=Base.norm, anorm=default_anorm(A)) where {T} = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm, anorm=anorm)
norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T} = expmv!(vec, t, A, vec; tol=tol, m=m, norm=norm, anorm=anorm)

function expmv!(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)) where {T}
tol::Real=1e-7, m::Int=min(30,size(A,1)), norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T}

if size(vec,1) != size(A,2)
error("dimension mismatch")
Expand All @@ -62,10 +62,10 @@ function expmv!(w::Vector{T}, t::Number, A, vec::Vector{T};
r = 1/m
fact = (((m+1)/exp(1.0))^(m+1))*sqrt(2.0*pi*(m+1))
tau = (1.0/anorm)*((fact*tol)/(4.0*beta*anorm))^r
tau = signif(tau, 2)
tau = round(tau, sigdigits=2)

# storage for Krylov subspace vectors
vm = Array{typeof(w)}(m+1)
vm = Array{typeof(w)}(undef,m+1)
for i=1:m+1
vm[i]=similar(w)
end
Expand All @@ -75,19 +75,19 @@ function expmv!(w::Vector{T}, t::Number, A, vec::Vector{T};
tsgn = sign(t)
tk = zero(tf)

copy!(w, vec)
copyto!(w, vec)
p = similar(w)
mx = m
while tk < tf
tau = min(tf-tk, tau)

# Arnoldi procedure
# vm[1] = v/beta
scale!(copy!(vm[1],w),1/beta)
rmul!(copyto!(vm[1],w),1/beta)
mx = m
for j=1:m
# p[:] = A*vm[j]
Base.A_mul_B!(p, A, vm[j])
mul!(p, A, vm[j])

for i=1:j
hm[i,j] = dot(vm[i], p)
Expand All @@ -102,7 +102,7 @@ function expmv!(w::Vector{T}, t::Number, A, vec::Vector{T};

# F = expm(tsgn*tau*hm[1:j,1:j])
# F = expm!(scale(tsgn*tau,view(hm,1:j,1:j)))
F = expm!(tsgn*tau*view(hm,1:j,1:j))
F = exp!(tsgn*tau*view(hm,1:j,1:j))

fill!(w, zero(T))
for k=1:j
Expand All @@ -115,18 +115,18 @@ function expmv!(w::Vector{T}, t::Number, A, vec::Vector{T};
hm[j+1,j] = s

# vm[j+1] = p/hm[j+1,j]
scale!(copy!(vm[j+1],p),1/hm[j+1,j])
rmul!(copyto!(vm[j+1],p),1/hm[j+1,j])
end
hm[m+2,m+1] = one(T)
(mx != m) || (avnorm = norm(Base.A_mul_B!(p,A,vm[m+1])))
(mx != m) || (avnorm = norm(mul!(p,A,vm[m+1])))

# propagate using adaptive step size
iter = 1
while (iter < maxiter) && (mx == m)

# F = expm(tsgn*tau*hm)
# F = expm!(scale(tsgn*tau,hm))
F = expm!(tsgn*tau*hm)
F = exp!(tsgn*tau*hm)

# local error estimation
err1 = abs( beta*F[m+1,1] )
Expand Down Expand Up @@ -154,7 +154,7 @@ function expmv!(w::Vector{T}, t::Number, A, vec::Vector{T};
break
end
tau = gamma * tau * (tau*tol/err_loc)^r # estimate new time-step
tau = signif(tau, 2) # round to 2 significant digits
tau = round(tau, sigdigits=2) # round to 2 significant digits
# to prevent numerical noise
iter = iter + 1
end
Expand All @@ -167,7 +167,7 @@ function expmv!(w::Vector{T}, t::Number, A, vec::Vector{T};
tk = tk + tau

tau = gamma * tau * (tau*tol/err_loc)^r # estimate new time-step
tau = signif(tau, 2) # round to 2 significant digits
tau = round(tau, sigdigits=2) # round to 2 significant digits
# to prevent numerical noise
err_loc = max(err_loc,rndoff)

Expand Down
10 changes: 5 additions & 5 deletions src/padm.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ function padm(A; p::Int64=6)
end

# scaling
normA = norm(A, Inf)
normA = opnorm(A, Inf)
s = 0
if normA > 0.5
s = max(0, round(Int64, log(normA)/log(2), RoundToZero) + 2)
Expand All @@ -65,8 +65,8 @@ function padm(A; p::Int64=6)

# Horner evaluation of the irreducible fraction
A2 = A * A
Q = c[p+1]*eye(A)
P = c[p]*eye(A)
Q = c[p+1]*Matrix{eltype(A)}(I, size(A))
P = c[p]*Matrix{eltype(A)}(I, size(A))
odd = 1
@inbounds begin
for k = p-1:-1:1
Expand All @@ -82,11 +82,11 @@ function padm(A; p::Int64=6)
if odd == 1
Q = Q * A
Q = Q - P
E = -(I + 2 * \(Q, full(P)))
E = -(I + 2 * \(Q, Matrix(P)))
else
P = P * A
Q = Q - P
E = I + 2 * \(Q, full(P))
E = I + 2 * \(Q, Matrix(P))
end

# squaring
Expand Down
25 changes: 12 additions & 13 deletions src/phimv.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,12 +31,11 @@ Calculate the solution of a nonhomogeneous linear ODE problem with constant inpu
EXPOKIT: Software Package for Computing Matrix Exponentials.
ACM - Transactions On Mathematical Software, 24(1):130-156, 1998
"""

function phimv( t::Number,
A, u::Vector{T}, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30, size(A, 1)),
norm=Base.norm, anorm=default_anorm(A)) where {T}
norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T}

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)
Expand All @@ -47,10 +46,10 @@ phimv!( t::Number,
A, u::Vector{T}, vec::Vector{T};
tol::Real=1e-7,
m::Int=min(30, size(A, 1)),
norm=Base.norm, anorm=default_anorm(A)) where {T} = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm, anorm=anorm)
norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T} = phimv!(vec, t, A, u, vec; tol=tol, m=m, norm=norm, anorm=anorm)

function phimv!( 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=default_anorm(A)) where {T}
tol::Real=1e-7, m::Int=min(30, size(A, 1)), norm=LinearAlgebra.norm, anorm=default_anorm(A)) where {T}

if size(vec, 1) != size(A, 2)
error("dimension mismatch")
Expand All @@ -70,10 +69,10 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
r = 1/m
fact = (((m+1)/exp(1))^(m+1))*sqrt(2*pi*(m+1))
tau = (1.0/anorm)*((fact*tol)/(4.0*beta*anorm))^r
tau = signif(tau, 2)
tau = round(tau, sigdigits=2)

# storage for Krylov subspace vectors
vm = Array{typeof(w)}(m+1)
vm = Array{typeof(w)}(undef,m+1)
for i = 1:m+1
vm[i] = similar(w)
end
Expand All @@ -83,7 +82,7 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
tsgn = sign(t)
tk = zero(tf)

copy!(w, vec)
copyto!(w, vec)
p = similar(w)
k1 = 3
mb = m
Expand All @@ -92,10 +91,10 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
tau = min(tf-tk, tau)

# Arnoldi procedure
scale!(copy!(vm[1], A*w+u), 1/beta)
rmul!(copyto!(vm[1], A*w+u), 1/beta)

for j = 1:m
Base.A_mul_B!(p, A, vm[j])
mul!(p, A, vm[j])

for i = 1:j
hm[i, j] = dot(vm[i], p)
Expand All @@ -110,7 +109,7 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
break
end
hm[j+1,j] = s
scale!(copy!(vm[j+1], p), 1/hm[j+1,j])
rmul!(copyto!(vm[j+1], p), 1/hm[j+1,j])
end

hm[1, mb+1] = one(T)
Expand All @@ -125,7 +124,7 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
iter = 0
while iter <= maxiter
mx = mb + max(1,k1)
F = expm!(tsgn*tau*view(hm, 1:mx, 1:mx))
F = exp!(tsgn*tau*view(hm, 1:mx, 1:mx))
if k1 == 0
err_loc = btol
break
Expand All @@ -151,7 +150,7 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
break
else
tau = gamma * tau * (tau*tol/err_loc)^r # estimate new time-step
tau = signif(tau, 2) # round to 2 significant digits
tau = round(tau, sigdigits=2) # round to 2 significant digits
# to prevent numerical noise
if iter == maxiter
error("Number of iterations exceeded $(maxiter). Requested tolerance might be too high.")
Expand All @@ -169,7 +168,7 @@ function phimv!( w::Vector{T}, t::Number, A, u::Vector{T}, vec::Vector{T};
tk = tk + tau

tau = gamma * tau * (tau*tol/err_loc)^r # estimate new time-step
tau = signif(tau, 2) # round to 2 significant digits
tau = round(tau, sigdigits=2) # round to 2 significant digits
# to prevent numerical noise

err_loc = max(err_loc, rndoff)
Expand Down
Loading

0 comments on commit 4f91a8a

Please sign in to comment.