`ExponentialUtilities.Stegr.StegrWork`

— Method`StegrWork(T, n)`

Allocate work arrays for diagonalization of real-symmetric tridiagonal matrices of sizes up to `n`

×`n`

.

`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.ExpMethodDiagonalization`

— Type`ExpMethodDiagonalization(enforce_real=true)`

Matrix exponential method corresponding to the diagonalization with `eigen`

possibly by removing imaginary part introduced by the numerical approximation.

`ExponentialUtilities.ExpMethodGeneric`

— Type```
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.ExpMethodHigham2005`

— Type```
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.ExpMethodHigham2005Base`

— Type`ExpMethodHigham2005Base()`

The same as `ExpMethodHigham2005`

but follows `Base.exp`

closer.

`ExponentialUtilities.ExpMethodNative`

— Type`ExpMethodNative()`

Matrix exponential method corresponding to calling `Base.exp`

.

`ExponentialUtilities.KrylovSubspace`

— Type`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.alloc_mem`

— Method`cache=alloc_mem(A,method)`

Pre-allocates memory associated with matrix exponential function `method`

and matrix `A`

. To be used in combination with `exponential!`

.

`ExponentialUtilities.arnoldi!`

— Method`arnoldi!(Ks,A,b[;tol,m,opnorm,iop,init]) -> Ks`

Non-allocating version of `arnoldi`

.

`ExponentialUtilities.arnoldi`

— Method`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.arnoldi_step!`

— Method`arnoldi_step!(j, iop, n, A, V, H)`

Take the `j`

'th step of the Lanczos iteration.

`ExponentialUtilities.coeff`

— Method`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.expv`

— Method`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_timestep!`

— Method`expv_timestep!(u,t,A,b[;kwargs]) -> u`

Non-allocating version of `expv_timestep`

.

`ExponentialUtilities.expv_timestep`

— Method`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.firststep!`

— Method`firststep!(Ks, V, H, b) -> nothing`

Compute the first step of Arnoldi or Lanczos iteration.

`ExponentialUtilities.firststep!`

— Method`firststep!(Ks, V, H, b, b_aug, t, mu, l) -> nothing`

Compute the first step of Arnoldi or Lanczos iteration of augmented system.

`ExponentialUtilities.kiops`

— Method`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.lanczos!`

— Method`lanczos!(Ks,A,b[;tol,m,opnorm]) -> Ks`

A variation of `arnoldi!`

that uses the Lanczos algorithm for Hermitian matrices.

`ExponentialUtilities.lanczos_step!`

— Method`lanczos_step!(j, m, n, A, V, H)`

Take the `j`

'th step of the Lanczos iteration.

`ExponentialUtilities.phi!`

— Method`phi!(out,A,k[;caches]) -> out`

Non-allocating version of `phi`

for non-diagonal matrix inputs.

`ExponentialUtilities.phi`

— Method`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 `Diagonal`

s

`ExponentialUtilities.phi`

— Method`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.phiv`

— Method`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 phi*0 through phi*k-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) -> [phi*0(tA)b phi*1(tA)b ... phi_k(tA)b][, errest]

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

`ExponentialUtilities.phiv_dense!`

— Method`phiv_dense!(w,A,v,k[;cache]) -> w`

Non-allocating version of `phiv_dense`

.

`ExponentialUtilities.phiv_dense`

— Method`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_timestep!`

— Method`phiv_timestep!(U,ts,A,B[;kwargs]) -> U`

Non-allocating version of `phiv_timestep`

.

`ExponentialUtilities.phiv_timestep`

— Method`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.

`ExponentialUtilities.realview`

— Method`realview(::Type, V) -> real view of `V``

`ExponentialUtilities.@diagview`

— Macro`@diagview(A,d) -> view of the `d`th diagonal of `A`.`

- 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).