arnoldi Interface

public interface arnoldi

Description

Given a square linear operator $A$, find matrices $X$ and $H$ such that

$$A X_k = X_k H_k + H_{k+1,k} X_{k+1} E_k^T,$$

where $X$ is an orthogonal basis, $H$ is upper Hessenberg, and for block Arnoldi with block size $p$, $H_{k+1,k}$ is a $p \times p$ matrix (reducing to a scalar $h_{k+1,k}$ when $p = 1$).

Algorithmic Features

• The operator $A$ only needs to be accessed through matrix-vector products.
• Constructs an orthonormal Krylov basis $X$ via the Gram-Schmidt process.
• Constructs an upper Hessenberg matrix $H$ whose eigenvalues approximates those of $A$.
• Checks for convergence and invariant subspaces.

References

• Y. Saad. "Iterative methods for sparse linear systems", SIAM 2nd edition, 2003. see Chapter 6.3: Arnoldi's method.

Syntax

    call arnoldi(A, X, H, info [, kstart] [, kend] [, tol] [, transpose] [, blksize])

Arguments

• A : Linear operator derived from one the base types provided by the AbstractLinops module. The operator needs to be square, i.e. the dimension of its domain and co-domain is the same. It is an intent(inout) argument.

• X : Array of types derived from one the base types provided by the AbstractVectors module. It needs to be consistent with the type of A. On exit, it contains the the computed Krylov vectors. The first entry X(1) is the starting vector for the Arnoldi factorization. Additionally, the maximum number of Arnoldi steps is equal to size(X) - 1. It is an intent(inout) argument.

-H : real or complex rank-2 array. On exit, it contains the $(k_{max} \times p + p) \times (k_{max} \times p)$ upper Hessenberg matrix computed from the Arnoldi factorization, where $k_{max} = size(X)/p - 1$ is the maximum number of Arnoldi steps and $p$ is the block size (default 1). It is an intent(inout) argument.

-info : integer variable. It is the LightKrylov information flag. On exit, if info > 0, the Arnoldi factorization experienced a lucky breakdown. The array of Krylov vectors X spans an $A$-invariant subpsace of dimension info. If info < 0, an input argument is illegal: info = -3 indicates that H is not large enough: it must have leading dimensions at least (size(X), p*kdim) with kdim = (size(X) - p) / p and block size p; info = -5 indicates an illegal kstart; info = -6 an illegal kend; info = -7 a negative tolerance; info = -9 an illegal blksize (blksize $\leq$ 0 or size(X) not divisible by blksize).

• kstart (optional) : integer value determining the index of the first Arnoldi step to be computed. By default, kstart = 1. It is an optional intent(in) argument.

• kend (optional) : integer value determining the index of the last Arnoldi step to be computed. It is an optional intent(in) argument. By default, kend = size(X) - 1.

• tol (optional) : Numerical tolerance below which a subspace is considered to be $A$-invariant. It is an optional intent(in) argument. By default tol = atol_sp or tol = atol_rp depending on the kind of A.

• transpose (optional) : logical flag determining whether the Arnoldi factorization is applied to $A$ or $A^H$. It is an optional intent(in) argument. Default transpose = .false.

• blksize (optional) : integer value determining the dimension of a block for the block Arnoldi factorization. It is an optional intent(in) argument. Default is blksize=1. The block size $p$ must divide size(X) evenly, and the maximum number of Arnoldi steps is $k_{max} = (size(X) - p) / p$.

Subroutines