StegrWork(T, n)

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

stegr!(α, β, sw)

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

ExpMethodDiagonalization(enforce_real=true)

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

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.

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.

ExpMethodHigham2005Base()

The same as ExpMethodHigham2005 but follows Base.exp closer.

ExpMethodNative()

Matrix exponential method corresponding to calling Base.exp.

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.

cache=alloc_mem(A,method)

Pre-allocates memory associated with matrix exponential function method and matrix A. To be used in combination with exponential!.

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

Non-allocating version of arnoldi.

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.

arnoldi_step!(j, iop, n, A, V, H)

Take the j 'th step of the Lanczos iteration.

coeff(::Type,α)

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

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

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

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.

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

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

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.

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

Non-allocating version of expv_timestep.

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

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

Compute the first step of Arnoldi or Lanczos iteration.

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

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

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

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

A variation of arnoldi! that uses the Lanczos algorithm for Hermitian matrices.

lanczos_step!(j, m, n, A, V, H)

Take the j'th step of the Lanczos iteration.

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

Non-allocating version of phi for non-diagonal matrix inputs.

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

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

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

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

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.

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

Non-allocating version of phiv_dense.

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

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

Non-allocating version of phiv_timestep.

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.

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

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