Matrix functions¶
Direct inversion¶
Triangular¶
Inverts a (possibly unitdiagonal) triangular matrix inplace.

void
TriangularInverse
(UpperOrLower uplo, UnitOrNonUnit diag, Matrix<F> &A)¶

void
TriangularInverse
(UpperOrLower uplo, UnitOrNonUnit diag, DistMatrix<F> &A)¶ Inverts the triangle of A specified by the parameter uplo; if diag is set to UNIT, then A is treated as unitdiagonal.
General¶
This routine computes the inplace inverse of a general fullypopulated (invertible) matrix \(A\) as
where \(PA=LU\) is the result of LU factorization with partial pivoting. The algorithm essentially factors \(A\), inverts \(U\) in place, solves against \(L\) one block column at a time, and then applies the row pivots in reverse order to the columns of the result.

void
Inverse
(Matrix<F> &A)¶

void
Inverse
(DistMatrix<F> &A)¶ Overwrites the general matrix A with its inverse.
Symmetric/Hermitian¶

void
SymmetricInverse
(UpperOrLower uplo, Matrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)¶

void
SymmetricInverse
(UpperOrLower uplo, DistMatrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)¶ Invert a symmetric or Hermitian matrix using a pivoted LDL factorization.

void
HermitianInverse
(UpperOrLower uplo, Matrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)¶

void
HermitianInverse
(UpperOrLower uplo, DistMatrix<F> &A, bool conjugate = false, LDLPivotType pivotType = BUNCH_KAUFMAN_A)¶ Invert a Hermitian matrix using a pivoted LDL factorization.
HPD¶
This routine uses a custom algorithm for computing the inverse of a Hermitian positivedefinite matrix \(A\) as
where \(L\) is the lower Cholesky factor of \(A\) (the upper Cholesky factor is computed in the case of uppertriangular storage). Rather than performing Cholesky factorization, triangular inversion, and then the Hermitian triangular outer product in sequence, this routine makes use of the singlesweep algorithm described in Bientinesi et al.’s “Families of algorithms related to the inversion of a symmetric positive definite matrix”, in particular, the variant 2 algorithm from Fig. 9.
If the matrix is found to not be HPD, then a NonHPDMatrixException
is thrown.

void
HPDInverse
(UpperOrLower uplo, Matrix<F> &A)¶

void
HPDInverse
(UpperOrLower uplo, DistMatrix<F> &A)¶ Overwrite the uplo triangle of the HPD matrix A with the same triangle of the inverse of A.
Hermitian functions¶
Reform the matrix with the eigenvalues modified by a userdefined function. When the userdefined function is realvalued, the result will remain Hermitian, but when the function is complexvalued, the result is best characterized as normal.
When the userdefined function, say \(f\), is analytic, we can say much more about the result: if the eigenvalue decomposition of the Hermitian matrix \(A\) is \(A=Z \Lambda Z^H\), then
Two important special cases are \(f(\lambda) = \exp(\lambda)\) and \(f(\lambda)=\exp(i \lambda)\), where the former results in a Hermitian matrix and the latter in a normal (in fact, unitary) matrix.
Note
Since Elemental currently depends on PMRRR for its tridiagonal eigensolver, only doubleprecision results are supported as of now.

void
RealHermitianFunction
(UpperOrLower uplo, Matrix<F> &A, const RealFunctor &f)¶

void
RealHermitianFunction
(UpperOrLower uplo, DistMatrix<F> &A, const RealFunctor &f)¶ Modifies the eigenvalues of the passedin Hermitian matrix by replacing each eigenvalue \(\lambda_i\) with \(f(\lambda_i) \in \mathbb{R}\).
RealFunctor
is any class which has the member functionReal operator()( Real omega ) const
.

void
ComplexHermitianFunction
(UpperOrLower uplo, Matrix<Complex<Real>> &A, const ComplexFunctor &f)¶

void
ComplexHermitianFunction
(UpperOrLower uplo, DistMatrix<Complex<Real>> &A, const ComplexFunctor &f)¶ Modifies the eigenvalues of the passedin complex Hermitian matrix by replacing each eigenvalue \(\lambda_i\) with \(f(\lambda_i) \in \mathbb{C}\).
ComplexFunctor
can be any class which has the member functionComplex<Real> operator()( Real omega ) const
.
TODO: A version of ComplexHermitianFunction which begins with a real matrix
Pseudoinverse¶

Pseudoinverse
(Matrix<F> &A, Base<F> tolerance = 0)¶

Pseudoinverse
(DistMatrix<F> &A, Base<F> tolerance = 0)¶ Computes the pseudoinverse of a general matrix through computing its SVD, modifying the singular values with the function
\[\begin{split}f(\sigma) = \left\{\begin{array}{cc} 1/\sigma, & \sigma \ge \epsilon \, n \, \ A \_2 \\ 0, & \mbox{otherwise} \end{array}\right.,\end{split}\]where \(\epsilon\) is the relative machine precision, \(n\) is the height of \(A\), and \(\ A \_2\) is the maximum singular value. If a nonzero value for tolerance was specified, it is used instead of \(\epsilon n \ A \_2\).

HermitianPseudoinverse
(UpperOrLower uplo, Matrix<F> &A, Base<F> tolerance = 0)¶

HermitianPseudoinverse
(UpperOrLower uplo, DistMatrix<F> &A, Base<F> tolerance = 0)¶ Computes the pseudoinverse of a Hermitian matrix through a customized version of
RealHermitianFunction()
which used the eigenvalue mapping function\[\begin{split}f(\omega) = \left\{\begin{array}{cc} 1/\omega, & \omega \ge \epsilon \, n \, \ A \_2 \\ 0, & \mbox{otherwise} \end{array}\right.,\end{split}\]where \(\epsilon\) is the relative machine precision, \(n\) is the height of \(A\), and \(\ A \_2\) can be computed as the maximum absolute value of the eigenvalues of \(A\). If a nonzero value for tolerance is specified, it is used instead of \(\epsilon n \ A \_2\).
Square root¶
A matrix \(B\) satisfying
is referred to as the squareroot of the matrix \(A\). Such a matrix is guaranteed to exist as long as \(A\) is diagonalizable: if \(A = X \Lambda X^{1}\), then we may put
where each eigenvalue \(\lambda = r e^{i\theta}\) maps to \(\sqrt{\lambda} = \sqrt{r} e^{i\theta/2}\).

void
SquareRoot
(Matrix<F> &A)¶

void
SquareRoot
(DistMatrix<F> &A)¶ Currently uses a Newton iteration to compute the general matrix squareroot. See
square_root::Newton
for the more detailed interface.

void
HPSDSquareRoot
(UpperOrLower uplo, Matrix<F> &A)¶

void
HPSDSquareRoot
(UpperOrLower uplo, DistMatrix<F> &A)¶ Computes the Hermitian EVD, squareroots the eigenvalues, and then reforms the matrix. If any of the eigenvalues were sufficiently negative, a
NonHPSDMatrixException
is thrown.
TODO: HermitianSquareRoot
square_root namespace¶

int
square_root::
Newton
(Matrix<F> &A, int maxIts = 100, Base<F> tol = 0)¶

int
square_root::
Newton
(DistMatrix<F> &A, int maxIts = 100, Base<F> tol = 0)¶ Performs at most
maxIts
Newton steps in an attempt to compute the matrix squareroot within the specified tolerance, which defaults to \(n \epsilon\), where \(n\) is the matrix height and \(\epsilon\) is the machine precision.
Sign¶
The matrix sign function can be written as
as long as \(A\) does not have any pureimaginary eigenvalues.

void
Sign
(Matrix<F> &A)¶

void
Sign
(DistMatrix<F> &A)¶

void
Sign
(Matrix<F> &A, Matrix<F> &N)¶

void
Sign
(DistMatrix<F> &A, DistMatrix<F> &N)¶ Compute the matrix sign through a globallyconvergent Newton iteration scaled with the Frobenius norm of the iterate and its inverse. Optionally return the full decomposition, \(A=S N\), where \(A\) is overwritten by \(S\).

void
HermitianSign
(UpperOrLower uplo, Matrix<F> &A)¶

void
HermitianSign
(UpperOrLower uplo, DistMatrix<F> &A)¶

void
HermitianSign
(UpperOrLower uplo, Matrix<F> &A, Matrix<F> &N)¶

void
HermitianSign
(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &N)¶ Compute the Hermitian EVD, replace the eigenvalues with their sign, and then reform the matrix. Optionally return the full decomposition, \(A=SN\), where \(A\) is overwritten by \(S\). Note that this will also be a polar decomposition.
sign namespace¶

type
sign::
Scaling
¶ An enum for specifying the scaling strategy to be used for the Newton iteration for the matrix sign function. It must be either
NONE
,DETERMINANT
, orFROB_NORM
(the default).

int
sign::
Newton
(Matrix<F> &A, sign::Scaling scaling = FROB_NORM, int maxIts = 100, Base<F> tol = 0)¶

int
sign::
Newton
(DistMatrix<F> &A, sign::Scaling scaling = FROB_NORM, int maxIts = 100, Base<F> tol = 0)¶ Runs a (scaled) Newton iteration for at most
maxIts
iterations with the specified tolerance, which, if undefined, is set to \(n \epsilon\), where \(n\) is the matrix dimension and \(\epsilon\) is the machine epsilon. The return value is the number of performed iterations.