linalg.solvers#

This module provides iterative solvers and preconditioners.

Functions#

inverse(A, solver, **kwargs)

A function to create objects of all InverseLinearOperator subclasses.

Classes#

Inheritance diagram of psydac.linalg.solvers

`BiConjugateGradient`(A, *[, x0, tol, ...])	Biconjugate Gradient (BiCG).
`BiConjugateGradientStabilized`(A, *[, x0, ...])	Biconjugate Gradient Stabilized (BiCGStab).
`ConjugateGradient`(A, *[, x0, tol, maxiter, ...])	Conjugate Gradient (CG).
`GMRES`(A, *[, x0, tol, maxiter, verbose, recycle])	Generalized Minimal Residual (GMRES).
`LSMR`(A, *[, x0, tol, atol, btol, maxiter, ...])	Least Squares Minimal Residual (LSMR).
`MinimumResidual`(A, *[, x0, tol, maxiter, ...])	Minimum Residual (MinRes).
`PBiConjugateGradientStabilized`(A, *[, pc, ...])	Preconditioned Biconjugate Gradient Stabilized (PBiCGStab).
`PConjugateGradient`(A, *[, pc, x0, tol, ...])	Preconditioned Conjugate Gradient (PCG).

Details#

This module provides iterative solvers and preconditioners.

inverse(A, solver, **kwargs)[source]#

A function to create objects of all InverseLinearOperator subclasses.

These are, as of June 06, 2023: ConjugateGradient, PConjugateGradient, BiConjugateGradient, BiConjugateGradientStabilized, MinimumResidual, LSMR, GMRES.

The kwargs given must be compatible with the chosen solver subclass.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (e.g. a matrix-vector product A*p).
solverstr: Preferred iterative solver. Options are: ‘cg’, ‘pcg’, ‘bicg’, ‘bicgstab’, ‘pbicgstab’, ‘minres’, ‘lsmr’, ‘gmres’.

Returns:

objpsydac.linalg.basic.InverseLinearOperator: A linear operator acting as the inverse of A, of the chosen subclass (for example psydac.linalg.solvers.ConjugateGradient).

class ConjugateGradient(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Conjugate Gradient (CG).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on the Conjugate gradient algorithm for solving linear system Ax=b. Implementation from [1], page 137.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for L2-norm of residual r = A*x - b.
maxiterint: Maximum number of iterations.
verbosebool: If True, L2-norm of residual r is printed at each iteration.
recyclebool: Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

solve(b, out=None)[source]#

Conjugate gradient algorithm for solving linear system Ax=b. Only working if A is an hermitian and positive-definite linear operator. Implementation from [1], page 137. Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`.

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system Ax = b. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector: Numerical solution of the linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class PConjugateGradient(A, *, pc=None, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Preconditioned Conjugate Gradient (PCG).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on a preconditioned conjugate gradient method. The Preconditioned Conjugate Gradient (PCG) algorithm solves the linear system A x = b where A is a symmetric and positive-definite matrix, i.e. A = A^T and y A y > 0 for any vector y. The preconditioner P is a matrix which approximates the inverse of A. The algorithm assumes that P is also symmetric and positive definite.

Since this is a matrix-free iterative method, both A and P are provided as LinearOperator objects which must implement the dot method.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of the linear system. This should be symmetric and positive definite.
pc: psydac.linalg.basic.LinearOperator: Preconditioner which should approximate the inverse of A (optional). Like A, the preconditioner should be symmetric and positive definite.
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for L2-norm of residual r = A x - b. (Default: 1e-6)
maxiter: int: Maximum number of iterations. (Default: 1000)
verbosebool: If True, the L2-norm of the residual r is printed at each iteration. (Default: False)
recyclebool: If True, a copy of the output is stored in x0 to speed up consecutive calculations of slightly altered linear systems. (Default: False)

solve(b, out=None)[source]#

Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte system Ax = b. It assumes that pc.dot(r) returns the solution to Ps = r, where P is positive definite. Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`.

Parameters:

bpsydac.linalg.stencil.StencilVector: Right-hand-side vector of linear system.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector: Numerical solution of the linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class BiConjugateGradient(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Biconjugate Gradient (BiCG).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on the Biconjugate gradient (BCG) algorithm for solving linear system Ax=b. Implementation from [1], page 175.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for 2-norm of residual r = A*x - b.
maxiter: int: Maximum number of iterations.
verbosebool: If True, 2-norm of residual r is printed at each iteration.
recyclebool: Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

solve(b, out=None)[source]#

Biconjugate gradient (BCG) algorithm for solving linear system Ax=b. Implementation from [1], page 175. Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`. ToDo: Add optional preconditioner

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector: Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class BiConjugateGradientStabilized(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Biconjugate Gradient Stabilized (BiCGStab).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on the Biconjugate gradient Stabilized (BCGSTAB) algorithm for solving linear system Ax=b. Implementation from [1], page 175.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for 2-norm of residual r = A*x - b.
maxiter: int: Maximum number of iterations.
verbosebool: If True, 2-norm of residual r is printed at each iteration.
recyclebool: Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

solve(b, out=None)[source]#

Biconjugate gradient stabilized method (BCGSTAB) algorithm for solving linear system Ax=b. Implementation from [1], page 175. ToDo: Add optional preconditioner

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector

Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

infodict

Dictionary containing convergence information:

‘niter’ = (int) number of iterations
‘success’ = (boolean) whether convergence criteria have been met
‘res_norm’ = (float) 2-norm of residual vector r = A*x - b.

References

[1] H. A. van der Vorst. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comp., 13(2):631–644, 1992.

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class PBiConjugateGradientStabilized(A, *, pc=None, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Preconditioned Biconjugate Gradient Stabilized (PBiCGStab).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on the preconditioned Biconjugate gradient Stabilized (PBCGSTAB) algorithm for solving linear system Ax=b. Implementation from [1], page 251.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
pc: psydac.linalg.basic.LinearOperator: Preconditioner for A, it should approximate the inverse of A (can be None).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for 2-norm of residual r = A*x - b.
maxiter: int: Maximum number of iterations.
verbosebool: If True, 2-norm of residual r is printed at each iteration.

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

solve(b, out=None)[source]#

Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b. Implementation from [1], page 251.

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector

Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

infodict

Dictionary containing convergence information:

‘niter’ = (int) number of iterations
‘success’ = (boolean) whether convergence criteria have been met
‘res_norm’ = (float) 2-norm of residual vector r = A*x - b.

References

[1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015.

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class MinimumResidual(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Minimum Residual (MinRes).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function Use MINimum RESidual iteration to solve Ax=b

MINRES minimizes norm(A*x - b) for a real symmetric matrix A. Unlike the Conjugate Gradient method, A can be indefinite or singular.

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for 2-norm of residual r = A*x - b.
maxiter: int: Maximum number of iterations.
verbosebool: If True, 2-norm of residual r is printed at each iteration.
recyclebool: Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems

Notes

This is an adaptation of the MINRES Solver in Scipy, where the method is modified to accept Psydac data structures, scipy/scipy

References

Solution of sparse indefinite systems of linear equations, C. C. Paige and M. A. Saunders (1975), SIAM J. Numer. Anal. 12(4), pp. 617-629. https://web.stanford.edu/group/SOL/software/minres/

solve(b, out=None)[source]#

Use MINimum RESidual iteration to solve Ax=b MINRES minimizes norm(A*x - b) for a real symmetric matrix A. Unlike the Conjugate Gradient method, A can be indefinite or singular. Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`.

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector: Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().
infodict: Dictionary containing convergence information: - ‘niter’ = (int) number of iterations - ‘success’ = (boolean) whether convergence criteria have been met - ‘res_norm’ = (float) 2-norm of residual vector r = A*x - b.

Notes

This is an adaptation of the MINRES Solver in Scipy, where the method is modified to accept Psydac data structures, scipy/scipy References ———- Solution of sparse indefinite systems of linear equations, C. C. Paige and M. A. Saunders (1975), SIAM J. Numer. Anal. 12(4), pp. 617-629. https://web.stanford.edu/group/SOL/software/minres/

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class LSMR(A, *, x0=None, tol=None, atol=None, btol=None, maxiter=1000, conlim=100000000.0, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Least Squares Minimal Residual (LSMR).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on the Iterative solver for least-squares problems. lsmr solves the system of linear equations Ax = b. If the system is inconsistent, it solves the least-squares problem min ||b - Ax||_2. A is a rectangular matrix of dimension m-by-n, where all cases are allowed: m = n, m > n, or m < n. b is a vector of length m. The matrix A may be dense or sparse (usually sparse).

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for 2-norm of residual r = A*x - b.
atolfloat: Absolute tolerance for 2-norm of residual r = A*x - b.
btolfloat: Relative tolerance for 2-norm of residual r = A*x - b.
maxiter: int: Maximum number of iterations.
conlimfloat: lsmr terminates if an estimate of cond(A) exceeds conlim.
verbosebool: If True, 2-norm of residual r is printed at each iteration.
recyclebool: Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems

Notes

This is an adaptation of the LSMR Solver in Scipy, where the method is modified to accept Psydac data structures, scipy/scipy

References

[1]

D. C.-L. Fong and M. A. Saunders, “LSMR: An iterative algorithm for sparse least-squares problems”, SIAM J. Sci. Comput., vol. 33, pp. 2950-2971, 2011. arxiv:1006.0758

[2]

LSMR Software, https://web.stanford.edu/group/SOL/software/lsmr/

get_success()[source]#

solve(b, out=None)[source]#

Iterative solver for least-squares problems. lsmr solves the system of linear equations Ax = b. If the system is inconsistent, it solves the least-squares problem min ||b - Ax||_2. A is a rectangular matrix of dimension m-by-n, where all cases are allowed: m = n, m > n, or m < n. b is a vector of length m. The matrix A may be dense or sparse (usually sparse). Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`.

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector: Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

Notes

This is an adaptation of the LSMR Solver in Scipy, where the method is modified to accept Psydac data structures, scipy/scipy

References

[1]

D. C.-L. Fong and M. A. Saunders, “LSMR: An iterative algorithm for sparse least-squares problems”, SIAM J. Sci. Comput., vol. 33, pp. 2950-2971, 2011. arxiv:1006.0758

[2]

LSMR Software, https://web.stanford.edu/group/SOL/software/lsmr/

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.

class GMRES(A, *, x0=None, tol=1e-06, maxiter=100, verbose=False, recycle=False)[source]#

Bases: InverseLinearOperator

Generalized Minimal Residual (GMRES).

A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on the generalized minimal residual algorithm for solving linear system Ax=b. Implementation from Wikipedia

Parameters:

Apsydac.linalg.basic.LinearOperator: Left-hand-side matrix A of linear system; individual entries A[i,j] can’t be accessed, but A has ‘shape’ attribute and provides ‘dot(p)’ function (i.e. matrix-vector product A*p).
x0psydac.linalg.basic.Vector: First guess of solution for iterative solver (optional).
tolfloat: Absolute tolerance for L2-norm of residual r = A*x - b.
maxiter: int: Maximum number of iterations.
verbosebool: If True, L2-norm of residual r is printed at each iteration.
recyclebool: Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems

References

[1] Y. Saad and M.H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems”, SIAM J. Sci. Stat. Comput., 7:856–869, 1986.

solve(b, out=None)[source]#

Generalized minimal residual algorithm for solving linear system Ax=b. Implementation from Wikipedia. Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`.

Parameters:

bpsydac.linalg.basic.Vector: Right-hand-side vector of linear system Ax = b. Individual entries b[i] need not be accessed, but b has ‘shape’ attribute and provides ‘copy()’ and ‘dot(p)’ functions (dot(p) is the vector inner product b*p ); moreover, scalar multiplication and sum operations are available.
outpsydac.linalg.basic.Vector | NoneType: The output vector, or None (optional).

Returns:

xpsydac.linalg.basic.Vector: Numerical solution of the linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

References

[1] Y. Saad and M.H. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems”, SIAM J. Sci. Stat. Comput., 7:856–869, 1986.

solve_triangular(T, d)[source]#

arnoldi(k, p)[source]#

apply_givens_rotation(k, sn, cn)[source]#

dot(b, out=None)[source]#: Apply linear operator to Vector v. Result is written to Vector out, if provided.