Saddle-point and Hermitian quasi-definite systems

TriCG

Krylov.tricg — Function

(x, y, stats) = tricg(A, b::AbstractVector{FC}, c::AbstractVector{FC};
                      M=I, N=I, ldiv::Bool=false,
                      spd::Bool=false, snd::Bool=false,
                      flip::Bool=false, τ::T=one(T),
                      ν::T=-one(T), atol::T=√eps(T),
                      rtol::T=√eps(T), itmax::Int=0,
                      timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                      callback=workspace->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

(x, y, stats) = tricg(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)

TriCG can be warm-started from initial guesses x0 and y0 where kwargs are the same keyword arguments as above.

Given a matrix A of dimension m × n, TriCG solves the Hermitian linear system

[ τE    A ] [ x ] = [ b ]
[  Aᴴ  νF ] [ y ]   [ c ],

of size (n+m) × (n+m) where τ and ν are real numbers, E = M⁻¹ ≻ 0 and F = N⁻¹ ≻ 0. TriCG could breakdown if τ = 0 or ν = 0. It's recommended to use TriMR in these cases.

By default, TriCG solves Hermitian and quasi-definite linear systems with τ = 1 and ν = -1.

TriCG is based on the preconditioned orthogonal tridiagonalization process and its relation with the preconditioned block-Lanczos process.

The matrix

[ M   0 ]
[ 0   N ]

indicates the weighted norm in which residuals are measured, here denoted ‖·‖. It's the Euclidean norm when M and N are identity operators.

TriCG stops when itmax iterations are reached or when ‖rₖ‖ ≤ atol + ‖r₀‖ * rtol. atol is an absolute tolerance and rtol is a relative tolerance.

Interface

To easily switch between Krylov methods, use the generic interface krylov_solve with method = :tricg.

For an in-place variant that reuses memory across solves, see tricg!.

Input arguments

A: a linear operator that models a matrix of dimension m × n;
b: a vector of length m;
c: a vector of length n.

Optional arguments

x0: a vector of length m that represents an initial guess of the solution x;
y0: a vector of length n that represents an initial guess of the solution y.

Warm-starting is supported only when M and N are either I or the corresponding coefficient (τ or ν) is zero.

Keyword arguments

M: linear operator that models a Hermitian positive-definite matrix of size m used for centered preconditioning of the partitioned system;
N: linear operator that models a Hermitian positive-definite matrix of size n used for centered preconditioning of the partitioned system;
ldiv: define whether the preconditioners use ldiv! or mul!;
spd: if true, set τ = 1 and ν = 1 for Hermitian and positive-definite linear system;
snd: if true, set τ = -1 and ν = -1 for Hermitian and negative-definite linear systems;
flip: if true, set τ = -1 and ν = 1 for another known variant of Hermitian quasi-definite systems;
τ and ν: diagonal scaling factors of the partitioned Hermitian linear system;
atol: absolute stopping tolerance based on the residual norm;
rtol: relative stopping tolerance based on the residual norm;
itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
timemax: the time limit in seconds;
verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
callback: function or functor called as callback(workspace) that returns true if the Krylov method should terminate, and false otherwise;
iostream: stream to which output is logged.

Output arguments

x: a dense vector of length m;
y: a dense vector of length n;
stats: statistics collected on the run in a SimpleStats structure.

Reference

A. Montoison and D. Orban, TriCG and TriMR: Two Iterative Methods for Symmetric Quasi-Definite Systems, SIAM Journal on Scientific Computing, 43(4), pp. 2502–2525, 2021.

source

Krylov.tricg! — Function

workspace = tricg!(workspace::TricgWorkspace, A, b, c; kwargs...)
workspace = tricg!(workspace::TricgWorkspace, A, b, c, x0, y0; kwargs...)

In these calls, kwargs are keyword arguments of tricg.

See TricgWorkspace for instructions on how to create the workspace.

For a more generic interface, you can use krylov_workspace with method = :tricg to allocate the workspace, and krylov_solve! to run the Krylov method in-place.

source

Krylov.TricgWorkspace — Type

Workspace for the in-place methods tricg! and krylov_solve!.

The following outer constructors can be used to initialize this workspace:

workspace = TricgWorkspace(m, n, Sm, Sn)
workspace = TricgWorkspace(m, n, S)
workspace = TricgWorkspace(A, b)
workspace = TricgWorkspace(A, b, c)
workspace = TricgWorkspace(kc::KrylovConstructor{Sm,Sn})

m and n denote the dimensions of the linear operator A passed to the in-place methods. Since tricg supports rectangular linear operators, m and n can differ. Sm and Sn are the storage types of the workspace vectors of length m and n, respectively. If the same storage type can be used for both, a single type S may be provided, such as Vector{Float64}.

KrylovConstructor facilitates the allocation of vectors in the workspace if Sm(undef, m) and Sn(undef, n) are not available.

source

TriMR

Krylov.trimr — Function

(x, y, stats) = trimr(A, b::AbstractVector{FC}, c::AbstractVector{FC};
                      M=I, N=I, ldiv::Bool=false,
                      spd::Bool=false, snd::Bool=false,
                      flip::Bool=false, sp::Bool=false,
                      τ::T=one(T), ν::T=-one(T), atol::T=√eps(T),
                      rtol::T=√eps(T), itmax::Int=0,
                      timemax::Float64=Inf, verbose::Int=0, history::Bool=false,
                      callback=workspace->false, iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

(x, y, stats) = trimr(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)

TriMR can be warm-started from initial guesses x0 and y0 where kwargs are the same keyword arguments as above.

Given a matrix A of dimension m × n, TriMR solves the symmetric linear system

[ τE    A ] [ x ] = [ b ]
[  Aᴴ  νF ] [ y ]   [ c ],

of size (n+m) × (n+m) where τ and ν are real numbers, E = M⁻¹ ≻ 0, F = N⁻¹ ≻ 0. TriMR handles saddle-point systems (τ = 0 or ν = 0) and adjoint systems (τ = 0 and ν = 0) without any risk of breakdown.

By default, TriMR solves symmetric and quasi-definite linear systems with τ = 1 and ν = -1.

TriMR is based on the preconditioned orthogonal tridiagonalization process and its relation with the preconditioned block-Lanczos process.

The matrix

[ M   0 ]
[ 0   N ]

indicates the weighted norm in which residuals are measured, here denoted ‖·‖. It's the Euclidean norm when M and N are identity operators.

TriMR stops when itmax iterations are reached or when ‖rₖ‖ ≤ atol + ‖r₀‖ * rtol. atol is an absolute tolerance and rtol is a relative tolerance.

Interface

To easily switch between Krylov methods, use the generic interface krylov_solve with method = :trimr.

For an in-place variant that reuses memory across solves, see trimr!.

Input arguments

A: a linear operator that models a matrix of dimension m × n;
b: a vector of length m;
c: a vector of length n.

Optional arguments

x0: a vector of length m that represents an initial guess of the solution x;
y0: a vector of length n that represents an initial guess of the solution y.

Warm-starting is supported only when M and N are either I or the corresponding coefficient (τ or ν) is zero.

Keyword arguments

M: linear operator that models a Hermitian positive-definite matrix of size m used for centered preconditioning of the partitioned system;
N: linear operator that models a Hermitian positive-definite matrix of size n used for centered preconditioning of the partitioned system;
ldiv: define whether the preconditioners use ldiv! or mul!;
spd: if true, set τ = 1 and ν = 1 for Hermitian and positive-definite linear system;
snd: if true, set τ = -1 and ν = -1 for Hermitian and negative-definite linear systems;
flip: if true, set τ = -1 and ν = 1 for another known variant of Hermitian quasi-definite systems;
sp: if true, set τ = 1 and ν = 0 for saddle-point systems;
τ and ν: diagonal scaling factors of the partitioned Hermitian linear system;
atol: absolute stopping tolerance based on the residual norm;
rtol: relative stopping tolerance based on the residual norm;
itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
timemax: the time limit in seconds;
verbose: additional details can be displayed if verbose mode is enabled (verbose > 0). Information will be displayed every verbose iterations;
history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
callback: function or functor called as callback(workspace) that returns true if the Krylov method should terminate, and false otherwise;
iostream: stream to which output is logged.

Output arguments

x: a dense vector of length m;
y: a dense vector of length n;
stats: statistics collected on the run in a SimpleStats structure.

Reference

A. Montoison and D. Orban, TriCG and TriMR: Two Iterative Methods for Symmetric Quasi-Definite Systems, SIAM Journal on Scientific Computing, 43(4), pp. 2502–2525, 2021.

source

Krylov.trimr! — Function

workspace = trimr!(workspace::TrimrWorkspace, A, b, c; kwargs...)
workspace = trimr!(workspace::TrimrWorkspace, A, b, c, x0, y0; kwargs...)

In these calls, kwargs are keyword arguments of trimr.

See TrimrWorkspace for instructions on how to create the workspace.

For a more generic interface, you can use krylov_workspace with method = :trimr to allocate the workspace, and krylov_solve! to run the Krylov method in-place.

source

Krylov.TrimrWorkspace — Type

Workspace for the in-place methods trimr! and krylov_solve!.

The following outer constructors can be used to initialize this workspace:

workspace = TrimrWorkspace(m, n, Sm, Sn)
workspace = TrimrWorkspace(m, n, S)
workspace = TrimrWorkspace(A, b)
workspace = TrimrWorkspace(A, b, c)
workspace = TrimrWorkspace(kc::KrylovConstructor{Sm,Sn})

m and n denote the dimensions of the linear operator A passed to the in-place methods. Since trimr supports rectangular linear operators, m and n can differ. Sm and Sn are the storage types of the workspace vectors of length m and n, respectively. If the same storage type can be used for both, a single type S may be provided, such as Vector{Float64}.

KrylovConstructor facilitates the allocation of vectors in the workspace if Sm(undef, m) and Sn(undef, n) are not available.

source

USYMLQR

Krylov.usymlqr — Function

(x, y, stats) = usymlqr(A, b::AbstractVector{FC}, c::AbstractVector{FC};
                        ls::Bool=true, ln::Bool=true, ldiv::Bool=false,
                        atol::T=√eps(T), rtol::T=√eps(T), itmax::Int=0,
                        timemax::Float64=Inf, verbose::Int=0,
                        history::Bool=false, callback=workspace->false,
                        iostream::IO=kstdout)

T is an AbstractFloat such as Float32, Float64 or BigFloat. FC is T or Complex{T}.

(x, y, stats) = usymlqr(A, b, c, x0::AbstractVector, y0::AbstractVector; kwargs...)

USYMLQR can be warm-started from initial guesses x0 and y0 where kwargs are the same keyword arguments as above.

Solve the symmetric saddle-point system

[ I   A ] [ x ] = [ b ]
[ Aᴴ    ] [ y ]   [ c ]

by way of the Saunders-Simon-Yip tridiagonalization using USYMLQ and USYMQR methods.

The method solves the least-squares problem when ls = true

[ I   A ] [ r ] = [ b ]
[ Aᴴ    ] [ s ]   [ 0 ]

and the least-norm problem when ln = true

[ I   A ] [ w ] = [ 0 ]
[ Aᴴ    ] [ z ]   [ c ]

and simply adds the solutions.

Interface

To easily switch between Krylov methods, use the generic interface krylov_solve with method = :usymlqr.

For an in-place variant that reuses memory across solves, see usymlqr!.

Input arguments

A: a linear operator that models a matrix of dimension m × n;
b: a vector of length m;
c: a vector of length n.

Optional arguments

x0: a vector of length m that represents an initial guess of the solution x;
y0: a vector of length n that represents an initial guess of the solution y.

Keyword arguments

ls: define whether the least-squares problem is solved;
ln: define whether the least-norm problem is solved;
ldiv: define whether the preconditioners use ldiv! or mul!;
atol: absolute stopping tolerance based on the residual norm;
rtol: relative stopping tolerance based on the residual norm;
itmax: the maximum number of iterations. If itmax=0, the default number of iterations is set to m+n;
timemax: the time limit in seconds;
verbose: additional details can be kdisplayed if verbose mode is enabled (verbose > 0). Information will be kdisplayed every verbose iterations;
history: collect additional statistics on the run such as residual norms, or Aᴴ-residual norms;
callback: function or functor called as callback(workspace) that returns true if the Krylov method should terminate, and false otherwise;
iostream: stream to which output is logged.

Output arguments

x: a dense vector of length m;
y: a dense vector of length n;
stats: statistics collected on the run in a SimpleStats structure.

References

M. A. Saunders, H. D. Simon, and E. L. Yip, Two Conjugate-Gradient-Type Methods for Unsymmetric Linear Equations, SIAM Journal on Numerical Analysis, 25(4), pp. 927–940, 1988.
A. Buttari, D. Orban, D. Ruiz and D. Titley-Peloquin, A tridiagonalization method for symmetric saddle-point and quasi-definite systems, SIAM Journal on Scientific Computing, 41(5), pp. 409–432, 2019.

source

Krylov.usymlqr! — Function

workspace = usymlqr!(workspace::UsymlqrWorkspace, A, b, c; kwargs...)
workspace = usymlqr!(workspace::UsymlqrWorkspace, A, b, c, x0, y0; kwargs...)

In these calls, kwargs are keyword arguments of usymlqr.

See UsymlqrWorkspace for instructions on how to create the workspace.

For a more generic interface, you can use krylov_workspace method = :usymlqr to allocate the workspace, and krylov_solve! to run the Krylov method in-place.

source

Krylov.UsymlqrWorkspace — Type

Workspace for the in-place methods usymlqr! and krylov_solve!.

The following outer constructors can be used to initialize this workspace:

workspace = UsymlqrWorkspace(m, n, Sm, Sn)
workspace = UsymlqrWorkspace(m, n, S)
workspace = UsymlqrWorkspace(A, b)
workspace = UsymlqrWorkspace(A, b, c)
workspace = UsymlqrWorkspace(kc::KrylovConstructor{Sm,Sn})

m and n denote the dimensions of the linear operator A passed to the in-place methods. Since usymlqr supports rectangular linear operators, m and n can differ. Sm and Sn are the storage types of the workspace vectors of length m and n, respectively. If the same storage type can be used for both, a single type S may be provided, such as Vector{Float64}.

KrylovConstructor facilitates the allocation of vectors in the workspace if Sm(undef, m) and Sn(undef, n) are not available.

source