Matrix decompositions

Hermitian eigensolver

Implementation

Test driver

Example driver

Elemental provides a collection of routines for both full and partial solution of the Hermitian eigenvalue problem

\[A Z = Z \Lambda,\]

where A is the given Hermitian matrix, and unitary Z and real diagonal \(\Lambda\) are sought. In particular, with the eigenvalues and corresponding eigenpairs labeled in non-decreasing order, the three basic modes are:

1. Compute all eigenvalues or eigenpairs, \(\{\lambda_i\}_{i=0}^{n-1}\) or \(\{(z_i,\lambda_i)\}_{i=0}^{n-1}\).

2. Compute the eigenvalues or eigenpairs with a given range of indices, say \(\{\lambda_i\}_{i=a}^b\) or \(\{(z_i,\lambda_i)\}_{i=a}^b\), with \(0 \le a \le b < n\).

3. Compute all eigenpairs (or just eigenvalues) with eigenvalues lying in a particular half-open interval, either \(\{\lambda_i \;|\; \lambda_i \in (a,b] \}\) or \(\{ (z_i,\lambda_i) \;|\; \lambda_i \in (a,b] \}\).

As of now, all three approaches start with Householder tridiagonalization (ala HermitianTridiag()) and then call Matthias Petschow and Paolo Bientinesi’s PMRRR for the tridiagonal eigenvalue problem.

Note

Unfortunately, PMRRR currently only supports double-precision problems, and so the parallel versions of these routines are limited to real and complex double-precision matrices.

template<>
class HermitianEigCtrl<Real>

HermitianTridiagCtrl tridiagCtrl

HermitianSdcCtrl<Real> sdcCtrl

bool useSdc

HermitianEigCtrl()

Initializes tridiagCtrl and sdcCtrl to their defaults and sets useSdc to false.

Note

Please see the Tuning parameters section for extensive information on maximizing the performance of Householder tridiagonalization.

Full spectrum computation

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the full set of eigenvalues of the Hermitian matrix A.

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Z, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the full set of eigenpairs of the Hermitian matrix A.

Index-based subset computation

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenvalues of a Hermitian matrix A with indices in the range \(a,a+1,...,b\).

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Z, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenpairs of a Hermitian matrix A with indices in the range \(a,a+1,...,b\).

Range-based subset computation

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F, STAR, STAR> &A, DistMatrix<Base<F>, STAR, STAR> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenvalues of a Hermitian matrix A lying in the half-open interval \((a,b]\).

void HermitianEig(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F, STAR, STAR> &A, DistMatrix<Base<F>, STAR, STAR> &w, DistMatrix<F, STAR, STAR> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianEig(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenpairs of a Hermitian matrix A with eigenvalues lying in the half-open interval \((a,b]\).

Spectral divide and conquer

Implementation

The primary references for this approach is Demmel et al.’s Fast linear algebra is stable and Nakatsukasa et al.’s Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue problem.

template<>
class HermitianSdcCtrl<Real>

int cutoff

int maxInnerIts

int maxOuterIts

Real tol

Real spreadFactor

bool random

bool progress

template<>
type HermitianSdcCtrl<Base<F>>

A particular case where the datatype is the base of the potentially complex type F.

void herm_eig::SDC(Matrix<F> &A, Matrix<Base<F>> &w, HermitianSdcCtrl<Base<F>> sdcCtrl = HermitianSdcCtrl<Base<F>>())

void herm_eig::SDC(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, HermitianSdcCtrl<Base<F>> sdcCtrl = HermitianSdcCtrl<Base<F>>())

Compute the eigenvalues of the matrix \(A\) via a QDWH-based spectral divide and conquer process.

void herm_eig::SDC(Matrix<F> &A, Matrix<Base<F>> &w, Matrix<F> &Q, HermitianSdcCtrl<Base<F>> sdcCtrl = HermitianSdcCtrl<Base<F>>())

void herm_eig::SDC(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Q, HermitianSdcCtrl<Base<F>> sdcCtrl = HermitianSdcCtrl<Base<F>>())

Attempt to also compute the eigenvectors.

Skew-Hermitian eigensolver

Essentially identical to the Hermitian eigensolver, HermitianEig(); for any skew-Hermitian matrix \(G\), \(iG\) is Hermitian, as

\[(iG)^H = -iG^H = iG.\]

This fact implies a fast method for solving skew-Hermitian eigenvalue problems:

1. Form \(iG\) in \(O(n^2)\) work (switching to complex arithmetic in the real case)

2. Run a Hermitian eigensolve on \(iG\), yielding \(iG=Z \Lambda Z^H\).

3. Recognize that \(G=Z (-i \Lambda) Z^H\) provides an EVD of \(G\).

Please see the HermitianEig() documentation for more details.

Note

Unfortunately, PMRRR currently only supports double-precision problems, and so the parallel versions of these routines are limited to real and complex double-precision matrices.

Implementation

Full spectrum computation

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the full set of eigenvalues of the skew-Hermitian matrix G.

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Matrix<Complex<Base<F>>> &Z, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, DistMatrix<Complex<Base<F>>> &Z, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the full set of eigenpairs of the skew-Hermitian matrix G.

Index-based subset computation

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenvalues of a skew-Hermitian matrix G with indices in the range \(a,a+1,...,b\).

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Matrix<Complex<Base<F>>> &Z, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, DistMatrix<Complex<Base<F>>> &Z, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenpairs of a skew-Hermitian matrix G with indices in the range \(a,a+1,...,b\).

Range-based subset computation

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenvalues of a skew-Hermitian matrix G lying in the half-open interval \((a,b]i\).

void SkewHermitianEig(UpperOrLower uplo, Matrix<F> &G, Matrix<Base<F>> &wImag, Matrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void SkewHermitianEig(UpperOrLower uplo, DistMatrix<F> &G, DistMatrix<Base<F>, VR, STAR> &wImag, DistMatrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenpairs of a skew-Hermitian matrix G with eigenvalues lying in the half-open interval \((a,b]i\).

Hermitian generalized-definite eigensolvers

The following Hermitian generalized-definite eigenvalue problems frequently appear, where both \(A\) and \(B\) are Hermitian, and \(B\) is additionally positive-definite.

enum HermitianGenDefiniteEigType

enumerator ABX

\[A B x = \lambda x\]

enumerator BAX

\[B A x = \lambda x\]

enumerator AXBX

\[A x = B x \lambda\]

Implementation

Test driver

Full spectrum computation

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the full set of eigenvalues of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Matrix<Base<F>> &Z, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<Base<F>> &Z, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the full set of eigenpairs of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

Index-based subset computation

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenvalues with indices in the range \(a,a+1,...,b\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Matrix<F> &Z, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, int a, int b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenpairs with indices in the range \(a,a+1,...,b\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

Range-based subset computation

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenvalues lying in the half-open interval \((a,b]\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, Matrix<F> &A, Matrix<F> &B, Matrix<Base<F>> &w, Matrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

void HermitianGenDefiniteEig(HermitianGenDefiniteEigType type, UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &B, DistMatrix<Base<F>, VR, STAR> &w, DistMatrix<F> &Z, Base<F> a, Base<F> b, SortType sort = UNSORTED, const HermitianEigCtrl<Base<F>> &ctrl = HermitianEigCtrl<Base<F>>())

Compute the eigenpairs whose eigenvalues lie in the half-open interval \((a,b]\) of a generalized EVP involving the Hermitian matrices A and B, where B is also positive-definite.

Unitary eigensolver

Not yet written, will likely be based on Ming Gu’s unitary Divide and Conquer algorithm for unitary Hessenberg matrices. The spectral divide and conquer technique described below should suffice in the meantime.

Normal eigensolver

Not yet written, will likely be based on Angelika Bunse-Gerstner et al.’s Jacobi-like method for simultaneous diagonalization of the commuting Hermitian and skew-Hermitian portions of the matrix. The spectral divide and conquer scheme described below should suffice in the meantime.

Schur decomposition

Elemental contains a native prototype implementation of a spectral divide and conquer scheme for the Schur decomposition, but it is not yet robust enough to handle general matrices. For local matrices, Elemental defaults to calling LAPACK’s Hessenberg QR algorithm (with Aggressive Early Deflation); if support for ScaLAPACK was detected during configuration, Elemental defaults to ScaLAPACK’s Hessenberg QR algorithm (without deflation), otherwise the Spectral Divide and Conquer approach is attempted.

Main header file

void Schur(Matrix<F> &A)

void Schur(DistMatrix<F> &A)

void Schur(Matrix<F> &A, Matrix<F> &Q)

void Schur(DistMatrix<F> &A, DistMatrix<F> &Q)

Hessenberg QR algorithm

QR algorithm header file

void schur::QR(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, bool fullTriangle = false)

void schur::QR(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, Matrix<F> &Q, bool fullTriangle = true)

Condense the matrix to upper-Hessenberg form and then use a sequential QR algorithm to compute a (partial) Schur decomposition. It is optional whether or not the full Schur factor is computed.

void schur::QR(DistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, bool fullTriangle = false, bool aed = false)

void schur::QR(BlockDistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, bool fullTriangle = false, bool aed = false)

void schur::QR(DistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, DistMatrix<F> &Q, bool fullTriangle = true, bool aed = false)

void schur::QR(BlockDistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, BlockDistMatrix<F> &Q, bool fullTriangle = true, bool aed = false)

Condense the matrix to upper-Hessenberg form and then use ScaLAPACK’s parallel QR algorithm to compute a (partial) Schur decomposition. It is optional whether or not the full Schur factor is computed, and Aggressive Early Deflation is also optional for real matrices (as of now, its usage is not recommended due to known bugs in the implementation).

Quasi-triangular manipulation

void schur::QuasiTriangEig(const Matrix<F> &U, Matrix<Complex<Base<F>>> &w)

void schur::QuasiTriangEig(const DistMatrix<F> &U, DistMatrix<Complex<Base<F>>, colDist, rowDist> &w)

Return the eigenvalues of the upper quasi-triangular matrix U in the vector w.

Matrix<Complex<Base<F>>> schur::QuasiTriangEig(const Matrix<F> &U)

DistMatrix<Complex<Base<F>>, VR, STAR> schur::QuasiTriangEig(const DistMatrix<F> &U)

Return the eigenvalues of the upper quasi-triangular matrix U.

void schur::QuasiTriangEig(const Matrix<F> &dMain, const Matrix<F> &dSub, const Matrix<F> &dSup, Matrix<Complex<Base<F>>> &w)

The underlying computation routine for computing the eigenvalues of quasi-triangular matrices. The vectors dMain, dSub, and dSup should respectively contain the main, sub, and super-diagonals of the upper quasi-triangular matrix.

void schur::RealToComplex(const Matrix<Real> &UQuasi, Matrix<Complex<Real>> &U)

void schur::RealToComplex(const DistMatrix<Real> &UQuasi, DistMatrix<Complex<Real>> &U)

Rotate a real upper quasi-triangular matrix into a complex upper triangular matrix.

void schur::CheckRealSchur(const Matrix<Real> &U, bool standardForm = false)

void schur::CheckRealSchur(const DistMatrix<Real> &U, bool standardForm = false)

Check whether or not the largest diagonal blocks of the upper quasi-triangular matrix are at most \(2 \times 2\) and, optionally, check if the \(2 \times 2\) diagonal blocks are in standard form (if so, their diagonal must be constant and the product of the off-diagonal entries should be negative).

Spectral divide and conquer

The primary reference for this approach is Demmel et al.’s Fast linear algebra is stable. While the current implementation needs a large number of algorithmic improvements, especially with respect to choosing the Mobius transformations, it tends to succeed on random matrices.

SDC header file

template<>
class SdcCtrl<Real>

int cutoff

int maxInnerIts

int maxOuterIts

Real tol

Real spreadFactor

bool random

bool progress

template<>
type SdcCtrl<Base<F>>

A particular case where the datatype is the base of the potentially complex type F.

void schur::SDC(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, bool formATR = false, SdcCtrl<Base<F>> sdcCtrl = SdcCtrl<Base<F>>())

void schur::SDC(DistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, bool formATR = false, SdcCtrl<Base<F>> sdcCtrl = SdcCtrl<Base<F>>())

Compute the eigenvalues of the matrix \(A\) via a spectral divide and conquer process. On exit, the eigenvalues of \(A\) will be stored on its diagonal, and, if formATR was set to true, the upper triangle of \(A\) will be its corresponding upper-triangular Schur factor.

void schur::SDC(Matrix<F> &A, Matrix<Complex<Base<F>>> &w, Matrix<F> &Q, bool formATR = true, SdcCtrl<Base<F>> sdcCtrl = SdcCtrl<Base<F>>())

void schur::SDC(DistMatrix<F> &A, DistMatrix<Complex<Base<F>>, VR, STAR> &w, DistMatrix<F> &Q, bool formATR = true, SdcCtrl<Base<F>> sdcCtrl = SdcCtrl<Base<F>>())

Attempt to also compute the Schur vectors.

Hermitian SVD

Given an eigenvalue decomposition of a Hermitian matrix \(A\), say

\[A = V \Lambda V^H,\]

where \(V\) is unitary and \(\Lambda\) is diagonal and real. Then an SVD of \(A\) can easily be computed as

\[A = U |\Lambda| V^H,\]

where the columns of \(U\) equal the columns of \(V\), modulo sign flips introduced by negative eigenvalues.

Implementation

void HermitianSVD(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &U, Matrix<F> &V)

void HermitianSVD(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &U, DistMatrix<F> &V)

Return a vector of singular values, \(s\), and the left and right singular vector matrices, \(U\) and \(V\), such that \(A=U \mathrm{diag}(s) V^H\).

void HermitianSVD(UpperOrLower uplo, Matrix<F> &A, Matrix<Base<F>> &s)

void HermitianSVD(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s)

Return the singular values of \(A\) in s. Note that the appropriate triangle of A is overwritten during computation.

Polar decomposition

Every matrix \(A\) can be written as \(A=QP\), where \(Q\) is unitary and \(P\) is Hermitian and positive semi-definite. This is known as the polar decomposition of \(A\) and can be constructed as \(Q := U V^H\) and \(P := V \Sigma V^H\), where \(A = U \Sigma V^H\) is the SVD of \(A\). Alternatively, it can be computed through (a dynamically-weighted) Halley iteration.

Implementation

SVD approach

QWDH approach

void Polar(Matrix<F> &A)

void Polar(DistMatrix<F> &A)

void Polar(Matrix<F> &A, Matrix<F> &P)

void Polar(DistMatrix<F> &A, DistMatrix<F> &P)

Compute the polar decomposition of \(A\), \(A=QP\), returning \(Q\) within A and \(P\) within P. The current implementation first computes the SVD.

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

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

void HermitianPolar(UpperOrLower uplo, Matrix<F> &A, Matrix<F> &P)

void HermitianPolar(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &P)

Compute the polar decomposition through a Hermitian EVD. Since this is equivalent to a Hermitian sign decomposition (if \(\text{sgn}(0)\) is set to 1), these routines are equivalent to HermitianSign.

polar namespace

int polar::QDWH(Matrix<F> &A, bool colPiv = false, int maxits = 20)

int polar::QDWH(DistMatrix<F> &A, bool colPiv = false, int maxIts = 20)

int hermitian_polar::QDWH(UpperOrLower uplo, Matrix<F> &A, bool colPiv = false, int maxits = 20)

int hermitian_polar::QDWH(UpperOrLower uplo, DistMatrix<F> &A, bool colPiv = false, int maxIts = 20)

Overwrites \(A\) with the \(Q\) from the polar decomposition using a QR-based dynamically weighted Halley iteration. The number of iterations used is returned upon completion. TODO: reference to Yuji’s paper

int polar::QDWH(Matrix<F> &A, Matrix<F> &P, bool colPiv = false, int maxits = 20)

int polar::QDWH(DistMatrix<F> &A, DistMatrix<F> &P, bool colPiv = false, int maxIts = 20)

int hermitian_polar::QDWH(UpperOrLower uplo, Matrix<F> &A, Matrix<F> &P, bool colPiv = false, int maxits = 20)

int hermitian_polar::QDWH(UpperOrLower uplo, DistMatrix<F> &A, DistMatrix<F> &P, bool colPiv = false, int maxIts = 20)

Return the full polar decomposition, where \(P\) is HPD.

SVD

Given a general matrix \(A\), the Singular Value Decomposition is the triplet \((U,\Sigma,V)\) such that

\[A = U \Sigma V^H,\]

where \(U\) and \(V\) are unitary, and \(\Sigma\) is diagonal with non-negative entries.

Implementation

Subroutines

void SVD(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V)

void SVD(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V)

Overwrites A with \(U\), s with the diagonal entries of \(\Sigma\), and V with \(V\).

void SVD(Matrix<F> &A, Matrix<Base<F>> &s)

void SVD(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s)

Forms the singular values of \(A\) in s. Note that A is overwritten in order to compute the singular values.

svd namespace

void svd::QRSVD(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V)

SVD which uses bidiagonal QR algorithm.

void svd::DivideAndConquerSVD(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V)

SVD which uses a bidiagonal divide-and-conquer algorithm.

void svd::Chan(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, double heightRatio = 1.2)

void svd::Chan(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V, double heightRatio = 1.5)

SVD which preprocesses with an initial QR decomposition if the matrix is sufficiently tall relative to its width.

void svd::GolubReinschUpper(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s)

void svd::GolubReinschUpper(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V)

Computes the singular values (and vectors) of a matrix which is taller than it is wide using the Golub-Reinsch algorithm, though DQDS is used when only the singular values are sought.

void svd::Thresholded(Matrix<F> &A, Matrix<Base<F>> &s, Matrix<F> &V, Base<F> tol = 0, bool relative = false)

void svd::Thresholded(DistMatrix<F> &A, DistMatrix<Base<F>, VR, STAR> &s, DistMatrix<F> &V, Base<F> tol = 0, bool relative = false)

Computes the singular triplets whose singular values are larger than a specified tolerance using the cross-product algorithm. This is often advantageous because tridiagonal eigensolvers tend to enjoy better parallel implementations than bidiagonal SVD’s.