From 60fe1e503ebef66cb629133ba04bb8b306928e31 Mon Sep 17 00:00:00 2001 From: jecs Date: Wed, 23 Oct 2024 13:18:09 -0400 Subject: [PATCH 1/4] Added StaticArrays extension for expv method * Added ext/ExponentialUtilitiesStaticArraysExt.jl, which computes the exponential matrix-vector product efficiently when using StaticArrays. * Updated authors list, version #, and uuid in Project.toml. * Changed == nothing and === nothing to isnothing. --- .gitignore | 3 +- Project.toml | 12 +- ext/ExponentialUtilitiesStaticArraysExt.jl | 143 +++++++++++++++++++++ src/krylov_phiv.jl | 8 +- src/phi.jl | 8 +- test/basictests.jl | 12 ++ 6 files changed, 174 insertions(+), 12 deletions(-) create mode 100644 ext/ExponentialUtilitiesStaticArraysExt.jl diff --git a/.gitignore b/.gitignore index c49ecce..75341ba 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ *.jl.mem deps/deps.jl Manifest.toml -.vscode \ No newline at end of file +.vscode +.DS_Store diff --git a/Project.toml b/Project.toml index 99bdc86..f689521 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ExponentialUtilities" -uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" -authors = ["Chris Rackauckas "] -version = "1.26.1" +uuid = "3927d2ad-0af0-4d29-8d72-bf799b39c66b" +authors = ["Chris Rackauckas ", "José E. Cruz Serrallés "] +version = "1.27.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -14,6 +14,12 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" libblastrampoline_jll = "8e850b90-86db-534c-a0d3-1478176c7d93" +[weakdeps] +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[extensions] +ExponentialUtilitiesStaticArraysExt = "StaticArrays" + [compat] Adapt = "3.4.0, 4" Aqua = "0.8" diff --git a/ext/ExponentialUtilitiesStaticArraysExt.jl b/ext/ExponentialUtilitiesStaticArraysExt.jl new file mode 100644 index 0000000..b82e261 --- /dev/null +++ b/ext/ExponentialUtilitiesStaticArraysExt.jl @@ -0,0 +1,143 @@ +#= +Based on +Al-Mohy, A. H., & Higham, N. J. (2011). "Computing the action of the matrix exponential, with +an application to exponential integrators". SIAM journal on scientific computing, 33(2), 488-511. +=# +module ExponentialUtilitiesStaticArraysExt + +export default_tolerance,theta,THETA32,THETA64 + +using StaticArrays +import Base: @propagate_inbounds +import LinearAlgebra: tr,I,opnorm,norm +import ExponentialUtilities + +# Look-Up Table Generation +default_tolerance(::Type{T}) where {T <: AbstractFloat} = eps(T)/2 +@inline function trexp(M::Integer,x::T) where {T} + y = T <: BigInt ? one(BigFloat) : T <: Integer ? one(Float64) : one(T) + for m ∈ M:-1:1 + y = 1+x/m*y + end + return y +end +h(M::Integer,x::Number) = log(exp(-x)*trexp(M,x)) +h̃(M::Integer,x::Number) = ifelse(isodd(M),-1,1)*h(M,-x) +function θf((M,ϵ)::Tuple{<:Integer,<:Number},x::Number) + return h̃(M+1,x)/x-ϵ +end +θf(M::Integer,ϵ::Number) = Base.Fix1(θf,(M,ϵ)) +function θfp(M::Integer,x::Number) + Tk = trexp(M+1,-x) + Tkm1 = trexp(M,-x) + return ifelse(isodd(M),-1,1)/x^2*(log(Tk)+x*Tkm1/Tk) +end +θfp(M::Integer) = Base.Fix1(θfp,M) + +function newton_find_zero(f::Function,dfdx::Function,x0::Real;xrtol::Real=eps(typeof(x0))/2,maxiter::Integer=100) + 0 ≤ xrtol ≤ 1 || throw(DomainError(xrtol,"relative tolerance in x must be in [0,1]")) + maxiter > 0 || throw(DomainError(maxiter,"maxiter should be a positive integer")) + x, xp = x0, typemax(x0) + for _ ∈ 1:maxiter + xp = x + x -= f(x)/dfdx(x) + if abs(x-xp) ≤ xrtol*max(x,xp) || !isfinite(x) + break + end + end + return x +end +function calc_thetas(m_max::Integer,::Type{T};tol::T=default_tolerance(T)) where {T <: AbstractFloat} + m_max > 0 || throw(DomainError(m_max,"argument m_max must be positive")) + ϵ = BigFloat(tol) + θ = Vector{T}(undef,m_max+1) + @inbounds θ[1] = eps(T) + @inbounds for m=1:m_max + θ[m+1] = newton_find_zero(θf(m,ϵ),θfp(m),big(θ[m]),xrtol=ϵ) + end + return θ +end + +const P_MAX = 8 +const M_MAX = 55 +const THETA32 = Tuple(calc_thetas(M_MAX,Float32)) +const THETA64 = Tuple(calc_thetas(M_MAX,Float64)) + +@propagate_inbounds theta(::Type{Float64},m::Integer) = THETA64[m] +@propagate_inbounds theta(::Type{Float32},m::Integer) = THETA32[m] +@propagate_inbounds theta(::Type{Complex{T}},m::Integer) where {T} = theta(T,m) +@propagate_inbounds theta(::Type{T},::Integer) where {T} = throw(DomainError(T,"type must be either Float32 or Float64")) +@propagate_inbounds theta(x::Number,m::Integer) = theta(typeof(x),m) + +# runtime parameter search +@propagate_inbounds @inline function calculate_s(α::T,m::I)::I where {T <: Number,I <: Integer} + return ceil(I,α/theta(T,m)) +end +@propagate_inbounds @inline function parameter_search(nA::Number,m::I)::I where {I <: Integer} + return m*calculate_s(nA,m) +end +@propagate_inbounds @inline function parameters(A::SMatrix{N,N,T})::Tuple{Int,Int} where {N,T} + 1 ≤ N ≤ 50 || throw(DomainError(N,"leading dimension of A must be ≤ 50; larger matrices require Higham's 1-norm estimation algorithm")) + nA = opnorm(A,1) + iszero(nA) && return (0,1) + @inbounds if nA ≤ 4theta(T,M_MAX)*P_MAX*(P_MAX+3)/(M_MAX*1) + mo = argmin(Base.Fix1(parameter_search,nA),1:M_MAX) + s = calculate_s(nA,mo) + return (mo,s) + else + Aᵐ = A*A + pη = √(opnorm(Aᵐ,1)) + (Cmo::Int,mo::Int) = (typemax(Int),1) + for p ∈ 2:P_MAX + Aᵐ *= A + η = opnorm(Aᵐ,1)^inv(p+1) + α = max(pη,η) + pη = η + (Cmp::Int,mp::Int) = findmin(Base.Fix1(parameter_search,α),p*(p-1)-1:M_MAX) + (Cmo,mo) = min((Cmp,mp),(Cmo,mo)) + end + s = max(Cmo÷mo,1) + return (mo,s) + end +end + +# exponential matrix-vector product for SArray types +""" + expv(t::Number,A::SMatrix{N,N},v::SVector{N};kwarg...) → exp(t*A)*v + + Computes the matrix-vector product exp(t*A)*v without forming exp(t*A) explicitly. + This implementation is based on the algorithm presented in Al-Mohy & Higham (2011). + Presently, the relative tolerance is fixed to eps(T)/2 where T is the type of the + output. +""" +@propagate_inbounds function ExponentialUtilities.expv(t::Number,A::SMatrix{N,N,T},v::SVector{N}; kwarg...) where {N,T} + Ti = promote_type(StaticArrays.arithmetic_closure(T),eltype(v)) + N ≤ 4 && return exp(t*A)*v + Ai::SMatrix{N,N,Ti} = A + + μ = tr(Ai)/N + Ai -= μ*I + Ai *= t + mo, s = parameters(Ai) + F = v + Ai /= s + η = exp(μ*t/s) + ϵ = default_tolerance(T) + for _ ∈ 1:s + c₁ = norm(v,Inf) + for j ∈ 1:mo + v = (Ai*v)/j + F += v + c₂ = norm(v,Inf) + c₁+c₂ ≤ ϵ*norm(F,Inf) && break + c₁ = c₂ + end + F *= η + v = F + all(isfinite,v) || break + end + + return F +end + +end diff --git a/src/krylov_phiv.jl b/src/krylov_phiv.jl index eb34145..e2fad7b 100644 --- a/src/krylov_phiv.jl +++ b/src/krylov_phiv.jl @@ -77,7 +77,7 @@ function expv!(w::AbstractVector{Tw}, t::Real, Ks::KrylovSubspace{T, U}; cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, T, U} m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks) @assert length(w)==size(V, 1) "Dimension mismatch" - if cache == nothing + if isnothing(cache) cache = Matrix{U}(undef, m, m) elseif isa(cache, ExpvCache) cache = get_cache(cache, m) @@ -105,7 +105,7 @@ function expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspac cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, Tt, T, U} m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks) @assert length(w)==size(V, 1) "Dimension mismatch" - if cache === nothing + if isnothing(cache) cache = Matrix{U}(undef, m, m) elseif isa(cache, ExpvCache) cache = get_cache(cache, m) @@ -135,7 +135,7 @@ function ExponentialUtilities.expv!(w::GPUArraysCore.AbstractGPUVector{Tw}, expmethod = ExpMethodHigham2005Base()) where {Tw, T, U} m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks) @assert length(w)==size(V, 1) "Dimension mismatch" - if cache === nothing + if isnothing(cache) cache = Matrix{U}(undef, m, m) elseif isa(cache, ExpvCache) cache = get_cache(cache, m) @@ -259,7 +259,7 @@ function phiv!(w::AbstractMatrix, t::Number, Ks::KrylovSubspace{T, U}, k::Intege m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks) @assert size(w, 1)==size(V, 1) "Dimension mismatch" @assert size(w, 2)==k + 1 "Dimension mismatch" - if cache === nothing + if isnothing(cache) cache = PhivCache(w, m, k) elseif !isa(cache, PhivCache) throw(ArgumentError("Cache must be a PhivCache")) diff --git a/src/phi.jl b/src/phi.jl index 342eb4c..3c5ce8b 100644 --- a/src/phi.jl +++ b/src/phi.jl @@ -20,7 +20,7 @@ Software (TOMS), 24(1), 130-156. Theorem 1). function phi(z::T, k::Integer; cache = nothing, expmethod = ExpMethodHigham2005Base()) where {T <: Number} # Construct the matrix - if cache == nothing + if isnothing(cache) cache = fill(zero(T), k + 1, k + 1) else fill!(cache, zero(T)) @@ -66,7 +66,7 @@ function phiv_dense!(w::AbstractMatrix{T}, A::AbstractMatrix{T}, @assert size(w, 2)==k + 1 "Dimension mismatch" m = length(v) # Construct the extended matrix - if cache == nothing + if isnothing(cache) cache = fill(zero(T), m + k, m + k) else @assert size(cache)==(m + k, m + k) "Dimension mismatch" @@ -121,7 +121,7 @@ function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches = expmethod = ExpMethodHigham2005Base()) where {T <: Number} m = size(A, 1) @assert length(out) == k + 1&&all(P -> size(P) == (m, m), out) "Dimension mismatch" - if caches == nothing + if isnothing(caches) e = Vector{T}(undef, m) W = Matrix{T}(undef, m, k + 1) C = Matrix{T}(undef, m + k, m + k) @@ -143,7 +143,7 @@ function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches = end function phi!(out::Vector{Diagonal{T, V}}, A::Diagonal{T, V}, k::Integer; caches = nothing) where {T <: Number, V <: AbstractVector{T}} - for i in 1:size(A, 1) + for i in axes(A,1) phiz = phi(A[i, i], k; cache = caches) for j in 1:(k + 1) out[j][i, i] = phiz[j] diff --git a/test/basictests.jl b/test/basictests.jl index 25062db..e1f0d23 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -219,6 +219,18 @@ end end end +hasStaticArraysExt = !isnothing(Base.get_extension(ExponentialUtilities,:ExponentialUtilitiesStaticArraysExt)) +if hasStaticArraysExt + @testset "Static Arrays" begin + Random.seed!(0) + for N in (3,4,6,8),t in (0.1,1.0,10.0) + A = I+randn(SMatrix{N,N,Float64})/3 + b = randn(SVector{N,Float64}) + @test expv(t,A,b) ≈ exp(t*A)*b + end + end +end + @testset "Arnoldi & Krylov" begin Random.seed!(0) n = 20 From 1ca02b687e87b40a0002629ba3e8047f28d0b73b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 24 Oct 2024 04:40:41 +0000 Subject: [PATCH 2/4] Update test/basictests.jl --- test/basictests.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/test/basictests.jl b/test/basictests.jl index e1f0d23..6cb644e 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -219,15 +219,12 @@ end end end -hasStaticArraysExt = !isnothing(Base.get_extension(ExponentialUtilities,:ExponentialUtilitiesStaticArraysExt)) -if hasStaticArraysExt - @testset "Static Arrays" begin - Random.seed!(0) - for N in (3,4,6,8),t in (0.1,1.0,10.0) - A = I+randn(SMatrix{N,N,Float64})/3 - b = randn(SVector{N,Float64}) - @test expv(t,A,b) ≈ exp(t*A)*b - end +@testset "Static Arrays" begin + Random.seed!(0) + for N in (3,4,6,8),t in (0.1,1.0,10.0) + A = I+randn(SMatrix{N,N,Float64})/3 + b = randn(SVector{N,Float64}) + @test expv(t,A,b) ≈ exp(t*A)*b end end From c3a79595ea322d4e5a4cd4e3619f24a8c2398e0b Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 24 Oct 2024 04:40:47 +0000 Subject: [PATCH 3/4] Update ext/ExponentialUtilitiesStaticArraysExt.jl --- ext/ExponentialUtilitiesStaticArraysExt.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/ext/ExponentialUtilitiesStaticArraysExt.jl b/ext/ExponentialUtilitiesStaticArraysExt.jl index b82e261..2a25edb 100644 --- a/ext/ExponentialUtilitiesStaticArraysExt.jl +++ b/ext/ExponentialUtilitiesStaticArraysExt.jl @@ -1,8 +1,3 @@ -#= -Based on -Al-Mohy, A. H., & Higham, N. J. (2011). "Computing the action of the matrix exponential, with -an application to exponential integrators". SIAM journal on scientific computing, 33(2), 488-511. -=# module ExponentialUtilitiesStaticArraysExt export default_tolerance,theta,THETA32,THETA64 From f788dbf30691ab34057cf7f3d32fd5bde0747b44 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 24 Oct 2024 04:41:09 +0000 Subject: [PATCH 4/4] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f689521..2f07565 100644 --- a/Project.toml +++ b/Project.toml @@ -1,5 +1,5 @@ name = "ExponentialUtilities" -uuid = "3927d2ad-0af0-4d29-8d72-bf799b39c66b" +uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" authors = ["Chris Rackauckas ", "José E. Cruz Serrallés "] version = "1.27.0"