LightKrylov_BaseKrylov Module

Description

Given an array $X$ and a permutation vector $p$, this function computes in-place the column-permuted matrix

$$X = X P^{-1}$$

where $P$ is the column-permutation matrix constructed from the permutation vector $p$ and $P^{-1}$ its inverse.

Syntax

    call apply_inverse_permutation_matrix(X, perm)

Arguments

Q : Array of vectors derived from the base types defined in the AbstractVectors module. On entry, it is the original array. On exit, it contains the column-permuted version computed in-place. It is an intent(inout) argument.

perm : Rank-1 array of integer corresponding to the desired permutation vector. It is an intent(in) argument.

Description

Given an array $X$ and a permutation vector $p$, this function computes in-place the column-permuted matrix

$$X = X P$$

where $P$ is the column-permutation matrix constructed from the permutation vector $p$.

Syntax

    call apply_permutation_matrix(X, perm)

Arguments

Q : Array of vectors derived from the base types defined in the AbstractVectors module. On entry, it is the original array. On exit, it contains the column-permuted version computed in-place. It is an intent(inout) argument.

perm : Rank-1 array of integer corresponding to the desired permutation vector. It is an intent(in) argument.

Description

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

$$AX_k = X_k H_k + h_{k+1, k} x_{k+1} e_k^T,$$

where $X$ is an orthogonal basis and $H$ is upper Hessenberg.

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(in) 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+1) \times k$ upper Hessenberg matrix computed from the Arnoldi factorization. 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.

kstart (optional): integer value determining the index of the first Arnoldi step to be computed. By default, kstart = 1.

kend (optional): integer value determining the index of the last Arnoldi step to be computed. By default, kend = size(X) - 1.

tol (optional): Numerical tolerance below which a subspace is considered to be $A$-invariant. 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$. Default transpose = .false.

blksize (optional): integer value determining the dimension of a block for the block Arnoldi factorization. Default is blksize=1.

Description

Given a general linear operator $A$, find matrices $U$, $V$ and $B$ such that

$$\begin{aligned} AV_k & = U_{k+1} B_k, \\ A^H U_{k+1} & = V_k B_k^T + b_{k+1} v_{k+1} e_{k+1}^T \end{aligned}$$

where $U$ and $V$ are orthogonal bases for the column span and row span of $A$, respectively, and $B$ is a bidiagonal matrix.

Algorithmic Features

• The operator $A$ only needs to be accessed through matrix-vector products.
• Constructs an orthonormal Krylov basis $U$ for the column span of $A$.
• Constructs an orthonormal Krylov basis $V$ for the row span of $A$.
• Constructs a bidiagonal matrix $B$ whose singular values approximates those of $A$.
• Checks for convergence and invariant subspaces.

References

• R. M. Larsen. "Lanczos bidiagonalization with partial reorthogonalization." Technical Report, 1998. (PDF)

Syntax

    call bidiagonalization(A, U, V, B, info [, kstart] [, kend] [, tol])

Arguments

A : Linear operator derived from one the base types provided by the AbstractLinops module. It is an intent(in) argument.

U : 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 for the column span of A. The first entry U(1) is the starting vector for the Lanczos factorization. Additionally, the maximum number of Lanczos steps is equal to size(X) - 1. It is an intent(inout) argument.

V : 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 for the row span of A. It is an intent(inout) argument.

B : real or complex rank-2 array. On exit, it contains the $(k+1) \times k$ bidiagonal matrix computed from the Lanczos factorization. It is an intent(inout) argument.

info : integer variable. It is the LightKrylov information flag. On exit, if info > 0, the Lanczos factorization experienced a lucky breakdown.

kstart (optional): integer value determining the index of the first Lanczos step to be computed. By default, kstart = 1.

kend (optional): integer value determining the index of the last Lanczos step to be computed. By default, kend = size(X) - 1.

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

Description

Given an array $X$ of abstract_vector and an abstract_vector (or array of abstract_vectors) $y$, this subroutine returns a modified vector $y$ orthogonal to all columns of $X$, i.e.

$$y \leftarrow \left( I - XX^H \right) y,$$

using a double Gram-Schmidt process. On exit, $y$ is orthogonal to $X$ albeit does not have unit-norm. Note moreover that $X$ is assumed to be an orthonormal set of vectors. The function can also return the projection coefficients $\beta = X^H y$.

Syntax

    call double_gram_schmidt_step(y, X, info [, if_chk_orthonormal] [, beta])

Arguments

y : abstract_vector (or array of abstract_vector) that needs to be orthogonalize in-place against $X$.

X : Array of abstract_vector against which $y$ needs to be orthogonalized. Note the function assumes that $X$ is an orthonormal set of vectors, i.e. $X^H X = I$. If it this is not the case, the result are meaningless.

info : integer Information flag.

if_chk_orthonormal (optional) : logical flag (default .true.) to check whether $X$ is an orthonormal set of vectors or not. If the orthonormality returns .false., the function throws an error. Note that this check is however computationally expensive and can be disable for the sake of performances.

beta (optional) : real or complex array containing the coefficients $\beta = X^H y$.

Description

Utility function to initialize a basis for a Krylov subspace.

Syntax

    call initialize_krylov_subspace(X [, X0])

Arguments

X : Array of vectors that needs to be initialized. It is an intent(inout) argument. Note that the first action in the subroutine is call zero_basis(X), effectively zeroing-out any data stored.

X0 (optional) : Collection of vectors which will form the first few Krylov vectors. Note that X0 need not be an orthonormal basis as this subroutine includes a call qr(X0).

Description

Utility function returning a logical .true. if the set of vectors stored in $X$ form an orthonormal set of vectors and .false. otherwise.

Syntax

    out = is_orthonormal(X)

Arguments

X : Array of derived types extended from the base types provided in the AbstractVectors module.

Description

Given a partial Krylov decomposition

$$AX_k = X_k H_k + h_{k+1, k} x_{k+1} e_k^H,$$

this subroutine implements the Krylov-Schur restarting strategy proposed by Stewart [1].

References

• G. W. Stewart. "A Krylov-Schur algorithm for large eigenproblems". SIAM Journal on Matrix Analysis and Applications, vol 23 (3), 2002.

Syntax

    call krylov_schur(n, X, H, select_eigs)

Arguments

n : Number of selected eigenvalues moved to the upper left-block of the Schur matrix. It is an intent(out) argument.

X : On entry, array of abstract_vector computed using the Arnoldi process. On exit, the first n columns form an orthonormal basis for the eigenspace associated with eigenvalues moved to the upper left-block of the Schur matrix. It is an intent(inout) argument.

H : On entry, real of complex upper Hessenberg matrix computed using the Arnoldi process. On exit, the leading $n \times n$ block contains the $S_{11}$ block of the re-ordered Schur matrix containing the selected eigenvalues. It is an intent(inout) argument.

select_eigs : Procedure to select which eigenvalues to move in the upper-left block. It is an intent(inout) argument.

Description

Given a symmetric or Hermitian linear operator $A$, find matrices $X$ and $T$ such that

$$AX_k = X_k T_k + t_{k+1, k} x_{k+1} e_k^T,$$

where $X$ is an orthogonal basis and $T$ is symmetric tridiagonal.

Algorithmic Features

• The operator $A$ only needs to be accessed through matrix-vector products.
• Constructs an orthonormal Krylov basis $X$ via the Lanczos process with full reorthogonalization.
• Constructs a symmetric tridiagonal matrix $T$ 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.6: The symmetric Lanczos algorithm.

Syntax

    call lanczos(A, X, T, info [, kstart] [, kend] [, tol])

Arguments

A : Symmetric or Hermitian linear operator derived from one the base types provided by the AbstractLinops module. It is an intent(in) 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 Lanczos factorization. Additionally, the maximum number of Lanczos steps is equal to size(X) - 1. It is an intent(inout) argument.

T : real or complex rank-2 array. On exit, it contains the $(k+1) \times k$ symmetric tridiagonal matrix computed from the Arnoldi factorization. It is an intent(inout) argument.

info : integer variable. It is the LightKrylov information flag. On exit, if info > 0, the Lanczos factorization experienced a lucky breakdown. The array of Krylov vectors X spans an $A$-invariant subpsace of dimension info.

kstart (optional): integer value determining the index of the first Lanczos step to be computed. By default, kstart = 1.

kend (optional): integer value determining the index of the last Lanczos step to be computed. By default, kend = size(X) - 1.

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

Description

Given an array $X$ of vectors, it computes an orthonormal basis for its column-span using the double_gram_schmidt process. All computations are done in-place.

Syntax

    call orthonormalize_basis(X)

Arguments

X : Array of abstract_vector to orthonormalize. Note that this process is done in-place. It is an intent(inout) argument.

Description

Given an array $X$ of types derived from abstract_vector, it computes the in-place QR factorization of $X$, i.e.

$$X = QR,$$

where $Q$ is an orthonormal arrays of vectors such that $Q^H Q = I$ and $R$ is upper triangular. Note that it can also perform the QR factorization with column pivoting

$$XP = QR$$

where $P$ is a permutation matrix ensuring that the diagonal entries of $R$ have non-increasing absolute values. This amounts to using the pivoting QR as a rank-revealing factorization.

References

• G. H. Golub & C. F. Van Loan. "Matrix Computations". 4th edition, The John Hopkins University Press, 2013. See Chapter 5.2.8: Modified Gram-Schmidt algorithm.

Syntax

    call qr(Q [, R] [, perm], info [, tol])

Arguments

Q: Array of types derived from one of the base types provided in the AbstractVectors module. On entry, it contains the original array. On exit, it is overwritten by the orthogonal basis for its span. It is an intent(inout) argument.

R: real or complex rank-2 array. On exit, its contains the upper triangular matrix resulting from the QR factorization. It is an intent(out) argument.

perm (optional): Rank-1 array of integer corresponding to the indices of permuted columns. If perm is absent, the naive QR factorization is being computed.

info: integer information flag.

tol (optional): Numerical tolerance to determine whether two vectors are colinear or not. Default tol = atol_sp or tol = atol_dp.