diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 9c79359..959ad88 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -1,2 +1,3 @@ style = "sciml" -format_markdown = true \ No newline at end of file +format_markdown = true +format_docstrings = true \ No newline at end of file diff --git a/docs/pages.jl b/docs/pages.jl index 8e44b5e..98e3b6a 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -4,5 +4,5 @@ pages = [ "Home" => "index.md", "matrix_exponentials.md", "expv.md", - "arnoldi.md", + "arnoldi.md" ] diff --git a/src/ExponentialUtilities.jl b/src/ExponentialUtilities.jl index 95630f7..1dd4919 100644 --- a/src/ExponentialUtilities.jl +++ b/src/ExponentialUtilities.jl @@ -38,11 +38,11 @@ include("krylov_phiv_error_estimate.jl") include("precompile.jl") export phi, phi!, KrylovSubspace, arnoldi, arnoldi!, lanczos!, ExpvCache, PhivCache, - expv, expv!, phiv, phiv!, kiops, expv_timestep, expv_timestep!, phiv_timestep, - phiv_timestep!, - StegrCache, get_subspace_cache, exponential! + expv, expv!, phiv, phiv!, kiops, expv_timestep, expv_timestep!, phiv_timestep, + phiv_timestep!, + StegrCache, get_subspace_cache, exponential! export ExpMethodHigham2005, - ExpMethodHigham2005Base, ExpMethodGeneric, ExpMethodNative, - ExpMethodDiagonalization + ExpMethodHigham2005Base, ExpMethodGeneric, ExpMethodNative, + ExpMethodDiagonalization end diff --git a/src/StegrWork.jl b/src/StegrWork.jl index 4a06455..b69062a 100644 --- a/src/StegrWork.jl +++ b/src/StegrWork.jl @@ -73,7 +73,7 @@ Allocate work arrays for diagonalization of real-symmetric tridiagonal matrices of sizes up to `n`×`n`. """ function StegrWork(::Type{T}, n::Integer, - jobz::Char = 'V', range::Char = 'A') where {T} + jobz::Char = 'V', range::Char = 'A') where {T} n = convert(BlasInt, n) dv = Array{T}(undef, n) ev = Array{T}(undef, n) diff --git a/src/arnoldi.jl b/src/arnoldi.jl index 236c330..c85396f 100644 --- a/src/arnoldi.jl +++ b/src/arnoldi.jl @@ -35,7 +35,7 @@ mutable struct KrylovSubspace{T, U, B, VType <: AbstractMatrix{T}, end function KrylovSubspace{T, U}(n::Integer, maxiter::Integer = 30, - augmented::Integer = false) where {T, U} + augmented::Integer = false) where {T, U} V = Matrix{T}(undef, n + augmented, maxiter + 1) H = fill(zero(U), maxiter + 1, maxiter + !iszero(augmented)) return KrylovSubspace{T, U, real(T), Matrix{T}, Matrix{U}}(maxiter, maxiter, augmented, @@ -101,7 +101,7 @@ the dimension of `Ks` is smaller than `m`. Springer, Cham. """ function arnoldi(A, b; m = min(30, size(A, 1)), ishermitian = LinearAlgebra.ishermitian(A), - kwargs...) + kwargs...) TA, Tb = eltype(A), eltype(b) T = promote_type(TA, Tb) n = length(b) @@ -112,7 +112,8 @@ function arnoldi(A, b; m = min(30, size(A, 1)), ishermitian = LinearAlgebra.ishe # H is a small dense array in tridiagonal form H = zeros(U, (m + 1, m)) - Ks = KrylovSubspace{T, U, real(T), typeof(V), typeof(H)}(m, m, false, zero(real(T)), false, V, H) + Ks = KrylovSubspace{T, U, real(T), typeof(V), typeof(H)}( + m, m, false, zero(real(T)), false, V, H) arnoldi!(Ks, A, b; m = m, ishermitian = ishermitian, kwargs...) end @@ -211,8 +212,8 @@ end Take the `j` 'th step of the Lanczos iteration. """ function arnoldi_step!(j::Integer, iop::Integer, A::AT, - V::AbstractMatrix{T}, H::AbstractMatrix{U}, - n::Int = -1, p::Int = -1) where {AT, T, U} + V::AbstractMatrix{T}, H::AbstractMatrix{U}, + n::Int = -1, p::Int = -1) where {AT, T, U} x, y = @view(V[:, j]), @view(V[:, j + 1]) applyA!(y, A, x, V, j, n, p) @@ -235,12 +236,11 @@ end Non-allocating version of `arnoldi`. """ function arnoldi!(Ks::KrylovSubspace{T1, U}, A::AT, b; - tol::Real = 1e-7, m::Int = min(Ks.maxiter, size(A, 1)), - ishermitian::Bool = LinearAlgebra.ishermitian(A isa Tuple ? first(A) : A), - opnorm = nothing, iop::Int = 0, - init::Int = 0, t::Number = NaN, mu::Number = NaN, - l::Int = -1) where {T1 <: Number, U <: Number, AT} - + tol::Real = 1e-7, m::Int = min(Ks.maxiter, size(A, 1)), + ishermitian::Bool = LinearAlgebra.ishermitian(A isa Tuple ? first(A) : A), + opnorm = nothing, iop::Int = 0, + init::Int = 0, t::Number = NaN, mu::Number = NaN, + l::Int = -1) where {T1 <: Number, U <: Number, AT} Ks.wasbreakdown = false ishermitian && @@ -277,10 +277,10 @@ end Take the `j`'th step of the Lanczos iteration. """ function lanczos_step!(j::Integer, A, - V::AbstractMatrix{T}, - u::AbstractVector{U}, - v::AbstractVector{B}, - n::Int = -1, p::Int = -1) where {B, T, U} + V::AbstractMatrix{T}, + u::AbstractVector{U}, + v::AbstractVector{B}, + n::Int = -1, p::Int = -1) where {B, T, U} x, y = @view(V[:, j]), @view(V[:, j + 1]) applyA!(y, A, x, V, j, n, p) α = u[j] = coeff(U, dot(x, y)) @@ -316,11 +316,10 @@ A variation of `arnoldi!` that uses the Lanczos algorithm for Hermitian matrices. """ function lanczos!(Ks::KrylovSubspace{T1, U, B}, A::AT, b; - tol = 1e-7, m = min(Ks.maxiter, size(A, 1)), - opnorm = nothing, - init::Int = 0, t::Number = NaN, mu::Number = NaN, - l::Int = -1) where {T1 <: Number, U <: Number, B, AT} - + tol = 1e-7, m = min(Ks.maxiter, size(A, 1)), + opnorm = nothing, + init::Int = 0, t::Number = NaN, mu::Number = NaN, + l::Int = -1) where {T1 <: Number, U <: Number, B, AT} Ks.wasbreakdown = false m > Ks.maxiter ? resize!(Ks, m) : Ks.m = m # might change if happy-breakdown occurs diff --git a/src/exp_baseexp.jl b/src/exp_baseexp.jl index b862528..cac99d0 100644 --- a/src/exp_baseexp.jl +++ b/src/exp_baseexp.jl @@ -5,7 +5,7 @@ The same as `ExpMethodHigham2005` but follows `Base.exp` closer. """ struct ExpMethodHigham2005Base end function alloc_mem(A::StridedMatrix{T}, - method::ExpMethodHigham2005Base) where {T <: LinearAlgebra.BlasFloat} + method::ExpMethodHigham2005Base) where {T <: LinearAlgebra.BlasFloat} n = LinearAlgebra.checksquare(A) A2 = Matrix{T}(undef, n, n) P = Matrix{T}(undef, n, n) @@ -21,7 +21,7 @@ end ## Non-allocating version of `LinearAlgebra.exp!`. Modifies `A` to ## become (approximately) `exp(A)`. function exponential!(A::StridedMatrix{T}, method::ExpMethodHigham2005Base, - cache = alloc_mem(A, method)) where {T <: LinearAlgebra.BlasFloat} + cache = alloc_mem(A, method)) where {T <: LinearAlgebra.BlasFloat} X = A n = LinearAlgebra.checksquare(A) # if ishermitian(A) diff --git a/src/exp_generic.jl b/src/exp_generic.jl index b3a7112..6f5f591 100644 --- a/src/exp_generic.jl +++ b/src/exp_generic.jl @@ -3,7 +3,7 @@ intlog2(x::T) where {T <: Integer} = T(8 * sizeof(T) - leading_zeros(x - one(T)) intlog2(x) = x > typemax(UInt) ? ceil(Int, log2(x)) : intlog2(ceil(UInt, x)) % Int function naivemul!(C::StridedMatrix{T}, A::StridedMatrix{T}, - B::StridedMatrix{T}) where {T <: LinearAlgebra.BlasFloat} + B::StridedMatrix{T}) where {T <: LinearAlgebra.BlasFloat} mul!(C, A, B) end function naivemul!(C, A, B) @@ -37,7 +37,7 @@ _const(A) = A _const(A::Array) = Base.Experimental.Const(A) # Separated to make it easier to test. @generated function naivemul!(C::AbstractMatrix{T}, A, B, Maxis, Naxis, ::Val{MU}, - ::Val{NU}) where {T, MU, NU} + ::Val{NU}) where {T, MU, NU} nrem_body = quote m = first(Maxis) - 1 while m < M - $(MU - 1) @@ -124,7 +124,7 @@ struct ExpMethodGeneric{T} end ExpMethodGeneric() = ExpMethodGeneric{Val(13)}(); function exponential!(x, method::ExpMethodGeneric{Vk}, - cache = alloc_mem(x, method)) where {Vk} + cache = alloc_mem(x, method)) where {Vk} nx = opnorm(x, 1) if isnan(nx) || nx > 4611686018427387904 # i.e. 2^62 since it would cause overflow in 2^s # This should (hopefully) ensure that the result is Inf or NaN depending on @@ -197,7 +197,7 @@ function exp_generic_core!(y1, y2, y3, x, ::Val{13}) end function exp_pade_p!(y1::AbstractMatrix{T}, y2::AbstractMatrix{T}, x::AbstractMatrix, - ::Val{13}, ::Val{13}) where {T} + ::Val{13}, ::Val{13}) where {T} N = size(x, 1) # check square is in `exp_generic` y1 .= x .* 1.5440497506703088e-17 for c in (2.8101705462199623e-15, 2.529153491597966e-13, 1.48377004840414e-11, diff --git a/src/kiops.jl b/src/kiops.jl index 5fca2c5..877e964 100644 --- a/src/kiops.jl +++ b/src/kiops.jl @@ -48,8 +48,8 @@ References: - Niesen, J. and Wright, W.M., 2012. Algorithm 919: A Krylov subspace algorithm for evaluating the ``φ``-functions appearing in exponential integrators. ACM Transactions on Mathematical Software (TOMS), 38(3), p.22 """ function kiops(tau_out, A, u; mmin::Int = 10, mmax::Int = 128, m::Int = min(mmin, mmax), - tol::Real = 1e-7, opnorm = LinearAlgebra.opnorm(A, Inf), iop::Int = 2, - ishermitian::Bool = LinearAlgebra.ishermitian(A), task1::Bool = false) + tol::Real = 1e-7, opnorm = LinearAlgebra.opnorm(A, Inf), iop::Int = 2, + ishermitian::Bool = LinearAlgebra.ishermitian(A), task1::Bool = false) n, ppo = size(u, 1), size(u, 2) p = ppo - 1 @@ -259,8 +259,8 @@ function kiops(tau_out, A, u; mmin::Int = 10, mmax::Int = 128, m::Int = min(mmin end Base.@propagate_inbounds function kiops_update_solution!(tau_now, tau, tau_out, w, l, V, F, - H, beta, j, n, step, numSteps, - reject, ireject) + H, beta, j, n, step, numSteps, + reject, ireject) # Yep, got the required tolerance update reject = reject + ireject step = step + 1 diff --git a/src/krylov_phiv_adaptive.jl b/src/krylov_phiv_adaptive.jl index c13917f..6b2ae6f 100644 --- a/src/krylov_phiv_adaptive.jl +++ b/src/krylov_phiv_adaptive.jl @@ -44,12 +44,12 @@ end Non-allocating version of `expv_timestep`. """ function expv_timestep!(u::AbstractVector{T}, t::tType, A, b::AbstractVector{T}; - kwargs...) where {T <: Number, tType <: Real} + kwargs...) where {T <: Number, tType <: Real} expv_timestep!(reshape(u, length(u), 1), [t], A, b; kwargs...) return u end function expv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, b::AbstractVector{T}; - kwargs...) where {T <: Number, tType <: Real} + kwargs...) where {T <: Number, tType <: Real} B = reshape(b, length(b), 1) phiv_timestep!(U, ts, A, B; kwargs...) end @@ -95,19 +95,19 @@ end Non-allocating version of `phiv_timestep`. """ function phiv_timestep!(u::AbstractVector{T}, t::tType, A, B::AbstractMatrix{T}; - kwargs...) where {T <: Number, tType <: Real} + kwargs...) where {T <: Number, tType <: Real} phiv_timestep!(reshape(u, length(u), 1), [t], A, B; kwargs...) return u end function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractMatrix{T}; - tau::Real = 0.0, - m::Int = min(10, size(A, 1)), tol::Real = 1e-7, - opnorm = LinearAlgebra.opnorm(A, Inf), iop::Int = 0, - correct::Bool = false, caches = nothing, adaptive = false, - delta::Real = 1.2, - ishermitian::Bool = LinearAlgebra.ishermitian(A), - gamma::Real = 0.8, NA::Int = 0, - verbose = false) where {T <: Number, tType <: Real} + tau::Real = 0.0, + m::Int = min(10, size(A, 1)), tol::Real = 1e-7, + opnorm = LinearAlgebra.opnorm(A, Inf), iop::Int = 0, + correct::Bool = false, caches = nothing, adaptive = false, + delta::Real = 1.2, + ishermitian::Bool = LinearAlgebra.ishermitian(A), + gamma::Real = 0.8, NA::Int = 0, + verbose = false) where {T <: Number, tType <: Real} # Choose initial timestep opnorm = opnorm isa Number ? opnorm : opnorm(A, Inf) # backward compatibility abstol = tol * opnorm @@ -246,7 +246,7 @@ function phiv_timestep!(U::AbstractMatrix{T}, ts::Vector{tType}, A, B::AbstractM end # Helper functions for phiv_timestep! function _phiv_timestep_adapt(m, tau, epsilon, m_old, tau_old, epsilon_old, q, kappa, - gamma, omega, maxtau, n, p, NA, iop, Hnorm, verbose) + gamma, omega, maxtau, n, p, NA, iop, Hnorm, verbose) # Compute new m and tau (Algorithm 4) if tau_old > tau q = log(tau / tau_old) / log(epsilon / epsilon_old) - 1 diff --git a/src/krylov_phiv_error_estimate.jl b/src/krylov_phiv_error_estimate.jl index feda44b..eff5d8f 100644 --- a/src/krylov_phiv_error_estimate.jl +++ b/src/krylov_phiv_error_estimate.jl @@ -84,7 +84,7 @@ function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T}, @. V[:, 1] = b / Ks.beta ε = atol + rtol * Ks.beta - verbose && @printf("Initial norm: β₀ %e, stopping threshold: %e\n", Ks.beta,ε) + verbose && @printf("Initial norm: β₀ %e, stopping threshold: %e\n", Ks.beta, ε) α = @diagview(H) β = @diagview(H, -1) @@ -99,7 +99,7 @@ function expv!(w::AbstractVector{T}, t::Number, A, b::AbstractVector{T}, # Saad, Y. (1992). Analysis of some Krylov subspace # approximations. SIAM Journal on Numerical Analysis. σ = β[j] * Ks.beta * abs(cache.v[j]) - verbose && @printf("iter %d, α[%d] %e, β[%d] %e, σ %e\n", j, j, α[j], j, β[j],σ) + verbose && @printf("iter %d, α[%d] %e, β[%d] %e, σ %e\n", j, j, α[j], j, β[j], σ) if σ < ε Ks.m = j break diff --git a/src/phi.jl b/src/phi.jl index 971834b..342eb4c 100644 --- a/src/phi.jl +++ b/src/phi.jl @@ -18,7 +18,7 @@ package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156. Theorem 1). """ function phi(z::T, k::Integer; cache = nothing, - expmethod = ExpMethodHigham2005Base()) where {T <: Number} + expmethod = ExpMethodHigham2005Base()) where {T <: Number} # Construct the matrix if cache == nothing cache = fill(zero(T), k + 1, k + 1) @@ -59,9 +59,9 @@ end Non-allocating version of `phiv_dense`. """ function phiv_dense!(w::AbstractMatrix{T}, A::AbstractMatrix{T}, - v::AbstractVector{T}, k::Integer; - cache = nothing, - expmethod = ExpMethodHigham2005Base()) where {T <: Number} + v::AbstractVector{T}, k::Integer; + cache = nothing, + expmethod = ExpMethodHigham2005Base()) where {T <: Number} @assert size(w, 1)==size(A, 1)==size(A, 2)==length(v) "Dimension mismatch" @assert size(w, 2)==k + 1 "Dimension mismatch" m = length(v) @@ -118,7 +118,7 @@ end Non-allocating version of `phi` for non-diagonal matrix inputs. """ function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches = nothing, - expmethod = ExpMethodHigham2005Base()) where {T <: Number} + 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 @@ -142,7 +142,7 @@ function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches = return out end function phi!(out::Vector{Diagonal{T, V}}, A::Diagonal{T, V}, k::Integer; - caches = nothing) where {T <: Number, V <: AbstractVector{T}} + caches = nothing) where {T <: Number, V <: AbstractVector{T}} for i in 1:size(A, 1) phiz = phi(A[i, i], k; cache = caches) for j in 1:(k + 1) diff --git a/src/precompile.jl b/src/precompile.jl index e83c0ee..ec31710 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -32,11 +32,11 @@ end precomp_ms = [ - #ExpMethodHigham2005(), - ExpMethodHigham2005Base(), - #ExpMethodGeneric(), - #ExpMethodNative(), - #ExpMethodDiagonalization(), + #ExpMethodHigham2005(), + ExpMethodHigham2005Base() + #ExpMethodGeneric(), + #ExpMethodNative(), + #ExpMethodDiagonalization(), ] Txs = [Float64] diff --git a/test/basictests.jl b/test/basictests.jl index 6f11828..25062db 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -1,7 +1,7 @@ using Test, LinearAlgebra, Random, SparseArrays, ExponentialUtilities using ExponentialUtilities: getH, getV, exponential!, ExpMethodNative, - ExpMethodDiagonalization, ExpMethodHigham2005, ExpMethodGeneric, - ExpMethodHigham2005Base, alloc_mem + ExpMethodDiagonalization, ExpMethodHigham2005, ExpMethodGeneric, + ExpMethodHigham2005Base, alloc_mem using ForwardDiff, StaticArrays @testset "exp!" begin @@ -15,7 +15,7 @@ using ForwardDiff, StaticArrays ExpMethodDiagonalization(), ExpMethodHigham2005(), ExpMethodHigham2005Base(), - ExpMethodGeneric(), + ExpMethodGeneric() ] for m in methodlist @@ -119,13 +119,13 @@ end if VERSION >= v"1.7" out = Pipe() res = redirect_stdout(out) do - expv_timestep(ts, reshape([1], (1, 1)), [1.]; verbose = true) + expv_timestep(ts, reshape([1], (1, 1)), [1.0]; verbose = true) end close(Base.pipe_writer(out)) @test occursin("Completed after 1 time step(s)", read(out, String)) else - res = expv_timestep(ts, reshape([1], (1, 1)), [1.]; verbose = true) + res = expv_timestep(ts, reshape([1], (1, 1)), [1.0]; verbose = true) end @test vec(res) ≈ exp.(ts) @@ -137,28 +137,28 @@ end # FiniteDiff.finite_difference_gradient(sum ∘ exp, A) @test ForwardDiff.gradient(sum ∘ exp_generic, A) ≈ [1057.6079767302406 1061.5362688188497 1065.4645657136762 1069.39286741472 1073.3211571002205 1077.249468413699 1081.1777677116343 1085.1060742188956 1089.0343494857439 1092.9626752178747 1096.890964903375 1100.8192569919843 1104.747561096137 1108.6758820220502 1112.6041596920072 1116.5324493775076 1120.460763094095 1124.3890623920304 1128.3173592868568 1132.245658584792 1136.173974704488 1140.1022643899885 1144.030546866163; - 1546.2686015296754 1551.8428067925008 1557.4170336833045 1562.9912269305864 1568.5654249840857 1574.1396398593458 1579.7138523314973 1585.2880503849965 1590.8622628571482 1596.4364681199734 1602.0106853983423 1607.584888258059 1613.159081505341 1618.7332915743837 1624.3075064496438 1629.881714115578 1635.455931393947 1641.0301270443374 1646.6043419195976 1652.1785327637708 1657.752740429705 1663.3269745298346 1668.9011653740079; - 2034.9292407477624 2042.1493471692606 2049.3694608000847 2056.5895816402353 2063.8097024803856 2071.02981611121 2078.2499273389253 2085.4700457759673 2092.6901642130088 2099.910270634507 2107.130393877766 2114.350493089938 2121.570616333197 2128.79073236713 2136.010845997954 2143.230966838105 2150.4510684533857 2157.6711844873184 2164.8913077305774 2172.1114261676194 2179.331554217096 2186.5516558323766 2193.771774269418; - 2523.5898631440887 2532.4558707242595 2541.3219119479522 2550.187929140558 2559.0539511393813 2567.9199659288784 2576.78599033081 2585.6520195389594 2594.518051150217 2603.384075552149 2612.2500975509724 2621.116121952904 2629.9821487579447 2638.848151531898 2647.7141975618083 2656.5802195606316 2665.446239156346 2674.3122491396257 2683.178283153992 2692.044300346598 2700.9103199423125 2709.776363569114 2718.642368746176; - 3012.2504783310887 3022.7624255196715 3033.27436309582 3043.786281447098 3054.298197395268 3064.8101421807423 3075.3220869662164 3085.8340173330384 3096.3459476998605 3106.8578804697913 3117.3698108366134 3127.881736397218 3138.393659554714 3148.905589921536 3159.4175154821405 3169.9294602676146 3180.441381022002 3190.9533210012587 3201.4652561742982 3211.9771937504465 3222.489109698616 3233.0010544840907 3243.512987254021; - 3500.911105533632 3513.0689562839966 3525.226797421926 3537.384628947421 3549.542472488459 3561.7003256419325 3573.8581523612097 3586.015995902248 3598.1738226215257 3610.331658953238 3622.489526525363 3634.647350841532 3646.805189576353 3658.9630114894126 3671.120879061538 3683.2787177963587 3695.4365589342883 3707.5943808473485 3719.7522267914956 3731.9100511076645 3744.067906664246 3756.225767027045 3768.3836033587577; - 3989.5717399455016 4003.3754942576475 4017.1792317480326 4030.9829932695047 4044.786757194085 4058.590473056492 4072.3942201593118 4086.197972068349 4100.001714364951 4113.805480692641 4127.609222989243 4141.412986913824 4155.216714791774 4169.02046910392 4182.824213803631 4196.627960906451 4210.431710412379 4224.235455112091 4238.039211827345 4251.8429445115125 4265.646701226768 4279.450453135804 4293.254188223081; - 4478.232371954263 4493.6820346344075 4509.1316828959 4524.581340769827 4540.030991434429 4555.480649308356 4570.930309585392 4586.379962653102 4601.829615720812 4617.279259176087 4632.728936274884 4648.1785749239425 4663.628242410305 4679.077907493558 4694.527553351942 4709.977206419651 4725.426871502905 4740.876519764398 4756.326184847651 4771.775830706035 4787.225483773745 4802.675148856998 4818.124801924709; - 4966.892999156807 4983.988577414276 5001.084124431332 5018.179695479476 5035.275256915185 5052.370806335351 5069.466386995929 5086.56194362542 5103.657507464238 5120.753083318599 5137.848637544982 5154.944198980691 5172.039770028834 5189.135324255217 5206.230890497143 5223.3264591421785 5240.422022980996 5257.517596432248 5274.613138643088 5291.708714497448 5308.804278336266 5325.899844578193 5342.99540841701; - 5455.553635971784 5474.29509135684 5493.036573176091 5511.778040576691 5530.519532008377 5549.26098979654 5568.002459600249 5586.743934210173 5605.485399207664 5624.226888236241 5642.968336411971 5661.70980381257 5680.451290438038 5699.1927674510725 5717.934239657889 5736.675716670922 5755.417176862195 5774.158646665904 5792.900109260285 5811.641574257776 5830.383070495679 5849.124537896278 5867.866010103095; - 5944.214258368111 5964.601636539817 5984.989024323959 6005.376400092557 6025.763761442502 6046.151166048405 6066.53852739835 6086.9259175856005 6107.313293354198 6127.70068113834 6148.088042488285 6168.475437481753 6188.862810847242 6209.250198631384 6229.637569593764 6250.024959781014 6270.412342758938 6290.7997161244275 6311.187099102352 6331.574472467841 6351.961841027113 6372.349233617471 6392.736621401613; - 6432.874883167546 6454.908169707251 6476.941465859391 6498.974752399096 6521.008041341911 6543.0413302847255 6565.074602405779 6587.10789615481 6609.141180291407 6631.17447163733 6653.20776778947 6675.241054329175 6697.2743288533375 6719.307610586826 6741.340925963836 6763.374205294215 6785.407489430812 6807.440792792278 6829.474062510223 6851.50734664682 6873.540640395851 6895.573917323122 6917.607213475262; - 6921.535507966981 6945.214717293336 6968.893909797933 6992.573109511854 7016.252299613341 7039.931506536589 7063.610687025642 7087.289881933346 7110.969084050376 7134.648286167407 7158.327473865786 7182.006666370381 7205.685868487411 7229.365068201333 7253.0442607059285 7276.723453210525 7300.402640908903 7324.081835816607 7347.761042739855 7371.440235244451 7395.119427749047 7418.798627462968 7442.477822370673; - 7410.1961447819585 7435.521238445227 7460.846356139583 7486.171457012178 7511.496562690989 7536.821680385345 7562.1467620330695 7587.471877324317 7612.796975793803 7638.122062247745 7663.447172732774 7688.772283217803 7714.097386493507 7739.422499381645 7764.747607463565 7790.072689111291 7815.397801999428 7840.722898065806 7866.048015760161 7891.37312624519 7916.698215102242 7942.023323184162 7967.348433669192; - 7898.856784000046 7925.8277980468565 7952.798792868798 7979.769828543587 8006.740820962419 8033.711822993687 8060.682839443606 8087.653858296635 8114.624865134119 8141.595881584039 8168.566883615306 8195.537909677661 8222.50890930582 8249.479918546413 8276.450947011876 8303.42195384936 8330.392951074411 8357.363974733657 8384.33498157114 8411.305986005516 8438.277016874088 8465.248026114683 8492.219037758385; - 8387.517403993263 8416.134326408073 8444.751251225991 8473.368171237693 8501.98507683074 8530.601992036225 8559.218914451034 8587.835829656518 8616.452759280653 8645.069684098573 8673.686587288512 8702.303519315758 8730.920446536784 8759.53735693605 8788.154281753968 8816.771192153235 8845.38812177737 8874.005029773529 8902.621940172796 8931.238874603148 8959.855797017957 8988.47271702966 9017.089646653794; - 8876.178026389589 8906.440864381724 8936.70368314899 8966.966497110037 8997.229342311497 9027.492173094306 9057.755011086441 9088.017827450598 9118.28065342719 9148.54347700067 9178.806319799023 9209.069126550745 9239.331950124228 9269.594790519472 9299.857633317823 9330.12043045711 9360.383270852353 9390.646101635162 9420.908927611754 9451.171768006996 9481.43458917737 9511.697424766397 9541.96024593677; - 9364.838670413894 9396.747399952266 9428.656122281312 9460.56486143212 9492.473617404688 9524.382334927519 9556.291074078326 9588.199798810481 9620.108533155071 9652.017279515205 9683.926013859795 9715.834753010602 9747.743477742757 9779.652214490456 9811.56096085059 9843.469690388963 9875.378434345988 9907.28718551234 9939.195893422733 9971.10465179841 10003.013395755435 10034.922110875155 10066.830857235289; - 9853.499285600894 9887.053928313482 9920.60857583229 9954.163220947987 9987.717863660575 10021.272501566948 10054.827149085753 10088.381801410778 10121.93643691404 10155.491101254607 10189.045729548545 10222.600367454916 10256.155007764397 10289.70966008942 10323.264307608226 10356.818945514598 10390.373585824078 10423.928238149103 10457.482880861691 10491.03752117117 10524.592161480652 10558.146825821219 10591.701454115157; - 10342.15992722209 10377.36047109335 10412.561017367721 10447.761582866962 10482.962136350658 10518.162670609485 10553.363216883856 10588.563779979988 10623.764331060574 10658.964884544272 10694.165430818643 10729.365986705447 10764.566544995361 10799.767081657297 10834.967647156536 10870.168195834016 10905.368749317713 10940.569285979649 10975.769858688214 11010.970392947042 11046.170956043174 11081.371499914434 11116.572070219892; - 10830.820552021523 10867.667009067003 10904.513456500046 10941.359930367285 10978.206382606546 11015.052854070675 11051.899306309937 11088.745756146089 11125.592222804002 11162.438677446371 11199.285141701175 11236.131596343545 11272.978065404566 11309.824522450044 11346.670989107955 11383.517436541 11420.363895989587 11457.210350631956 11494.056812483652 11530.903281544674 11567.749738590152 11604.596198038738 11641.442674309084; - 11319.48117441785 11357.973542234437 11396.465905244804 11434.958285076933 11473.450640877976 11511.94301109767 11550.435376511146 11588.927746730842 11627.420116950536 11665.91248717023 11704.404838165055 11742.897220400293 11781.389576201334 11819.881958436572 11858.374326253159 11896.866686860418 11935.359049870785 11973.851427299805 12012.343804728827 12050.836150917434 12089.328530749563 12127.820886550606 12166.313261576517; - 11808.141794411067 11848.280085014305 11888.418346780238 11928.556642189691 11968.694903955624 12008.833187349534 12048.971465937228 12089.109730106267 12129.248003887744 12169.386292087873 12209.524549047586 12249.662842053933 12289.801132657169 12329.939380004449 12370.07766339836 12410.215941986053 12450.354210961312 12490.492489549004 12530.630772942915 12570.769034708848 12610.907310893432 12651.045589481126 12691.183868068818] + 1546.2686015296754 1551.8428067925008 1557.4170336833045 1562.9912269305864 1568.5654249840857 1574.1396398593458 1579.7138523314973 1585.2880503849965 1590.8622628571482 1596.4364681199734 1602.0106853983423 1607.584888258059 1613.159081505341 1618.7332915743837 1624.3075064496438 1629.881714115578 1635.455931393947 1641.0301270443374 1646.6043419195976 1652.1785327637708 1657.752740429705 1663.3269745298346 1668.9011653740079; + 2034.9292407477624 2042.1493471692606 2049.3694608000847 2056.5895816402353 2063.8097024803856 2071.02981611121 2078.2499273389253 2085.4700457759673 2092.6901642130088 2099.910270634507 2107.130393877766 2114.350493089938 2121.570616333197 2128.79073236713 2136.010845997954 2143.230966838105 2150.4510684533857 2157.6711844873184 2164.8913077305774 2172.1114261676194 2179.331554217096 2186.5516558323766 2193.771774269418; + 2523.5898631440887 2532.4558707242595 2541.3219119479522 2550.187929140558 2559.0539511393813 2567.9199659288784 2576.78599033081 2585.6520195389594 2594.518051150217 2603.384075552149 2612.2500975509724 2621.116121952904 2629.9821487579447 2638.848151531898 2647.7141975618083 2656.5802195606316 2665.446239156346 2674.3122491396257 2683.178283153992 2692.044300346598 2700.9103199423125 2709.776363569114 2718.642368746176; + 3012.2504783310887 3022.7624255196715 3033.27436309582 3043.786281447098 3054.298197395268 3064.8101421807423 3075.3220869662164 3085.8340173330384 3096.3459476998605 3106.8578804697913 3117.3698108366134 3127.881736397218 3138.393659554714 3148.905589921536 3159.4175154821405 3169.9294602676146 3180.441381022002 3190.9533210012587 3201.4652561742982 3211.9771937504465 3222.489109698616 3233.0010544840907 3243.512987254021; + 3500.911105533632 3513.0689562839966 3525.226797421926 3537.384628947421 3549.542472488459 3561.7003256419325 3573.8581523612097 3586.015995902248 3598.1738226215257 3610.331658953238 3622.489526525363 3634.647350841532 3646.805189576353 3658.9630114894126 3671.120879061538 3683.2787177963587 3695.4365589342883 3707.5943808473485 3719.7522267914956 3731.9100511076645 3744.067906664246 3756.225767027045 3768.3836033587577; + 3989.5717399455016 4003.3754942576475 4017.1792317480326 4030.9829932695047 4044.786757194085 4058.590473056492 4072.3942201593118 4086.197972068349 4100.001714364951 4113.805480692641 4127.609222989243 4141.412986913824 4155.216714791774 4169.02046910392 4182.824213803631 4196.627960906451 4210.431710412379 4224.235455112091 4238.039211827345 4251.8429445115125 4265.646701226768 4279.450453135804 4293.254188223081; + 4478.232371954263 4493.6820346344075 4509.1316828959 4524.581340769827 4540.030991434429 4555.480649308356 4570.930309585392 4586.379962653102 4601.829615720812 4617.279259176087 4632.728936274884 4648.1785749239425 4663.628242410305 4679.077907493558 4694.527553351942 4709.977206419651 4725.426871502905 4740.876519764398 4756.326184847651 4771.775830706035 4787.225483773745 4802.675148856998 4818.124801924709; + 4966.892999156807 4983.988577414276 5001.084124431332 5018.179695479476 5035.275256915185 5052.370806335351 5069.466386995929 5086.56194362542 5103.657507464238 5120.753083318599 5137.848637544982 5154.944198980691 5172.039770028834 5189.135324255217 5206.230890497143 5223.3264591421785 5240.422022980996 5257.517596432248 5274.613138643088 5291.708714497448 5308.804278336266 5325.899844578193 5342.99540841701; + 5455.553635971784 5474.29509135684 5493.036573176091 5511.778040576691 5530.519532008377 5549.26098979654 5568.002459600249 5586.743934210173 5605.485399207664 5624.226888236241 5642.968336411971 5661.70980381257 5680.451290438038 5699.1927674510725 5717.934239657889 5736.675716670922 5755.417176862195 5774.158646665904 5792.900109260285 5811.641574257776 5830.383070495679 5849.124537896278 5867.866010103095; + 5944.214258368111 5964.601636539817 5984.989024323959 6005.376400092557 6025.763761442502 6046.151166048405 6066.53852739835 6086.9259175856005 6107.313293354198 6127.70068113834 6148.088042488285 6168.475437481753 6188.862810847242 6209.250198631384 6229.637569593764 6250.024959781014 6270.412342758938 6290.7997161244275 6311.187099102352 6331.574472467841 6351.961841027113 6372.349233617471 6392.736621401613; + 6432.874883167546 6454.908169707251 6476.941465859391 6498.974752399096 6521.008041341911 6543.0413302847255 6565.074602405779 6587.10789615481 6609.141180291407 6631.17447163733 6653.20776778947 6675.241054329175 6697.2743288533375 6719.307610586826 6741.340925963836 6763.374205294215 6785.407489430812 6807.440792792278 6829.474062510223 6851.50734664682 6873.540640395851 6895.573917323122 6917.607213475262; + 6921.535507966981 6945.214717293336 6968.893909797933 6992.573109511854 7016.252299613341 7039.931506536589 7063.610687025642 7087.289881933346 7110.969084050376 7134.648286167407 7158.327473865786 7182.006666370381 7205.685868487411 7229.365068201333 7253.0442607059285 7276.723453210525 7300.402640908903 7324.081835816607 7347.761042739855 7371.440235244451 7395.119427749047 7418.798627462968 7442.477822370673; + 7410.1961447819585 7435.521238445227 7460.846356139583 7486.171457012178 7511.496562690989 7536.821680385345 7562.1467620330695 7587.471877324317 7612.796975793803 7638.122062247745 7663.447172732774 7688.772283217803 7714.097386493507 7739.422499381645 7764.747607463565 7790.072689111291 7815.397801999428 7840.722898065806 7866.048015760161 7891.37312624519 7916.698215102242 7942.023323184162 7967.348433669192; + 7898.856784000046 7925.8277980468565 7952.798792868798 7979.769828543587 8006.740820962419 8033.711822993687 8060.682839443606 8087.653858296635 8114.624865134119 8141.595881584039 8168.566883615306 8195.537909677661 8222.50890930582 8249.479918546413 8276.450947011876 8303.42195384936 8330.392951074411 8357.363974733657 8384.33498157114 8411.305986005516 8438.277016874088 8465.248026114683 8492.219037758385; + 8387.517403993263 8416.134326408073 8444.751251225991 8473.368171237693 8501.98507683074 8530.601992036225 8559.218914451034 8587.835829656518 8616.452759280653 8645.069684098573 8673.686587288512 8702.303519315758 8730.920446536784 8759.53735693605 8788.154281753968 8816.771192153235 8845.38812177737 8874.005029773529 8902.621940172796 8931.238874603148 8959.855797017957 8988.47271702966 9017.089646653794; + 8876.178026389589 8906.440864381724 8936.70368314899 8966.966497110037 8997.229342311497 9027.492173094306 9057.755011086441 9088.017827450598 9118.28065342719 9148.54347700067 9178.806319799023 9209.069126550745 9239.331950124228 9269.594790519472 9299.857633317823 9330.12043045711 9360.383270852353 9390.646101635162 9420.908927611754 9451.171768006996 9481.43458917737 9511.697424766397 9541.96024593677; + 9364.838670413894 9396.747399952266 9428.656122281312 9460.56486143212 9492.473617404688 9524.382334927519 9556.291074078326 9588.199798810481 9620.108533155071 9652.017279515205 9683.926013859795 9715.834753010602 9747.743477742757 9779.652214490456 9811.56096085059 9843.469690388963 9875.378434345988 9907.28718551234 9939.195893422733 9971.10465179841 10003.013395755435 10034.922110875155 10066.830857235289; + 9853.499285600894 9887.053928313482 9920.60857583229 9954.163220947987 9987.717863660575 10021.272501566948 10054.827149085753 10088.381801410778 10121.93643691404 10155.491101254607 10189.045729548545 10222.600367454916 10256.155007764397 10289.70966008942 10323.264307608226 10356.818945514598 10390.373585824078 10423.928238149103 10457.482880861691 10491.03752117117 10524.592161480652 10558.146825821219 10591.701454115157; + 10342.15992722209 10377.36047109335 10412.561017367721 10447.761582866962 10482.962136350658 10518.162670609485 10553.363216883856 10588.563779979988 10623.764331060574 10658.964884544272 10694.165430818643 10729.365986705447 10764.566544995361 10799.767081657297 10834.967647156536 10870.168195834016 10905.368749317713 10940.569285979649 10975.769858688214 11010.970392947042 11046.170956043174 11081.371499914434 11116.572070219892; + 10830.820552021523 10867.667009067003 10904.513456500046 10941.359930367285 10978.206382606546 11015.052854070675 11051.899306309937 11088.745756146089 11125.592222804002 11162.438677446371 11199.285141701175 11236.131596343545 11272.978065404566 11309.824522450044 11346.670989107955 11383.517436541 11420.363895989587 11457.210350631956 11494.056812483652 11530.903281544674 11567.749738590152 11604.596198038738 11641.442674309084; + 11319.48117441785 11357.973542234437 11396.465905244804 11434.958285076933 11473.450640877976 11511.94301109767 11550.435376511146 11588.927746730842 11627.420116950536 11665.91248717023 11704.404838165055 11742.897220400293 11781.389576201334 11819.881958436572 11858.374326253159 11896.866686860418 11935.359049870785 11973.851427299805 12012.343804728827 12050.836150917434 12089.328530749563 12127.820886550606 12166.313261576517; + 11808.141794411067 11848.280085014305 11888.418346780238 11928.556642189691 11968.694903955624 12008.833187349534 12048.971465937228 12089.109730106267 12129.248003887744 12169.386292087873 12209.524549047586 12249.662842053933 12289.801132657169 12329.939380004449 12370.07766339836 12410.215941986053 12450.354210961312 12490.492489549004 12530.630772942915 12570.769034708848 12610.907310893432 12651.045589481126 12691.183868068818] for n in 1:48 A = rand(n, n) @@ -287,7 +287,7 @@ end Hermitian(rand(ComplexF64, n, n)), Hermitian(rand(n, n)), rand(ComplexF64, n, n), - rand(n, n), + rand(n, n) ] for b in [rand(ComplexF64, n), rand(n)], t in [1e-2, 1e-2im, 1e-2 + 1e-2im] @test exp(t * A) * b ≈ expv(t, A, b; m = m) diff --git a/test/gpu/gputests.jl b/test/gpu/gputests.jl index 3150018..29c5d54 100644 --- a/test/gpu/gputests.jl +++ b/test/gpu/gputests.jl @@ -3,8 +3,8 @@ using SparseArrays using CUDA using CUDA.CUSPARSE using ExponentialUtilities: inplace_add!, - exponential!, ExpMethodHigham2005, expv, - expv_timestep + exponential!, ExpMethodHigham2005, expv, + expv_timestep using Random: Xoshiro