krypy.linsys
 Linear Algebraic Systems Solver¶
The linsys module provides functions for the solution of linear algebraic systems.

class
krypy.linsys.
LinearSystem
(A, b, M=None, Minv=None, Ml=None, Mr=None, ip_B=None, normal=None, self_adjoint=False, positive_definite=False, exact_solution=None)¶ Bases:
object
Representation of a (preconditioned) linear system.
Represents a linear system
\[Ax=b\]or a preconditioned linear system
\[M M_l A M_r y = M M_l b \quad\text{with}\quad x=M_r y.\]Parameters:  A – a linear operator on \(\mathbb{C}^N\) (has to be
compatible with
get_linearoperator()
).  b – the right hand side in \(\mathbb{C}^N\), i.e.,
b.shape == (N, 1)
.  M – (optional) a selfadjoint and positive definite
preconditioner, linear operator on \(\mathbb{C}^N\) with respect
to the inner product defined by
ip_B
. This preconditioner changes the inner product to \(\langle x,y\rangle_M = \langle Mx,y\rangle\) where \(\langle \cdot,\cdot\rangle\) is the inner product defined by the parameterip_B
. Defaults to the identity.  Minv – (optional) the inverse of the preconditioner provided by
M
. This operator is needed, e.g., for orthonormalizing vectors for the computation of Ritz vectors in deflated methods.  Ml – (optional) left preconditioner, linear operator on \(\mathbb{C}^N\). Defaults to the identity.
 Mr – (optional) right preconditioner, linear operator on \(\mathbb{C}^N\). Defaults to the identity.
 ip_B – (optional) defines the inner product, see
inner()
.  normal – (bool, optional) Is \(M_l A M_r\) normal
in the inner product defined by
ip_B
? Defaults toFalse
.  self_adjoint – (bool, optional) Is \(M_l A M_r\) selfadjoint
in the inner product defined by
ip_B
?self_adjoint=True
also setsnormal=True
. Defaults toFalse
.  positive_definite – (bool, optional) Is \(M_l A M_r\)
positive (semi)definite with respect to the inner product defined by
ip_B
? Defaults toFalse
.  exact_solution – (optional) If an exact solution \(x\) is
known, it can be provided as a
numpy.array
withexact_solution.shape == (N,1)
. Then error norms can be computed (for debugging or research purposes). Defaults toNone
.

MMlb_norm
= None¶ Norm of the right hand side.
\[\M M_l b\_{M^{1}}\]

N
= None¶ Dimension \(N\) of the space \(\mathbb{C}^N\) where the linear system is defined.

get_ip_Minv_B
()¶ Returns the inner product that is implicitly used with the positive definite preconditioner
M
.

get_residual
(z, compute_norm=False)¶ Compute residual.
For a given \(z\in\mathbb{C}^N\), the residual
\[r = M M_l ( b  A z )\]is computed. If
compute_norm == True
, then also the absolute residual norm\[\ M M_l (bAz)\_{M^{1}}\]is computed.
Parameters:  z – approximate solution with
z.shape == (N, 1)
.  compute_norm – (bool, optional) pass
True
if also the norm of the residual should be computed.
 z – approximate solution with
 A – a linear operator on \(\mathbb{C}^N\) (has to be
compatible with

class
krypy.linsys.
Cg
(linear_system, **kwargs)¶ Bases:
krypy.linsys._KrylovSolver
Preconditioned CG method.
The preconditioned conjugate gradient method can be used to solve a system of linear algebraic equations where the linear operator is selfadjoint and positive definite. Let the following linear algebraic system be given:
\[M M_l A M_r y = M M_l b,\]where \(x=M_r y\) and \(M_l A M_r\) is selfadjoint and positive definite with respect to the inner product \(\langle \cdot,\cdot \rangle\) defined by
ip_B
. The preconditioned CG method then computes (in exact arithmetics!) iterates \(x_k \in x_0 + M_r K_k\) with \(K_k:= K_k(M M_l A M_r, r_0)\) such that\[\x  x_k\_A = \min_{z \in x_0 + M_r K_k} \x  z\_A.\]The Lanczos alorithm is used with the operator \(M M_l A M_r\) and the inner product defined by \(\langle x,y \rangle_{M^{1}} = \langle M^{1}x,y \rangle\). The initial vector for Lanczos is \(r_0 = M M_l (b  Ax_0)\)  note that \(M_r\) is not used for the initial vector.
Memory consumption is:
 if
store_arnoldi==False
: 3 vectors or 6 vectors if \(M\) is used.  if
store_arnoldi==True
: about maxiter+1 vectors for the Lanczos basis. If \(M\) is used the memory consumption is 2*(maxiter+1).
Caution: CG’s convergence may be delayed significantly due to roundoff errors, cf. chapter 5.9 in [LieS13].
All parameters of
_KrylovSolver
are valid in this solver. Note the restrictions onM
,Ml
,A
,Mr
andip_B
above.
static
operations
(nsteps)¶ Returns the number of operations needed for nsteps of CG
 if

class
krypy.linsys.
Minres
(linear_system, ortho='lanczos', **kwargs)¶ Bases:
krypy.linsys._KrylovSolver
Preconditioned MINRES method.
The preconditioned minimal residual method can be used to solve a system of linear algebraic equations where the linear operator is selfadjoint. Let the following linear algebraic system be given:
\[M M_l A M_r y = M M_l b,\]where \(x=M_r y\) and \(M_l A M_r\) is selfadjoint with respect to the inner product \(\langle \cdot,\cdot \rangle\) defined by
inner_product
. The preconditioned MINRES method then computes (in exact arithmetics!) iterates \(x_k \in x_0 + M_r K_k\) with \(K_k:= K_k(M M_l A M_r, r_0)\) such that\[\M M_l(b  A x_k)\_{M^{1}} = \min_{z \in x_0 + M_r K_k} \M M_l (b  A z)\_{M^{1}}.\]The Lanczos alorithm is used with the operator \(M M_l A M_r\) and the inner product defined by \(\langle x,y \rangle_{M^{1}} = \langle M^{1}x,y \rangle\). The initial vector for Lanczos is \(r_0 = M M_l (b  Ax_0)\)  note that \(M_r\) is not used for the initial vector.
Memory consumption is:
 if
store_arnoldi==False
: 3 vectors or 6 vectors if \(M\) is used.  if
store_arnoldi==True
: about maxiter+1 vectors for the Lanczos basis. If \(M\) is used the memory consumption is 2*(maxiter+1).
Caution: MINRES’ convergence may be delayed significantly or even stagnate due to roundoff errors, cf. chapter 5.9 in [LieS13].
In addition to the attributes described in
_KrylovSolver
, the following attributes are available in an instance of this solver:lanczos
: the Lanczos relation (an instance ofArnoldi
).
All parameters of
_KrylovSolver
are valid in this solver. Note the restrictions onM
,Ml
,A
,Mr
andip_B
above.
static
operations
(nsteps)¶ Returns the number of operations needed for nsteps of MINRES
 if

class
krypy.linsys.
Gmres
(linear_system, ortho='mgs', **kwargs)¶ Bases:
krypy.linsys._KrylovSolver
Preconditioned GMRES method.
The preconditioned generalized minimal residual method can be used to solve a system of linear algebraic equations. Let the following linear algebraic system be given:
\[M M_l A M_r y = M M_l b,\]where \(x=M_r y\). The preconditioned GMRES method then computes (in exact arithmetics!) iterates \(x_k \in x_0 + M_r K_k\) with \(K_k:= K_k(M M_l A M_r, r_0)\) such that
\[\M M_l(b  A x_k)\_{M^{1}} = \min_{z \in x_0 + M_r K_k} \M M_l (b  A z)\_{M^{1}}.\]The Arnoldi alorithm is used with the operator \(M M_l A M_r\) and the inner product defined by \(\langle x,y \rangle_{M^{1}} = \langle M^{1}x,y \rangle\). The initial vector for Arnoldi is \(r_0 = M M_l (b  Ax_0)\)  note that \(M_r\) is not used for the initial vector.
Memory consumption is about maxiter+1 vectors for the Arnoldi basis. If \(M\) is used the memory consumption is 2*(maxiter+1).
If the operator \(M_l A M_r\) is selfadjoint then consider using the MINRES method
Minres
.All parameters of
_KrylovSolver
are valid in this solver.
static
operations
(nsteps)¶ Returns the number of operations needed for nsteps of GMRES

static

class
krypy.linsys.
_KrylovSolver
(linear_system, x0=None, tol=1e05, maxiter=None, explicit_residual=False, store_arnoldi=False, dtype=None)¶ Prototype of a Krylov subspace method for linear systems.
Init standard attributes and perform checks.
All Krylov subspace solvers in this module are applied to a
LinearSystem
. The specific methods may impose further restrictions on the operatorsParameters:  linear_system – a
LinearSystem
.  x0 – (optional) the initial guess to use. Defaults to zero
vector. Unless you have a good reason to use a nonzero initial guess
you should use the zero vector, cf. chapter 5.8.3 in Liesen,
Strakos. Krylov subspace methods. 2013. See also
hegedus()
.  tol –
(optional) the tolerance for the stopping criterion with respect to the relative residual norm:
\[\frac{ \ M M_l (bA (x_0+M_r y_k))\_{M^{1}} } { \M M_l b\_{M^{1}}} \leq \text{tol}\]  maxiter – (optional) maximum number of iterations. Defaults to N.
 explicit_residual – (optional)
if set to
False
(default), the updated residual norm from the used method is used in each iteration. If set toTrue
, the residual is computed explicitly in each iteration and thus requires an additional application ofM
,Ml
,A
andMr
in each iteration.  store_arnoldi – (optional)
if set to
True
then the computed Arnoldi basis and the Hessenberg matrix are set as attributesV
andH
on the returned object. IfM
is notNone
, then alsoP
is set whereV=M*P
. Defaults toFalse
. If the method is based on the Lanczos method (e.g.,Cg
orMinres
), thenH
is real, symmetric and tridiagonal.  dtype – (optional) an optional dtype that is used to determine the dtype for the Arnoldi/Lanczos basis and matrix.
Upon convergence, the instance contains the following attributes:
xk
: the approximate solution \(x_k\).resnorms
: relative residual norms of all iterations, see parametertol
.errnorms
: the error norms of all iterations ifexact_solution
was provided.V
,H
andP
ifstore_arnoldi==True
, seestore_arnoldi
If the solver does not converge, a
ConvergenceError
is thrown which can be used to examine the misconvergence.
errnorms
= None¶ Error norms.

iter
= None¶ Iteration number.

static
operations
(nsteps)¶ Returns the number of operations needed for nsteps of the solver.
Parameters: nsteps – number of steps. Returns: a dictionary with the same keys as the timings parameter. Each value is the number of operations of the corresponding type for nsteps
iterations of the method.

resnorms
= None¶ Relative residual norms as described for parameter
tol
.

xk
= None¶ Approximate solution.
 linear_system – a