Linear solvers

HPD solve

Implementation

Solves \(AX=B\), \(A^T X = B\), or \(A^H X=B\) for \(X\) given Hermitian positive-definite (HPD) \(A\) and right-hand side matrix \(B\) (note that these options are all identical except for when \(A\). is complex). The solution is computed by first finding the Cholesky factorization of \(A\) and then performing two successive triangular solves against \(B\).

void HPDSolve(UpperOrLower uplo, Orientation orientation, Matrix<F> &A, Matrix<F> &B)

void HPDSolve(UpperOrLower uplo, Orientation orientation, DistMatrix<F> &A, DistMatrix<F> &B)

Overwrite B with the solution to \(AX=B\) or \(A^T X=B\), where A is Hermitian positive-definite and only the triangle of A specified by uplo is accessed.

Symmetric solve

Implementation

Solve \(AX=B\), \(A^T X = B\), or \(A^H X = B\) for \(X\) given a symmetric or Hermitian matrix \(A\) and a right-hand side matrix \(B\) using Bunch-Kaufman.

void SymmetricSolve(UpperOrLower uplo, Orientation orientation, Matrix<F> &A, Matrix<F> &B, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

void SymmetricSolve(UpperOrLower uplo, Orientation orientation, DistMatrix<F> &A, DistMatrix<F> &B, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

Overwrites \(B\) with the solution to the linear system. \(A\) is assumed symmetric if conjugate is false, and Hermitian otherwise.

Note

Only the lower-storage case is currently supported.

Hermitian solve

Implementation

Solve \(AX=B\), \(A^T X = B\), or \(A^H X = B\) for \(X\) given a Hermitian matrix \(A\) and a right-hand side matrix \(B\) using Bunch-Kaufman.

void HermitianSolve(UpperOrLower uplo, Orientation orientation, Matrix<F> &A, Matrix<F> &B, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

void HermitianSolve(UpperOrLower uplo, Orientation orientation, DistMatrix<F> &A, DistMatrix<F> &B, LDLPivotType pivotType = BUNCH_KAUFMAN_A)

Overwrites \(B\) with the solution to the linear system.

Note

Only the lower-storage case is currently supported, as this is a simple wrapper around SymmetricSolve().

Gaussian elimination

Implementation

Solves \(AX=B\) for \(X\) given a general square nonsingular matrix \(A\) and right-hand side matrix \(B\). The solution is computed through (partially pivoted) Gaussian elimination.

void GaussianElimination(Matrix<F> &A, Matrix<F> &B)

void GaussianElimination(DistMatrix<F> &A, DistMatrix<F> &B)

Upon completion, \(A\) will have been overwritten with Gaussian elimination and \(B\) will be overwritten with \(X\).

Least Squares

Implementation

Given \(A \in \mathbb{F}^{m \times n}\) and a right-hand side \(b \in \mathbb{F}^m\), a least-squares method solves \(A x \approx b\) differently depending upon whether \(m \ge n\).

When \(m \ge n\), there are at least as many constraints as degrees of freedom, and so a solution is sought for

\[\min_x \| A x - b \|_2\]

This problem is solved through the use of QR().

When \(m < n\), the problem is under-constrained and a solution is sought for the problem

\[\min_x \| x \|_2 \;\;\; \text{such that } A x = b.\]

This problem is solved through the use of LQ().

The above optimization problems can be readily generalized to multiple right-hand sides by switching to Frobenius norms.

void LeastSquares(Orientation orientation, Matrix<F> &A, const Matrix<F> &B, Matrix<F> &X)

void LeastSquares(Orientation orientation, DistMatrix<F> &A, const DistMatrix<F> &B, DistMatrix<F> &X)

If orientation is set to NORMAL, then solve \(AX=B\), otherwise orientation must be equal to ADJOINT and \(A^H X=B\) will be solved. Upon completion, \(A\) is overwritten with its QR or LQ factorization, and each column of \(X\) is overwritten with a solution vector.

General (Gauss-Markov) Linear Model (GLM)

Implementation

Example driver

\[\min_{X,Y} \| Y \|_F \;\;\; \text{subject to } A X + B Y = D.\]

void GLM(Matrix<F> &A, Matrix<F> &B, Matrix<F> &D, Matrix<F> &Y)

void GLM(DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<F> &D, DistMatrix<F> &Y)

Equality-constrained Least Squares (LSE)

Implementation

Example driver

\[\min_X \| A X - C \|_F \;\;\; \text{subject to } B X = D.\]

void LSE(Matrix<F> &A, Matrix<F> &B, Matrix<F> &C, Matrix<F> &D, Matrix<F> &X, bool computeResidual = false)

void LSE(DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<F> &C, DistMatrix<F> &D, DistMatrix<F> &X, bool computeResidual = false)

Multi-shift Hessenberg solves

Implementation

Solve for \(X\) in the system

\[H^\# X - X D^\# = Y\]

where \(H\) is Hessenberg, \(D\) is diagonal, and \(A^\#\) is defined to be one of \(\{A,A^T,A^H\}\).

void MultiShiftHessSolve(UpperOrLower uplo, Orientation orientation, F alpha, const Matrix<F> &H, const Matrix<F> &shifts, Matrix<F> &X)

void MultiShiftHessSolve(UpperOrLower uplo, Orientation orientation, F alpha, const DistMatrix<F, UH, VH> &H, const DistMatrix<F, VX, STAR> &shifts, DistMatrix<F, STAR, VX> &X)

Overwrite the columns of X with the solutions to shifted linear systems.

Note

Only a few subcases are currently supported, as this was added as part of HessenbergPseudospectrum()