Hermitian eigensolvers

Implementation

Test driver

Example driver

Elemental provides a collection of routines for both full and partial solutions 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.

Computing eigenvalues

Compute the set of eigenvalues of the Hermitian matrix A determined by the subset structure (the default is to compute all eigenvalues). By default, Elemental is used to transform the problem to and from real symmetric tridiagonal form and PMRRR is used to solve the tridiagonal EVP. Alternatively, the ctrl structure may be altered to request the usage of a (prototype) Spectral Divide and Conquer algorithm.

Note

Strict subset computation is not currently supported with Spectral Divide and Conquer.

C++ API

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

void HermitianEig(UpperOrLower uplo, AbstractDistMatrix<F> &A, AbstractDistMatrix<Base<F>> &w, SortType sort = ASCENDING, const HermitianEigSubset<Base<F>> subset = HermitianEigSubset<Base<F>>(), const HermitianEigCtrl<Base<F>> ctrl = HermitianEigCtrl<Base<F>>())

C API

ElError ElHermitianEig_s(ElUpperOrLower uplo, ElMatrix_s A, ElMatrix_s w, ElSortType sort)

ElError ElHermitianEig_d(ElUpperOrLower uplo, ElMatrix_d A, ElMatrix_d w, ElSortType sort)

ElError ElHermitianEig_c(ElUpperOrLower uplo, ElMatrix_c A, ElMatrix_s w, ElSortType sort)

ElError ElHermitianEig_z(ElUpperOrLower uplo, ElMatrix_z A, ElMatrix_d w, ElSortType sort)

ElError ElHermitianEigDist_s(ElUpperOrLower uplo, ElDistMatrix_s A, ElDistMatrix_s w, ElSortType sort)

ElError ElHermitianEigDist_d(ElUpperOrLower uplo, ElDistMatrix_d A, ElDistMatrix_d w, ElSortType sort)

ElError ElHermitianEigDist_c(ElUpperOrLower uplo, ElDistMatrix_c A, ElDistMatrix_s w, ElSortType sort)

ElError ElHermitianEigDist_z(ElUpperOrLower uplo, ElDistMatrix_z A, ElDistMatrix_d w, ElSortType sort)

Return the full set of eigenvalues

TODO: An expert interface which used HermitianEigCtrl

ElError ElHermitianEigPartial_s(ElUpperOrLower uplo, ElMatrix_s A, ElMatrix_s w, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPartial_d(ElUpperOrLower uplo, ElMatrix_d A, ElMatrix_d w, ElSortType sort, ElHermitianEigSubset_d subset)

ElError ElHermitianEigPartial_c(ElUpperOrLower uplo, ElMatrix_c A, ElMatrix_s w, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPartial_z(ElUpperOrLower uplo, ElMatrix_z A, ElMatrix_d w, ElSortType sort, ElHermitianEigSubset_d subset)

ElError ElHermitianEigPartialDist_s(ElUpperOrLower uplo, ElDistMatrix_s A, ElDistMatrix_s w, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPartialDist_d(ElUpperOrLower uplo, ElDistMatrix_d A, ElDistMatrix_d w, ElSortType sort, ElHermitianEigSubset_d subset)

ElError ElHermitianEigPartialDist_c(ElUpperOrLower uplo, ElDistMatrix_c A, ElDistMatrix_s w, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPartialDist_z(ElUpperOrLower uplo, ElDistMatrix_z A, ElDistMatrix_d w, ElSortType sort, ElHermitianEigSubset_d subset)

Return a subset of the eigenvalues

Computing eigenpairs

Compute the set of eigenpairs of the Hermitian matrix A determined by the subset structure (the default is to compute all eigenpairs). By default, Elemental is used to transform the problem to and from real symmetric tridiagonal form and PMRRR is used to solve the tridiagonal EVP. Alternatively, the ctrl structure may be altered to request and tune the usage of a (prototype) Spectral Divide and Conquer algorithm.

Note

Strict subset computation is not currently supported with Spectral Divide and Conquer.

C++ API

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

void HermitianEig(UpperOrLower uplo, AbstractDistMatrix<F> &A, AbstractDistMatrix<Base<F>> &w, AbstractDistMatrix<F> &Z, SortType sort = ASCENDING, const HermitianEigSubset<Base<F>> subset = HermitianEigSubset<Base<F>>(), const HermitianEigCtrl<Base<F>> ctrl = HermitianEigCtrl<Base<F>>())

C API

ElError ElHermitianEigPair_s(ElUpperOrLower uplo, ElMatrix_s A, ElMatrix_s w, ElMatrix_s Z, ElSortType sort)

ElError ElHermitianEigPair_d(ElUpperOrLower uplo, ElMatrix_d A, ElMatrix_d w, ElMatrix_d Z, ElSortType sort)

ElError ElHermitianEigPair_c(ElUpperOrLower uplo, ElMatrix_c A, ElMatrix_s w, ElMatrix_c Z, ElSortType sort)

ElError ElHermitianEigPair_z(ElUpperOrLower uplo, ElMatrix_z A, ElMatrix_d w, ElMatrix_z Z, ElSortType sort)

ElError ElHermitianEigPairDist_s(ElUpperOrLower uplo, ElDistMatrix_s A, ElDistMatrix_s w, ElDistMatrix_s Z, ElSortType sort)

ElError ElHermitianEigPairDist_d(ElUpperOrLower uplo, ElDistMatrix_d A, ElDistMatrix_d w, ElDistMatrix_d Z, ElSortType sort)

ElError ElHermitianEigPairDist_c(ElUpperOrLower uplo, ElDistMatrix_c A, ElDistMatrix_s w, ElDistMatrix_c Z, ElSortType sort)

ElError ElHermitianEigPairDist_z(ElUpperOrLower uplo, ElDistMatrix_z A, ElDistMatrix_d w, ElDistMatrix_z Z, ElSortType sort)

Return the full eigenvalue decomposition.

TODO: An expert interface which used HermitianEigCtrl

ElError ElHermitianEigPairPartial_s(ElUpperOrLower uplo, ElMatrix_s A, ElMatrix_s w, ElMatrix_s Z, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPairPartial_d(ElUpperOrLower uplo, ElMatrix_d A, ElMatrix_d w, ElMatrix_d Z, ElSortType sort, ElHermitianEigSubset_d subset)

ElError ElHermitianEigPairPartial_c(ElUpperOrLower uplo, ElMatrix_c A, ElMatrix_s w, ElMatrix_c Z, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPairPartial_z(ElUpperOrLower uplo, ElMatrix_z A, ElMatrix_d w, ElMatrix_z Z, ElSortType sort, ElHermitianEigSubset_d subset)

ElError ElHermitianEigPairPartialDist_s(ElUpperOrLower uplo, ElDistMatrix_s A, ElDistMatrix_s w, ElDistMatrix_s Z, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPairPartialDist_d(ElUpperOrLower uplo, ElDistMatrix_d A, ElDistMatrix_d w, ElDistMatrix_d Z, ElSortType sort, ElHermitianEigSubset_d subset)

ElError ElHermitianEigPairPartialDist_c(ElUpperOrLower uplo, ElDistMatrix_c A, ElDistMatrix_s w, ElDistMatrix_c Z, ElSortType sort, ElHermitianEigSubset_s subset)

ElError ElHermitianEigPairPartialDist_z(ElUpperOrLower uplo, ElDistMatrix_z A, ElDistMatrix_d w, ElDistMatrix_z Z, ElSortType sort, ElHermitianEigSubset_d subset)

Return a subset of the eigenpairs

Algorithmic options

The default approach starts with Householder tridiagonalization (ala HermitianTridiag()) and then calls Matthias Petschow and Paolo Bientinesi’s PMRRR for the tridiagonal eigenvalue problem.

Note

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

However, it is also possible to use a (prototype) Spectral Divide and Conquer algorithm (see, for example, Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue problem, Nakatsukasa et al., and Fast linear algebra is stable, Demmel et al.). In order to do so, the HermitianEigCtrl<Real> structure should be modified so that useSDC is true. The Spectral Divide and Conquer algorithm (if selected) is controlled via the HermitianSDCCtrl<Real> member structure.

C++ API

template<>
type HermitianEigCtrl<Real>

HermitianTridiagCtrl tridiagCtrl

HermitianSDCCtrl<Real> sdcCtrl

bool useSDC

HermitianEigCtrl()

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

template<>
type HermitianEigCtrl<Base<F>>

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

template<>
type HermitianSDCCtrl<Real>

Int cutoff

Int maxInnerIts

Int maxOuterIts

Real tol

Real spreadFactor

bool random

bool progress

HermitianSDCCtrl()

Defaults to using a sequential Schur decomposition for problems of size 256 or smaller (cutoff=256), a maximum of two random rank-revealing attempts per spectral split trials (maxInnerIts=2), and a maximum of ten trial spectral split trials before failure (maxOuterIts=10). The tolerance for accepting a split is kept at its default (via tol=0), the relative spreading factor for perturbing the estimate of the projected spectral median is set to :1e-6 (spreadFactor=1e-6), and progress is not displayed by default (progress=false).

template<>
type HermitianSDCCtrl<Base<F>>

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

C API

HermitianEigCtrl_s

ElHermitianTridiagCtrl tridiagCtrl

ElHermitianSDCCtrl_s sdcCtrl

bool useSDC

HermitianEigCtrl_d

ElHermitianTridiagCtrl tridiagCtrl

ElHermitianSDCCtrl_d sdcCtrl

bool useSDC

void ElHermitianEigCtrlDefault_s(ElHermitianEigCtrl_s* ctrl)

void ElHermitianEigCtrlDefault_d(ElHermitianEigCtrl_d* ctrl)

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

ElHermitianSDCCtrl_s

ElInt cutoff

ElInt maxInnerIts

ElInt maxOuterIts

float tol

float spreadFactor

bool random

bool progress

ElHermitianSDCCtrl_d

ElInt cutoff

ElInt maxInnerIts

ElInt maxOuterIts

double tol

double spreadFactor

bool random

bool progress

ElHermitianSDCCtrlDefault_s(ElHermitianSDCCtrl_s* ctrl)

ElHermitianSDCCtrlDefault_d(ElHermitianSDCCtrl_d* ctrl)

Defaults to using a sequential Schur decomposition for problems of size 256 or smaller (cutoff=256), a maximum of two random rank-revealing attempts per spectral split trials (maxInnerIts=2), and a maximum of ten trial spectral split trials before failure (maxOuterIts=10). The tolerance for accepting a split is kept at its default (via tol=0), the relative spreading factor for perturbing the estimate of the projected spectral median is set to :1e-6 (spreadFactor=1e-6), and progress is not displayed by default (progress=false).