Add missing documentation for exponentiation methods

src/exp.jl

@@ -13,6 +13,15 @@ end
 @deprecate _exp! exponential!
 @deprecate exp_generic exponential!
 exponential!(A) = exponential!(A, ExpMethodHigham2005(A));
+
+"""
+    exponential!(A::GPUArraysCore.AbstractGPUArray)
+
+Computes the matrix exponential for GPU arrays. Automatically disables balancing
+for GPU arrays as it is not supported on GPU devices.
+
+See also: [`exponential!`](@ref)
+"""
 function exponential!(A::GPUArraysCore.AbstractGPUArray)
     exponential!(A, ExpMethodHigham2005(false))
 end;

src/exp_generated/exp_1.jl

@@ -1,5 +1,30 @@
 using LinearAlgebra
 
+"""
+    exp_gen!(cache, A, ::Val{s})
+
+Internal function implementing the matrix exponential using Padé approximation of order `s`.
+
+This is part of the Higham (2005) scaling and squaring algorithm implementation. Each
+`exp_gen!` function (for s=1 to 13) implements a different order Padé approximation for
+computing the matrix exponential. The function modifies `A` in-place to contain exp(A).
+
+# Arguments
+
+  - `cache`: Pre-allocated workspace arrays for intermediate computations
+  - `A`: Input matrix that will be overwritten with exp(A)
+  - `::Val{s}`: Order of the Padé approximation (s ∈ {1,2,...,13})
+
+# Notes
+
+  - This is an internal implementation detail and should not be called directly
+  - Higher orders provide better accuracy but require more computation
+  - The choice of order is made automatically by the ExpMethodHigham2005 algorithm
+  - Generated code from the algorithm in: Higham, N. J. (2005). "The scaling and squaring
+    method for the matrix exponential revisited." SIAM J. Matrix Anal. Appl. 26(4), 1179-1193.
+
+See also: [`exponential!`](@ref), [`ExpMethodHigham2005`](@ref)
+"""
 function exp_gen!(cache, A, ::Val{1})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=4

src/exp_generated/exp_10.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 10
 function exp_gen!(cache, A, ::Val{10})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_11.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 11
 function exp_gen!(cache, A, ::Val{11})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_12.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 12
 function exp_gen!(cache, A, ::Val{12})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_13.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 13
 function exp_gen!(cache, A, ::Val{13})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_2.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 2
 function exp_gen!(cache, A, ::Val{2})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=5

src/exp_generated/exp_3.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 3
 function exp_gen!(cache, A, ::Val{3})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_4.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 4
 function exp_gen!(cache, A, ::Val{4})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_5.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 5
 function exp_gen!(cache, A, ::Val{5})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_6.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 6
 function exp_gen!(cache, A, ::Val{6})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_7.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 7
 function exp_gen!(cache, A, ::Val{7})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_8.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 8
 function exp_gen!(cache, A, ::Val{8})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_generated/exp_9.jl

@@ -1,5 +1,7 @@
 using LinearAlgebra
 
+# See exp_gen!(cache, A, ::Val{1}) for documentation
+# This implements the Padé approximation of order 9
 function exp_gen!(cache, A, ::Val{9})
     T = promote_type(eltype(A), Float64) # Make it work for many 'bigger' types (matrices and scalars)
     # max_memslots=6

src/exp_noalloc.jl

@@ -19,13 +19,20 @@ function alloc_mem(A, ::ExpMethodHigham2005)
 end
 
 # Import the generated code
+# exp_gen! functions are internal implementation details of the ExpMethodHigham2005 algorithm.
+# They implement different Padé approximation orders (1-13) for computing matrix exponentials.
+# These functions are automatically generated and should not be called directly by users.
 for i in 1:13
     include("exp_generated/exp_$i.jl")
 end
 
+# Internal helper function for accessing cache memory slots (used by generated code)
 function getmem(cache, k) # Called from generated code
     return cache[k - 1]
 end
+
+# Internal helper function for in-place left division (used by generated code)
+# Computes C = A\B, storing result in C
 function ldiv_for_generated!(C, A, B) # C=A\B. Called from generated code
     F = lu!(A) # This allocation is unavoidable, due to the interface of LinearAlgebra
     ldiv!(F, B) # Result stored in B

src/krylov_phiv.jl

@@ -42,6 +42,7 @@ function expv(t, A, b; mode = :happy_breakdown, kwargs...)
         throw(ArgumentError("Unknown Krylov iteration termination mode, $(mode)"))
     end
 end
+# Internal helper function for happy-breakdown mode
 function _expv_hb(t::Tt, A, b;
         expmethod = ExpMethodHigham2005Base(),
         cache = nothing, kwargs_arnoldi...) where {Tt}
@@ -50,6 +51,8 @@ function _expv_hb(t::Tt, A, b;
     w = similar(b, promote_type(Tt, eltype(A), eltype(b)))
     expv!(w, t, Ks; cache = cache, expmethod = expmethod)
 end
+
+# Internal helper function for error-estimate mode
 function _expv_ee(t::Tt, A, b; m = min(30, size(A, 1)), tol = 1e-7, rtol = √(tol),
         ishermitian::Bool = LinearAlgebra.ishermitian(A),
         expmethod = ExpMethodHigham2005Base()) where {Tt}
@@ -62,6 +65,19 @@ function _expv_ee(t::Tt, A, b; m = min(30, size(A, 1)), tol = 1e-7, rtol = √(t
     expv!(w, t, A, b, Ks, get_subspace_cache(Ks); atol = tol, rtol = rtol,
         ishermitian = ishermitian, expmethod = expmethod)
 end
+"""
+    expv(t, Ks::KrylovSubspace; expmethod=ExpMethodHigham2005(), kwargs...) -> exp(tA)b
+
+Compute the matrix-exponential-vector product using a precomputed Krylov subspace `Ks`.
+
+# Arguments
+
+  - `t`: Time parameter for the exponential
+  - `Ks`: Precomputed Krylov subspace containing the approximation data
+  - `expmethod`: Method for computing the matrix exponential (default: ExpMethodHigham2005())
+
+See also: [`expv!`](@ref), [`arnoldi`](@ref)
+"""
 function expv(t::Tt, Ks::KrylovSubspace{T, U}; expmethod = ExpMethodHigham2005(),
         kwargs...) where {Tt, T, U}
     n = size(getV(Ks), 1)
@@ -98,9 +114,28 @@ function expv!(w::AbstractVector{Tw}, t::Real, Ks::KrylovSubspace{T, U};
     lmul!(beta, mul!(w, @view(V[:, 1:m]), expHe)) # exp(A) ≈ norm(b) * V * exp(H)e
 end
 
-# NOTE: Tw can be Float64, while t is ComplexF64 and T is Float32
-#       or Tw can be Float64, while t is ComplexF32 and T is Float64
-#       thus they can not share the same TypeVar.
+"""
+    expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspace{T, U}; 
+          cache=nothing, expmethod=ExpMethodHigham2005Base()) -> w
+
+Non-allocating version of `expv` for complex time parameter `t` and complex output vector `w`.
+
+# Arguments
+
+  - `w`: Pre-allocated output vector for the result
+  - `t`: Complex time parameter for the exponential
+  - `Ks`: Precomputed Krylov subspace containing the approximation data
+  - `cache`: Optional pre-allocated workspace matrix
+  - `expmethod`: Method for computing the matrix exponential
+
+# Notes
+
+  - Tw can be Float64, while t is ComplexF64 and T is Float32
+  - or Tw can be Float64, while t is ComplexF32 and T is Float64
+  - thus they can not share the same TypeVar.
+
+See also: [`expv`](@ref)
+"""
 function expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspace{T, U};
         cache = nothing, expmethod = ExpMethodHigham2005Base()) where {Tw, Tt, T, U}
     m, beta, V, H = Ks.m, Ks.beta, getV(Ks), getH(Ks)
@@ -129,6 +164,22 @@ function expv!(w::AbstractVector{Complex{Tw}}, t::Complex{Tt}, Ks::KrylovSubspac
     lmul!(beta, mul!(w, @view(V[:, 1:m]), compatible_multiplicative_operand(V, expHe))) # exp(A) ≈ norm(b) * V * exp(H)e
 end
 
+"""
+    expv!(w::GPUArraysCore.AbstractGPUVector{Tw}, t::Real, Ks::KrylovSubspace{T, U};
+          cache=nothing, expmethod=ExpMethodHigham2005Base()) -> w
+
+Non-allocating GPU-optimized version of `expv` that uses precomputed Krylov subspace `Ks`.
+
+# Arguments
+
+  - `w`: Pre-allocated GPU output vector for the result
+  - `t`: Time parameter for the exponential
+  - `Ks`: Precomputed Krylov subspace containing the approximation data
+  - `cache`: Optional pre-allocated workspace matrix
+  - `expmethod`: Method for computing the matrix exponential
+
+See also: [`expv`](@ref), [`expv!`](@ref)
+"""
 function ExponentialUtilities.expv!(w::GPUArraysCore.AbstractGPUVector{Tw},
         t::Real, Ks::KrylovSubspace{T, U};
         cache = nothing,
@@ -243,6 +294,19 @@ function phiv(t, A, b, k; cache = nothing, correct = false, errest = false,
     w = Matrix{eltype(b)}(undef, length(b), k + 1)
     phiv!(w, t, Ks, k; cache = cache, correct = correct, errest = errest)
 end
+"""
+    phiv(t, Ks::KrylovSubspace{T, U}, k; kwargs...) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b]
+
+Compute the matrix-phi-vector products using a precomputed Krylov subspace `Ks`.
+
+# Arguments
+
+  - `t`: Time parameter for the phi functions
+  - `Ks`: Precomputed Krylov subspace containing the approximation data
+  - `k`: Order of phi functions to compute (k >= 1)
+
+See also: [`phiv!`](@ref), [`arnoldi`](@ref)
+"""
 function phiv(t, Ks::KrylovSubspace{T, U}, k; kwargs...) where {T, U}
     n = size(getV(Ks), 1)
     w = Matrix{T}(undef, n, k + 1)