LinearAlgebra.LAPACK.stegr!Method
stegr!(α, β, sw)

Diagonalize the real-symmetric tridiagonal matrix with α on the diagonal and β on the super-/subdiagonal, using the workspaces allocated in sw.

ExponentialUtilities.ExpMethodDiagonalizationType
ExpMethodDiagonalization(enforce_real=true)

Matrix exponential method corresponding to the diagonalization with eigen possibly by removing imaginary part introduced by the numerical approximation.

ExponentialUtilities.ExpMethodGenericType
struct ExpMethodGeneric{T}
ExpMethodGeneric()=ExpMethodGeneric{Val{13}}();

Generic exponential implementation of the method ExpMethodHigham2005, for any exp argument x for which the functions LinearAlgebra.opnorm, +, *, ^, and / (including addition with UniformScaling objects) are defined. The type T is used to adjust the number of terms used in the Pade approximants at compile time.

See "The Scaling and Squaring Method for the Matrix Exponential Revisited" by Higham, Nicholas J. in 2005 for algorithm details.

ExponentialUtilities.ExpMethodHigham2005Type
ExpMethodHigham2005(A::AbstractMatrix);
ExpMethodHigham2005(b::Bool=true);

Computes the matrix exponential using the algorithm Higham, N. J. (2005). "The scaling and squaring method for the matrix exponential revisited." SIAM J. Matrix Anal. Appl.Vol. 26, No. 4, pp. 1179–1193" based on generated code. If a matrix is specified, balancing is determined automatically.

ExponentialUtilities.KrylovSubspaceType
KrylovSubspace{T}(n,[maxiter=30]) -> Ks

Constructs an uninitialized Krylov subspace, which can be filled by arnoldi!.

The dimension of the subspace, Ks.m, can be dynamically altered but should be smaller than maxiter, the maximum allowed arnoldi iterations.

getV(Ks) -> V
getH(Ks) -> H

Access methods for the (extended) orthonormal basis V and the (extended) Gram-Schmidt coefficients H. Both methods return a view into the storage arrays and has the correct dimensions as indicated by Ks.m.

resize!(Ks, maxiter) -> Ks

Resize Ks to a different maxiter, destroying its contents.

This is an expensive operation and should be used scarcely.

ExponentialUtilities.arnoldiMethod
arnoldi(A,b[;m,tol,opnorm,iop]) -> Ks

Performs m arnoldi iterations to obtain the Krylov subspace K_m(A,b).

The n x (m + 1) basis vectors getV(Ks) and the (m + 1) x m upper Hessenberg matrix getH(Ks) are related by the recurrence formula

v_1=b,\quad Av_j = \sum_{i=1}^{j+1}h_{ij}v_i\quad(j = 1,2,\ldots,m)

iop determines the length of the incomplete orthogonalization procedure [1]. The default value of 0 indicates full Arnoldi. For symmetric/Hermitian A, iop will be ignored and the Lanczos algorithm will be used instead.

Refer to KrylovSubspace for more information regarding the output.

Happy-breakdown occurs whenever norm(v_j) < tol * opnorm, in this case, the dimension of Ks is smaller than m.

ExponentialUtilities.coeffMethod
coeff(::Type,α)

Helper functions that returns the real part if that is what is required (for Hermitian matrices), otherwise returns the value as-is.

ExponentialUtilities.expT!Method
expT!(α, β, t, cache)

Calculate the subspace exponential exp(t*T) for a tridiagonal subspace matrix T with α on the diagonal and β on the super-/subdiagonal, diagonalizing via stegr!.

ExponentialUtilities.exponential!Function
E=exponential!(A,[method [cache]])

Computes the matrix exponential with the method specified in method. The contents of A are modified, allowing for fewer allocations. The method parameter specifies the implementation and implementation parameters, e.g. ExpMethodNative, ExpMethodDiagonalization, ExpMethodGeneric, ExpMethodHigham2005. Memory needed can be preallocated and provided in the parameter cache such that the memory can be recycled when calling exponential! several times. The preallocation is done with the command alloc_mem: cache=alloc_mem(A,method).

Example

julia> A = randn(50, 50);

julia> Acopy = B * 2;

julia> method = ExpMethodHigham2005();

julia> cache = alloc_mem(A, method); # Main allocation done here

julia> E1 = exponential!(A, method, cache) # Very little allocation here

julia> E2 = exponential!(B, method, cache) # Very little allocation here

ExponentialUtilities.expv!Method
expv!(w, t, A, b, Ks, cache)

Alternative interface for calculating the action of exp(t*A) on the vector b, storing the result in w. The Krylov iteration is terminated when an error estimate for the matrix exponential in the generated subspace is below the requested tolerance. Ks is a KrylovSubspace and typeof(cache)<:HermitianSubspaceCache, the exact type decides which algorithm is used to compute the subspace exponential.

ExponentialUtilities.expv!Method
expv!(w,t,Ks[;cache]) -> w

Non-allocating version of expv that uses precomputed Krylov subspace Ks.

ExponentialUtilities.expvMethod
expv(t,A,b; kwargs) -> exp(tA)b

Compute the matrix-exponential-vector product using Krylov.

A Krylov subspace is constructed using arnoldi and exp! is called on the Hessenberg matrix. Consult arnoldi for the values of the keyword arguments. An alternative algorithm, where an error estimate generated on-the-fly is used to terminate the Krylov iteration, can be employed by setting the kwarg mode=:error_estimate.

expv(t,Ks; cache) -> exp(tA)b

Compute the expv product using a pre-constructed Krylov subspace.

ExponentialUtilities.expv_timestepMethod
expv_timestep(ts,A,b[;adaptive,tol,kwargs...]) -> U

Evaluates the matrix exponentiation-vector product using time stepping

$$$u = \exp(tA)b$$$

ts is an array of time snapshots for u, with U[:,j] ≈ u(ts[j]). ts can also be just one value, in which case only the end result is returned and U is a vector.

The time stepping formula of Niesen & Wright is used [1]. If the time step tau is not specified, it is chosen according to (17) of Niesen & Wright. If adaptive==true, the time step and Krylov subspace size adaptation scheme of Niesen & Wright is used, the relative tolerance of which can be set using the keyword parameter tol. The delta and gamma parameters of the adaptation scheme can also be adjusted.

Set verbose=true to print out the internal steps (for debugging). For the other keyword arguments, consult arnoldi and phiv, which are used internally.

Note that this function is just a special case of phiv_timestep with a more intuitive interface (vector b instead of a n-by-1 matrix B).

ExponentialUtilities.kiopsMethod
kiops(tstops, A, u; kwargs...) -> (w, stats)

Evaluate a linear combination of the $φ$ functions evaluated at $tA$ acting on vectors from $u$, that is

$$$w(i) = φ_0(t[i] A) u[:, 1] + φ_1(t[i] A) u[:, 2] + φ_2(t[i] A) u[:, 3] + ...$$$

The size of the Krylov subspace is changed dynamically during the integration. The Krylov subspace is computed using the incomplete orthogonalization method.

Arguments:

• tstops - Array of tstop
• A - the matrix argument of the $φ$ functions
• u - the matrix with columns representing the vectors to be multiplied by the $φ$ functions

Keyword arguments:

• tol - the convergence tolerance required (default: 1e-7)
• mmin, mmax - let the Krylov size vary between mmin and mmax (default: 10, 128)
• m - an estimate of the appropriate Krylov size (default: mmin)
• iop - length of incomplete orthogonalization procedure (default: 2)
• ishermitian - whether $A$ is Hermitian (default: ishermitian(A))
• opnorm - the operator norm of $A$ (default: opnorm(A, Inf))
• task1 - if true, divide the result by 1/T^p

Returns:

• w - the linear combination of the $φ$ functions evaluated at $tA$ acting on the vectors from $u$
• stats[1] - number of substeps
• stats[2] - number of rejected steps
• stats[3] - number of Krylov steps
• stats[4] - number of matrix exponentials
• stats[5] - the Krylov size of the last substep

n is the size of the original problem p is the highest index of the $φ$ functions

References:

• Gaudreault, S., Rainwater, G. and Tokman, M., 2018. KIOPS: A fast adaptive Krylov subspace solver for exponential integrators. Journal of Computational Physics. Based on the PHIPM and EXPMVP codes (http://www1.maths.leeds.ac.uk/~jitse/software.html). https://gitlab.com/stephane.gaudreault/kiops.
• Niesen, J. and Wright, W.M., 2011. A Krylov subspace method for option pricing. SSRN 1799124
• 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
ExponentialUtilities.phiMethod
phi(A,k[;cache]) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the matrix phi functions for all orders up to k. k >= 1.

The phi functions are defined as

$$$\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{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 Diagonals

ExponentialUtilities.phiMethod
phi(z,k[;cache]) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the scalar phi functions for all orders up to k.

The phi functions are defined as

$$$\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{z}$$$

Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156. Theorem 1).

ExponentialUtilities.phiv!Method
phiv!(w,t,Ks,k[;cache,correct,errest]) -> w[,errest]

Non-allocating version of 'phiv' that uses precomputed Krylov subspace Ks.

ExponentialUtilities.phivMethod
phiv(t,A,b,k;correct,kwargs) -> [phi_0(tA)b phi_1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using Krylov. k >= 1.

The phi functions are defined as

$$$\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{z}$$$

A Krylov subspace is constructed using arnoldi and phiv_dense is called on the Hessenberg matrix. If correct=true, then phi0 through phik-1 are updated using the last Arnoldi vector v_m+1 [1]. If errest=true then an additional error estimate for the second-to-last phi is also returned. For the additional keyword arguments, consult arnoldi.

phiv(t,Ks,k;correct,kwargs) -> [phi0(tA)b phi1(tA)b ... phi_k(tA)b][, errest]

Compute the matrix-phi-vector products using a pre-constructed Krylov subspace.

ExponentialUtilities.phiv_denseMethod
phiv_dense(A,v,k[;cache]) -> [phi_0(A)v phi_1(A)v ... phi_k(A)v]

Compute the matrix-phi-vector products for small, dense A. k >= 1.

The phi functions are defined as

$$$\varphi_0(z) = \exp(z),\quad \varphi_{k+1}(z) = \frac{\varphi_k(z) - 1}{z}$$$

Instead of using the recurrence relation, which is numerically unstable, a formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software package for computing matrix exponentials. ACM Transactions on Mathematical Software (TOMS), 24(1), 130-156. Theorem 1).

ExponentialUtilities.phiv_timestepMethod
phiv_timestep(ts,A,B[;adaptive,tol,kwargs...]) -> U

Evaluates the linear combination of phi-vector products using time stepping

$$$u = \varphi_0(tA)b_0 + t\varphi_1(tA)b_1 + \cdots + t^p\varphi_p(tA)b_p$$$

ts is an array of time snapshots for u, with U[:,j] ≈ u(ts[j]). ts can also be just one value, in which case only the end result is returned and U is a vector.

The time stepping formula of Niesen & Wright is used [1]. If the time step tau is not specified, it is chosen according to (17) of Niesen & Wright. If adaptive==true, the time step and Krylov subspace size adaptation scheme of Niesen & Wright is used, the relative tolerance of which can be set using the keyword parameter tol. The delta and gamma parameters of the adaptation scheme can also be adjusted.

When encountering a happy breakdown in the Krylov subspace construction, the time step is set to the remainder of the time interval since time stepping is no longer necessary.

Set verbose=true to print out the internal steps (for debugging). For the other keyword arguments, consult arnoldi and phiv`, which are used internally.

• 1Koskela, A. (2015). Approximating the matrix exponential of an advection-diffusion operator using the incomplete orthogonalization method. In Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353). Springer, Cham.
• 1Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631.
• 1Niesen, J., & Wright, W. (2009). A Krylov subspace algorithm for evaluating the φ-functions in exponential integrators. arXiv preprint arXiv:0907.4631. Formula (10).