Linear algebra

Psydac solvers

This module provides iterative solvers and preconditioners.

psydac.linalg.solvers.inverse(A, solver, **kwargs)

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:

• A (psydac.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).

• solver (str) – Preferred iterative solver. Options are: ‘cg’, ‘pcg’, ‘bicg’, ‘bicgstab’, ‘pbicgstab’, ‘minres’, ‘lsmr’, ‘gmres’.

Returns:

obj – A linear operator acting as the inverse of A, of the chosen subclass (for example psydac.linalg.solvers.ConjugateGradient).

Return type:

psydac.linalg.basic.InverseLinearOperator

class psydac.linalg.solvers.ConjugateGradient(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for L2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• verbose (bool) – If True, L2-norm of residual r is printed at each iteration.

• recycle (bool) – 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)

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:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

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

Return type:

psydac.linalg.basic.Vector

References

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

dot(b, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.PConjugateGradient(A, *, pc=None, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)

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:

• A (psydac.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.

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for L2-norm of residual r = A x - b. (Default: 1e-6)

• maxiter (int) – Maximum number of iterations. (Default: 1000)

• verbose (bool) – If True, the L2-norm of the residual r is printed at each iteration. (Default: False)

• recycle (bool) – 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)

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:

• b (psydac.linalg.stencil.StencilVector) – Right-hand-side vector of linear system.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

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

Return type:

psydac.linalg.basic.Vector

dot(b, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.BiConjugateGradient(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for 2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• verbose (bool) – If True, 2-norm of residual r is printed at each iteration.

• recycle (bool) – 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)

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:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

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

Return type:

psydac.linalg.basic.Vector

References

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

dot(b, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.BiConjugateGradientStabilized(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for 2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• verbose (bool) – If True, 2-norm of residual r is printed at each iteration.

• recycle (bool) – 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)

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

Parameters:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

• x (psydac.linalg.basic.Vector) – Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

• info (dict) –

  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)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.PBiConjugateGradientStabilized(A, *, pc=None, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for 2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• verbose (bool) – 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)

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

Parameters:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

• x (psydac.linalg.basic.Vector) – Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

• info (dict) –

  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)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.MinimumResidual(A, *, x0=None, tol=1e-06, maxiter=1000, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for 2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• verbose (bool) – If True, 2-norm of residual r is printed at each iteration.

• recycle (bool) – 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)

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:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

• x (psydac.linalg.basic.Vector) – Numerical solution of linear system. To check the convergence of the solver, use the method InverseLinearOperator.get_info().

• info (dict) – 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)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.LSMR(A, *, x0=None, tol=None, atol=None, btol=None, maxiter=1000, conlim=100000000.0, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for 2-norm of residual r = A*x - b.

• atol (float) – Absolute tolerance for 2-norm of residual r = A*x - b.

• btol (float) – Relative tolerance for 2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• conlim (float) – lsmr terminates if an estimate of cond(A) exceeds conlim.

• verbose (bool) – If True, 2-norm of residual r is printed at each iteration.

• recycle (bool) – 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()

solve(b, out=None)

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:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

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

Return type:

psydac.linalg.basic.Vector

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)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

class psydac.linalg.solvers.GMRES(A, *, x0=None, tol=1e-06, maxiter=100, verbose=False, recycle=False)

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:

• A (psydac.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).

• x0 (psydac.linalg.basic.Vector) – First guess of solution for iterative solver (optional).

• tol (float) – Absolute tolerance for L2-norm of residual r = A*x - b.

• maxiter (int) – Maximum number of iterations.

• verbose (bool) – If True, L2-norm of residual r is printed at each iteration.

• recycle (bool) – 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)

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:

• b (psydac.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 ‘inner(p)’ functions (b.inner(p) is the vector inner product b*p); moreover, scalar multiplication and sum operations are available.

• out (psydac.linalg.basic.Vector | NoneType) – The output vector, or None (optional).

Returns:

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

Return type:

psydac.linalg.basic.Vector

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)

arnoldi(k, p)

apply_givens_rotation(k, sn, cn)

dot(b, out=None)

Apply the LinearOperator self to the Vector v.

The result is written to the Vector out, if provided.

Parameters:

• v (Vector) – The vector to which the linear operator (self) is applied. It must belong to the domain of self.

• out (Vector) – The vector in which the result of the operation is stored. It must belong to the codomain of self. If out is None, a new vector is created and returned.

Returns:

The result of the operation. If out is None, a new vector is returned. Otherwise, the result is stored in out and out is returned.

Return type:

Vector

Block systems

class struphy.linear_algebra.schur_solver.SchurSolver(A: LinearOperator, BC: LinearOperator, solver_name: str, precond=None, solver_params: SolverParameters | None = None)

Solves for \(x^{n+1}\) in the block system

\[\begin{split}\left( \matrix{ A & \Delta t B \cr \Delta t C & \\text{Id} } \\right) \left( \matrix{ x^{n+1} \cr y^{n+1} } \\right) = \left( \matrix{ A & - \Delta t B \cr - \Delta t C & \\text{Id} } \\right) \left( \matrix{ x^n \cr y^n } \\right)\end{split}\]

using the Schur complement \(S = A - \Delta t^2 BC\), where Id is the identity matrix and \((x^n, y^n)\) is given. The solution is given by

\[\begin{split}x^{n+1} = S^{-1} \left[ (A + \Delta t^2 BC) \, x^n - 2 \Delta t B \, y^n \\right] \,.\end{split}\]

Parameters:

• A (LinearOperator) – Upper left block from [[A B], [C Id]].

• BC (LinearOperator) – Product from [[A B], [C Id]].

• solver_name (str) – See [psydac.linalg.solvers](pyccel/psydac) for possible names.

• **solver_params – Must correspond to the chosen solver.

property A

Upper left block from [[A B], [C Id]].

property BC

Product from [[A B], [C Id]].

class struphy.linear_algebra.schur_solver.SchurSolverFull(M, solver_name, **solver_params)

Solves the block system

\[\begin{split}\left( \matrix{ A & B \cr C & \\text{Id} } \\right) \left( \matrix{ x \cr y } \\right) = \left( \matrix{ b_x \cr b_y } \\right)\end{split}\]

using the Schur complement \(S = A - BC\), where Id is the identity matrix and \((b_x, b_y)^T\) is given. The solution is given by

\[ \begin{align}\begin{aligned}x &= S^{-1} \, (b_x - B b_y ) \,,\\y &= b_y - C x \,.\end{aligned}\end{align} \]

Parameters:

• M (BlockLinearOperator) – Matrix [[A B], [C Id]].

• solver_name (str) – See [psydac.linalg.solvers](pyccel/psydac) for possible names.

• **solver_params – Must correspond to the chosen solver.

dot(v, out=None)

Solves the 2x2 block matrix linear system.

Parameters:

• v (psydac.linalg.basic.Vector) – Left hand side of the system.

• out (psydac.linalg.basic.Vector, optional) – If given, the converged solution will be written into this vector (in-place).

Returns:

• out (psydac.linalg.block.BLockVector) – Converged solution.

• info (dict) – Convergence information.

class struphy.linear_algebra.schur_solver.SchurSolverFull3(M, solver_name, **solver_params)

Solves the block system

\[\begin{split}\left( \matrix{ A & B & D \cr C & \\text{Id} & 0 \cr E & 0 & \\text{Id} \cr } \\right) \left( \matrix{ x \cr y \cr z } \\right) = \left( \matrix{ b_x \cr b_y \cr b_z } \\right)\end{split}\]

using the Schur complement \(S = A - BC - DE\), where Id is the identity matrix and \((b_x, b_y, b_z)^T\) is given. The solution is given by

\[ \begin{align}\begin{aligned}x &= S^{-1} \, (b_x - B b_y - D b_z) \,,\\y &= b_y - C x \,,\\z &= b_z - E x \,.\end{aligned}\end{align} \]

Parameters:

• M (BlockLinearOperator) – Matrix [[A B D], [C Id 0], [E 0 Id]].

• solver_name (str) – See [psydac.linalg.solvers](pyccel/psydac) for possible names.

• **solver_params – Must correspond to the chosen solver.

dot(v, out=None)

Solves the 3x3 block matrix linear system.

Parameters:

• v (psydac.linalg.basic.Vector) – Left hand side of the system.

• out (psydac.linalg.basic.Vector, optional) – If given, the converged solution will be written into this vector (in-place).

Returns:

• out (psydac.linalg.block.BLockVector) – Converged solution.

• info (dict) – Convergence information.

class struphy.linear_algebra.saddle_point.SaddlePointSolver(A: list | LinearOperator | BlockLinearOperator, B: list | LinearOperator | BlockLinearOperator, F: list | Vector | BlockVector, Apre: list | None = None, method_to_solve: str = 'DirectNPInverse', preconditioner: bool = False, spectralanalysis: bool = False, dimension: str = '2D', solver_name: str = 'GMRES', tol: float = 1e-08, max_iter: int = 1000, **solver_params)

Solves for \((x, y)\) in the saddle point problem

\[\left( \matrix{ A & B^{\top} \cr B & 0 } \right) \left( \matrix{ x \cr y } \right) = \left( \matrix{ f \cr 0 } \right)\]

using either the Uzawa iteration \(BA^{-1}B^{\top} y = BA^{-1} f\) or using on of the solvers given in psydac.linalg.solvers. The prefered solver is GMRES. The decission which variant to use is given by the type of A. If A is of type list of xp.ndarrays or sc.sparse.csr_matrices, then this class uses the Uzawa algorithm. If A is of type LinearOperator or BlockLinearOperator, a solver is used for the inverse. Using the Uzawa algorithm, solution is given by:

\[y = \left[ B A^{-1} B^{\top}\right]^{-1} B A^{-1} f \,, \qquad x = A^{-1} \left[ f - B^{\top} y \right] \,.\]

Parameters:

• A (list, LinearOperator or BlockLinearOperator) – Upper left block. Either the entries on the diagonals of block A are given as list of xp.ndarray or sc.sparse.csr_matrix. Alternative: Give whole matrice A as LinearOperator or BlockLinearOperator. list: Uzawa algorithm is used. LinearOperator: A solver given in psydac.linalg.solvers is used. Specified by solver_name. BlockLinearOperator: A solver given in psydac.linalg.solvers is used. Specified by solver_name.

• B (list, LinearOperator or BlockLinearOperator) – Lower left block. Uzwaw Algorithm: All entries of block B are given either as list of xp.ndarray or sc.sparse.csr_matrix. Solver: Give whole B as LinearOperator or BlocklinearOperator

• F (list) – Right hand side of the upper block. Uzawa: Given as list of xp.ndarray or sc.sparse.csr_matrix. Solver: Given as LinearOperator or BlockLinearOperator

• Apre (list) – The non-inverted preconditioner for entries on the diagonals of block A are given as list of xp.ndarray or sc.sparse.csr_matrix. Only required for the Uzawa algorithm.

• method_to_solve (str) – Method for the inverses. Choose from ‘DirectNPInverse’, ‘ScipySparse’, ‘InexactNPInverse’ ,’SparseSolver’. Only required for the Uzawa algorithm.

• preconditioner (bool) – Wheter to use preconditioners given in Apre or not. Only required for the Uzawa algorithm.

• spectralanalysis (bool) – Do the spectralanalyis for the matrices in A and if preconditioner given, compare them to the preconditioned matrices. Only possible if A is given as list.

• dimension (str) – Which of the predefined manufactured solutions to use (‘1D’, ‘2D’ or ‘Restelli’)

• tol (float) – Convergence tolerance for the potential residual.

• max_iter (int) – Maximum number of iterations allowed.

property A

Upper left block.

property B

Lower left block.

property F

Right hand side vector.

property Apre

Preconditioner for upper left block A.

Products

struphy.linear_algebra.linalg_kernels.matrix_vector(a: float[:, :], b: float[:], c: float[:])

Performs the matrix-vector product of a 3x3 matrix with a vector.

Parameters:

• a (array[float]) – The input array (matrix) of shape (3,3).

• b (array[float]) – The input array (vector) of shape (3,).

• c (array[float]) – The output array (vector) of shape (3,) which is the result of the matrix-vector product a.dot(b).

struphy.linear_algebra.linalg_kernels.matrix_matrix(a: float[:, :], b: float[:, :], c: float[:, :])

Performs the matrix-matrix product of a 3x3 matrix with another 3x3 matrix.

Parameters:

• a (array[float]) – The first input array (matrix) of shape (3,3).

• b (array[float]) – The second input array (matrix) of shape (3,3).

• c (array[float]) – The output array (matrix) of shape (3,3) which is the result of the matrix-matrix product a.dot(b).

struphy.linear_algebra.linalg_kernels.transpose(a: float[:, :], b: float[:, :])

Assembles the transposed of a 3x3 matrix.

Parameters:

• a (array[float]) – The input array (matrix) of shape (3,3).

• b (array[float]) – The output array (matrix) of shape (3,3).

struphy.linear_algebra.linalg_kernels.scalar_dot(a: Final[float[:]], b: Final[float[:]]) → float

Computes scalar (dot) product of two vectors of length 3.

Parameters:

• a (array[float]) – The first input array (vector) of shape (3,).

• b (array[float]) – The second input array (vector) of shape (3,).

Returns:

value – The scalar poduct of the two input vectors a and b.

Return type:

float

struphy.linear_algebra.linalg_kernels.det(a: Final[float[:, :]]) → float

Computes the determinant of a 3x3 matrix.

Parameters:

a (array[float]) – The input array (matrix) of shape (3,3) of which the determinant shall be computed.

Returns:

det_a – The determinant of the 3x3 matrix a.

Return type:

float

struphy.linear_algebra.linalg_kernels.cross(a: float[:], b: float[:], c: float[:])

Computes the vector (cross) product of two vectors of length 3.

Parameters:

• a (array[float]) – The first input array (vector) of shape (3,).

• b (array[float]) – The second input array (vector) of shape (3,).

• c (array[float]) – The output array (vector) of shape (3,) which is the vector product a x b.

struphy.linear_algebra.linalg_kernels.outer(a: float[:], b: float[:], c: float[:, :])

Computes the outer product of two vectors of length 3.

Parameters:

• a (array[float]) – The first input array (vector) of shape (3,).

• b (array[float]) – The second input array (vector) of shape (3,).

• c (array[float]) – The output array (matrix) of shape (3, 3) which is the outer product c_ij = a_i*b_j.

struphy.linear_algebra.linalg_kernels.matrix_inv(a: float[:, :], b: float[:, :])

Computes the inverse of a 3x3 matrix.

Parameters:

• a (array[float]) – The input array (matrix) of shape (3,3).

• b (array[float]) – The output array (matrix) of shape (3,3).

struphy.linear_algebra.linalg_kernels.matrix_inv_with_det(a: float[:, :], det_a: float, b: float[:, :])

Computes the inverse of a 3x3 matrix for the case that the determinant is already known such that its extra compuation can be avoided.

Parameters:

• a (array[float]) – The input array (matrix) of shape (3,3).

• det_a (float) – The determinant of the input array (matrix) a.

• b (array[float]) – The output array (matrix) of shape (3,3).

struphy.linear_algebra.linalg_kernels.matrix_vector4(a: float[:, :], b: float[:], c: float[:])

Performs the matrix-vector product of a 4x4 matrix with a vector.

Parameters:

• a (array[float]) – The input array (matrix) of shape (4,4).

• b (array[float]) – The input array (vector) of shape (4,).

• c (array[float]) – The output array (vector) of shape (4,) which is the result of the matrix-vector product a.dot(b).

struphy.linear_algebra.linalg_kernels.matrix_matrix4(a: float[:, :], b: float[:, :], c: float[:, :])

Performs the matrix-matrix product of a 4x4 matrix with another 4x4 matrix.

Parameters:

• a (array[float]) – The first input array (matrix) of shape (4,4).

• b (array[float]) – The second input array (matrix) of shape (4,4).

• c (array[float]) – The output array (matrix) of shape (4,4) which is the result of the matrix-matrix product a.dot(b).

struphy.linear_algebra.linalg_kernels.det4(a: float[:, :]) → float

Computes the determinant of a 4x4 matrix.

Parameters:

a (array[float]) – The input array (matrix) of shape (4,4) of which the determinant shall be computed.

Returns:

det_a – The determinant of the 3x3 matrix a.

Return type:

float

Matrix-vector product and solvers for matrices with a 3D Kronecker product structure.

M_ijkmno = A_im * B_jn * C_ko

where matrices A, B, C stem from 1D problems.

COMMENT: the reshape of a matrix can be viewed as ravel+reshape.

Let r = (r_ijk) be a 3D matrix of size M*N*O. ravel(r) = [r_111, r112, … , r_MNO] (row major always –> last index runs fastest) reshape(ravel(r), (M, N*O)) = [[r_111, r112, … , r_1NO],

[r_211, r212, … , r_2NO], …, [r_M11, rM12, … , r_MNO]]

struphy.linear_algebra.linalg_kron.kron_matvec_2d(kmat, vec2d)

2D Kronecker matrix-vector product.

res_ij = (A_ik * B_jl) * vec2d_kl

implemented as two matrix-matrix multiplications with intermediate transpose:

step1(v1,k0) = ( kmat[0](k0,v0) * vec2d(v0,v1) )^T step2(k0,k1) = ( kmat[1](k1,v1) * step1(v1,k0) )^T res = step2(k0,k1)

Parameters:

• kmat (2 sparse matrices for each direction of size (k0,v0) and (k1,v1))

• vec2d (2d array of size (v0,v1))

Returns:

res

Return type:

2d array of size (k0,k1)

struphy.linear_algebra.linalg_kron.kron_matvec_3d(kmat, vec3d)

3D Kronecker matrix-vector product.

res_ijk = (A_im * B_jn * C_ko) * vec3d_mno

implemented as three matrix-matrix multiplications with intermediate reshape and transpose. step1(v1*v2,k0) <= ( kmat[0](k0,v0) * reshaped_vec3d(v0,v1*v2) )^T step2(v2*k0,k1) <= ( kmat[1](k1,v1) * reshaped_step1(v1,v2*k0) )^T step3(k0*k1*k2) <= ( kmat[2](k2,v2) * reshaped_step2(v2,k0*k1) )^T res <= reshaped_step3(k0,k1,k2)

no overhead of numpy reshape command, as they do NOT copy the data.

Parameters:

• kmat (3 sparse matrices for each direction,) – of size (k0,v0),(k1,v1),(k2,v2)

• vec3d (3d array of size (v0,v1,v2))

Returns:

res

Return type:

3d array of size (k0,k1,k2)

struphy.linear_algebra.linalg_kron.kron_matvec_3d_1(kmat, vec3d)

struphy.linear_algebra.linalg_kron.kron_matvec_3d_2(kmat, vec3d)

struphy.linear_algebra.linalg_kron.kron_matvec_3d_3(kmat, vec3d)

struphy.linear_algebra.linalg_kron.kron_matvec_3d_23(kmat, vec3d)

struphy.linear_algebra.linalg_kron.kron_matvec_3d_13(kmat, vec3d)

struphy.linear_algebra.linalg_kron.kron_matvec_3d_12(kmat, vec3d)

struphy.linear_algebra.linalg_kron.kron_matmat_fft_3d(a_vec, b_vec)

matrix-matrix product between 3d kronecker circulant matrices.

res_ijk,qrs = (A_im * B_jn * C_ko) * (D_mq * E_nr * F_os)

the 1d matrix-matrix product is computed as A*D = ifft(fft(a)*fft(d))

This relies on the fact that the eigenvalues of a circulant matrix are given by

lambda = fft(c) = ifft(c)*N for c real

Parameters:

• a_vec (3 vectors for each direction defining the circulant matrices of the first matrix)

• b_vec (3 vectors for each direction defining the circulant matrices of the second matrix)

Returns:

c_vec

Return type:

3 vectors for each direction defining the resulting circulant matrices

struphy.linear_algebra.linalg_kron.kron_lusolve_3d(kmatlu, rhs)

3D Kronecker LU solver.

solve for x: (A_im * B_jn * C_ko) * x_mno = rhs_ijk

implemented as three matrix-matrix solve with intermediate reshape and transpose. step1(r1*r2,r0) <= ( A(r0,r0)^-1 * reshaped_rhs(r0,r1*r2) )^T step2(r2*r0,r1) <= ( B(r1,r1)^-1 * reshaped_step1(r1,r2*r0) )^T step3(r0*r1*r2) <= ( C(r2,r2)^-1 * reshaped_step2(r2,r0*r1) )^T res <= reshaped_step3(r0,r1,r2)

no overhead of numpy reshape command, as they do NOT copy the data.

Parameters:

• kmatlu (3 already LU decompositions of sparse matrices for each direction,) – of size (r0,r0),(r1,r1),(r2,r2)

• rhs (3d array of size (r0,r1,r2), right-hand size)

Returns:

res

Return type:

3d array of size (r0,r1,r2), solution

struphy.linear_algebra.linalg_kron.kron_solve_3d(kmat, rhs)

3D Kronecker solver.

solve for x: (A_im * B_jn * C_ko) * x_mno = rhs_ijk

implemented as three matrix-matrix solve with intermediate reshape and transpose. step1(r1*r2,r0) <= ( A(r0,r0)^-1 * reshaped_rhs(r0,r1*r2) )^T step2(r2*r0,r1) <= ( B(r1,r1)^-1 * reshaped_step1(r1,r2*r0) )^T step3(r0*r1*r2) <= ( C(r2,r2)^-1 * reshaped_step2(r2,r0*r1) )^T res <= reshaped_step3(r0,r1,r2)

no overhead of numpy reshape command, as they do NOT copy the data.

Parameters:

• kmat (3 sparse matrices for each direction,) – of size (r0,r0),(r1,r1),(r2,r2)

• rhs (3d array of size (r0,r1,r2), right-hand size)

Returns:

res

Return type:

3d array of size (r0,r1,r2), solution

struphy.linear_algebra.linalg_kron.kron_fftsolve_3d(cvec, rhs)

3D Kronecker fft solver for circulant matrices.

solve for x: (A_im * B_jn * C_ko) * x_mno = rhs_ijk

implemented as three matrix-matrix solve with intermediate reshape and transpose. step1(r1*r2,r0) <= ( A(r0,r0)^-1 * reshaped_rhs(r0,r1*r2) )^T step2(r2*r0,r1) <= ( B(r1,r1)^-1 * reshaped_step1(r1,r2*r0) )^T step3(r0*r1*r2) <= ( C(r2,r2)^-1 * reshaped_step2(r2,r0*r1) )^T res <= reshaped_step3(r0,r1,r2)

no overhead of numpy reshape command, as they do NOT copy the data.

Parameters

cvec : 3 vectors of size (r0),(r1),(r2) defining 3 circulant matrices for each direction,

rhs : 3d array of size (r0,r1,r2), right-hand size

res : 3d array of size (r0,r1,r2), solution

struphy.linear_algebra.linalg_kron.kron_lusolve_2d(kmatlu, rhs)

2D Kronecker LU solver.

solve for x: (A_im * B_jn) * x_mn = rhs_ij

Parameters:

• kmatlu (2 already LU decompositions of sparse matrices for each direction,) – of size (r0,r0),(r1,r1)

• rhs (2d array of size (r0,r1), right-hand size)

Returns:

res

Return type:

2d array of size (r0,r1), solution

struphy.linear_algebra.linalg_kron.kron_fftsolve_2d(A_LU, b_vec, rhs)

ALERT: docstring wrong.

Solve for 3d vector, matrix would be a 3d kronecker circulant matrix,

but system is only solved in each direction.

solve for x: (A_im * B_jn * C_ko) * x_mno = rhs_ijk

implemented as three matrix-matrix solve with intermediate reshape and transpose. step1(r1*r2,r0) <= ( A(r0,r0)^-1 * reshaped_rhs(r0,r1*r2) )^T step2(r2*r0,r1) <= ( B(r1,r1)^-1 * reshaped_step1(r1,r2*r0) )^T step3(r0*r1*r2) <= ( C(r2,r2)^-1 * reshaped_step2(r2,r0*r1) )^T res <= reshaped_step3(r0,r1,r2)

no overhead of numpy reshape command, as they do NOT copy the data.

COMMENT: the reshape of a matrix can be viewed as ravel+reshape. Let r = (r_ijk) be a 3D matrix of size M*N*O. ravel(r) = [r_111, r112, … , r_MNO] (row major always –> last index runs fastest) reshape(ravel(r), (M, N*O)) = [[r_111, r112, … , r_1NO],

[r_211, r212, … , r_2NO], …, [r_M11, rM12, … , r_MNO]]

Parameters

cvec : 3 vectors of size (r0),(r1),(r2) defining 3 circulant matrices for each direction,

rhs : 3d array of size (r0,r1,r2), right-hand size

res : 3d array of size (r0,r1,r2), solution