Add O(p m³) Al-Mohy–Liu scaling-and-recovering algorithm for the matrix φ-functions

Project.toml

@@ -1,13 +1,14 @@
 name = "ExponentialUtilities"
 uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "José E. Cruz Serrallés <Jose.CruzSerralles@nyulangone.org>"]
-version = "1.30.2"
+version = "1.31.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -25,10 +26,11 @@ Adapt = "3.4.0, 4"
 Aqua = "0.8"
 ArrayInterface = "7.1"
 DoubleFloats = "1.1.14"
-ForwardDiff = "0.10.13"
+ForwardDiff = "0.10.13, 1"
 GPUArraysCore = "0.1, 0.2"
 GenericSchur = "0.5.3"
 LinearAlgebra = "1.10"
+LinearSolve = "4"
 PrecompileTools = "1"
 Printf = "1.10"
 Random = "1"

docs/src/matrix_exponentials.md

@@ -3,6 +3,8 @@
 ```@docs
 exponential!
 phi
+phi!
+PhiPadeCache
 ```
 
 ## Methods

src/ExponentialUtilities.jl

@@ -29,6 +29,7 @@ include("exp_noalloc.jl")
 include("exp_generic.jl")
 include("exp_sparse.jl")
 include("phi.jl")
+include("phi_almohy.jl")
 include("arnoldi.jl")
 include("krylov_phiv.jl")
 include("krylov_phiv_adaptive.jl")
@@ -41,7 +42,7 @@ 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!
+    StegrCache, PhiPadeCache, get_subspace_cache, exponential!
 export ExpMethodHigham2005,
     ExpMethodHigham2005Base, ExpMethodGeneric, ExpMethodNative,
     ExpMethodDiagonalization

src/exp_generic.jl

@@ -222,7 +222,12 @@ function exponential!(
         # be all Infs when there are both Infs and zero since Inf*0===NaN
         return x * sum(x * exp(nx))
     end
-    s = iszero(nx) ? 0 : intlog2(nx)
+    # `nx > 1` (not `iszero(nx)`) so zero-valued Duals with nonzero partials
+    # take the s = 0 branch: ForwardDiff >= 1 defines `iszero(::Dual)` to also
+    # require zero partials, which sent Dual(0, 1) into `intlog2` where the
+    # partials were dropped (and 2^intlog2(0) = 2^64 overflows to 0). For real
+    # nx in (0, 1], intlog2(nx) was already 0, so this changes nothing else.
+    s = nx > 1 ? intlog2(nx) : 0
     (Vk === Val{13}() && x isa AbstractMatrix && ismutable(x)) &&
         return exp_generic_mutable(x, s, Val{13}())
     if s >= 1

src/phi.jl

@@ -103,9 +103,12 @@ The phi functions are defined as
 \\varphi_0(z) = \\exp(z),\\quad \\varphi_{k+1}(z) = \\frac{\\varphi_k(z) - \\varphi_k(0)}{z}
 ```
 
-Calls `phiv_dense` on each of the basis vectors to obtain the answer. If A
-is `Diagonal`, instead calls the scalar `phi` on each diagonal element and the
-return values are also `Diagonal`s
+For `Float64`/`ComplexF64` matrices this uses the scaling-and-recovering
+algorithm of Al-Mohy and Liu (arXiv:2506.01193), which computes all of
+`phi_0(A), ..., phi_k(A)` simultaneously in `O(k m^3)` operations. Other element
+types fall back to calling `phiv_dense` on each basis vector (`O(m (m+k)^3)`). If
+`A` is `Diagonal`, the scalar `phi` is instead applied to each diagonal element
+and the return values are also `Diagonal`s.
 """
 function phi(A::AbstractMatrix{T}, k; kwargs...) where {T <: Number}
     m = size(A, 1)
@@ -119,14 +122,56 @@ end
 """
     phi!(out,A,k[;caches]) -> out
 
-Non-allocating version of `phi` for non-diagonal matrix inputs.
+Non-allocating version of `phi` for non-diagonal matrix inputs, writing
+`phi_j(A)` into `out[j+1]`.
+
+For dense `Float64`/`ComplexF64` matrices, pass a reusable
+[`PhiPadeCache`](@ref) as `caches` to make repeated evaluations of the same
+size and order allocation-free (the linear solve runs through an embedded
+LinearSolve.jl cache):
+
+```julia
+cache = PhiPadeCache(A, k)
+phi!(out, A, k; caches = cache)
+```
+
+Numerical failure (a singular Padé denominator or non-finite result, only
+possible for pathological inputs such as matrices containing `NaN`/`Inf`) does
+not throw: the outputs are filled with `NaN` and, when a `PhiPadeCache` is
+used, `cache.info[]` is set to a nonzero return code (`0` on success). This
+lets adaptive integrators reject the step instead of aborting.
+
+For the legacy basis-vector algorithm, `caches` is instead the tuple
+`(Vector{T}(undef, m), Matrix{T}(undef, m, k+1), Matrix{T}(undef, m+k, m+k))`;
+supplying it forces that code path.
 """
 function phi!(
         out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches = nothing,
         expmethod = ExpMethodHigham2005Base()
     ) where {T <: Number}
     m = size(A, 1)
     @assert length(out) == k + 1&&all(P -> size(P) == (m, m), out) "Dimension mismatch"
+    # The scaling-and-recovering algorithm of Al-Mohy and Liu (arXiv:2506.01193)
+    # computes phi_0..phi_k simultaneously in O(k m^3), versus O(m (m+k)^3) for
+    # the basis-vector approach below. Its Pade tables are tuned for double
+    # precision, so it is used only for Float64/ComplexF64. The legacy path is
+    # kept for other element types and when a caller supplies the legacy
+    # `caches` tuple.
+    if k >= 1 && T <: Union{Float64, ComplexF64} && !(caches isa Tuple)
+        if A isa StridedMatrix && (isnothing(caches) || caches isa PhiPadeCache)
+            cache = caches isa PhiPadeCache ? caches : PhiPadeCache(A, k)
+            return _phi_almohy!(out, A, k, cache)
+        elseif isnothing(caches) && !ismutable(A)
+            # Container-preserving path for immutable dense matrices (e.g.
+            # `SMatrix`); mutable non-strided types (e.g. sparse) fall through to
+            # the legacy path below.
+            Rm = _phi_almohy_generic(A, k)
+            for j in 1:(k + 1)
+                copyto!(out[j], Rm[j])
+            end
+            return out
+        end
+    end
     if isnothing(caches)
         e = Vector{T}(undef, m)
         W = Matrix{T}(undef, m, k + 1)