LDL factorization

Implementation

Subroutines

Test driver

Example driver

Symmetric pivoting

The following routines return a pivoted LDL factorization, where the vector \(p\) contains the column indices of the nonzero entries of the permutation matrix \(P\) such that \(PAP^T\) equals either \(LDL^T\) or \(LDL^H\), where \(D\) is quasi-diagonal. The Bunch-Kaufman pivoting rules are used within a higher-performance blocked algorithm, whereas the Bunch-Parlett strategy uses an unblocked, but typically more accurate, algorithm.

Factorization

C++ API

type LDLPivotType

An enum for specifying the symmetric pivoting strategy. The current (not yet all supported) options include:

  • BUNCH_KAUFMAN_A

  • BUNCH_KAUFMAN_C (not yet supported)

  • BUNCH_KAUFMAN_D

  • BUNCH_KAUFMAN_BOUNDED (not yet supported)

  • BUNCH_PARLETT

  • LDL_WITHOUT_PIVOTING

template<>
Real LDLPivotConstant<Real>(LDLPivotType pivotType)

Maps various LDL pivotings schemes to their optimal threshold constant.

template<>
type LDLPivotCtrl<Real>
LDLPivotType pivotType

The type of pivoting to perform (by default, BUNCH_KAUFMAN_A)

Real gamma

Pivot tolerance (by default, set to LDLPivotConstant<Real>(pivotType))

LDLPivotCtrl(LDLPivotType piv = BUNCH_KAUFMAN_A)
type LDLPivot
Int nb

Whether the pivot is 1x1 or 2x2.

Int from[2]

The source indices of the row or rows to swap with for the 1x1 or 2x2 pivot.

void LDL(Matrix<F> &A, Matrix<F> &dSub, Matrix<int> &p, bool conjugate = false, const LDLPivotCtrl<Base<F>> &ctrl = LDLPivotCtrl<Base<F>>())
void LDL(AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &dSub, AbstractDistMatrix<int> &p, bool conjugate = false, const LDLPivotCtrl<Base<F>> &ctrl = LDLPivotCtrl<Base<F>>())

C API

ElLDLPivotType

An enum for specifying the symmetric pivoting strategy. The current (not yet all supported) options include:

  • EL_BUNCH_KAUFMAN_A

  • EL_BUNCH_KAUFMAN_C (not yet supported)

  • EL_BUNCH_KAUFMAN_D

  • EL_BUNCH_KAUFMAN_BOUNDED (not yet supported)

  • EL_BUNCH_PARLETT

ElLDLPivotCtrl_s
ElLDLPivotType pivotType
float gamma
ElLDLPivotCtrl_d
ElLDLPivotType pivotType
double gamma
ElLDLPivot
Int nb

Whether the pivot is 1x1 or 2x2.

Int from[2]

The source indices of the row or rows to swap with for the 1x1 or 2x2 pivot.

ElError ElLDLPiv_s(ElMatrix_s A, ElMatrix_s dSub, ElMatrix_i p)
ElError ElLDLPiv_d(ElMatrix_d A, ElMatrix_d dSub, ElMatrix_i p)
ElError ElLDLPiv_c(ElMatrix_c A, ElMatrix_c dSub, ElMatrix_i p, bool conjugate)
ElError ElLDLPiv_z(ElMatrix_z A, ElMatrix_z dSub, ElMatrix_i p, bool conjugate)
ElError ElLDLPivDist_s(ElDistMatrix_s A, ElDistMatrix_s dSub, ElDistMatrix_i p)
ElError ElLDLPivDist_d(ElDistMatrix_d A, ElDistMatrix_d dSub, ElDistMatrix_i p)
ElError ElLDLPivDist_c(ElDistMatrix_c A, ElDistMatrix_c dSub, ElDistMatrix_i p, bool conjugate)
ElError ElLDLPivDist_z(ElDistMatrix_z A, ElDistMatrix_z dSub, ElDistMatrix_i p, bool conjugate)

TODO:: Document expert versions

Python API

TODO: Document

Solving linear systems after factorization

C++ API

void ldl::SolveAfter(const Matrix<F> &A, const Matrix<F> &dSub, const Matrix<int> &p, Matrix<F> &B, bool conjugated = false)
void ldl::SolveAfter(const AbstractDistMatrix<F> &A, const AbstractDistMatrix<F> &dSub, const AbstractDistMatrix<int> &p, AbstractDistMatrix<F> &B, bool conjugated = false)

C API

ElError ElSolveAfterLDLPiv_s(ElConstMatrix_s A, ElConstMatrix_s dSub, ElConstMatrix_i p, ElMatrix_s B)
ElError ElSolveAfterLDLPiv_d(ElConstMatrix_d A, ElConstMatrix_d dSub, ElConstMatrix_i p, ElMatrix_d B)
ElError ElSolveAfterLDLPiv_c(ElConstMatrix_c A, ElConstMatrix_c dSub, ElConstMatrix_i p, ElMatrix_c B, bool conjugate)
ElError ElSolveAfterLDLPiv_z(ElConstMatrix_z A, ElConstMatrix_z dSub, ElConstMatrix_i p, ElMatrix_z B, bool conjugate)
ElError ElSolveAfterLDLPivDist_s(ElConstDistMatrix_s A, ElConstDistMatrix_s dSub, ElConstDistMatrix_i p, ElDistMatrix_s B)
ElError ElSolveAfterLDLPivDist_d(ElConstDistMatrix_d A, ElConstDistMatrix_d dSub, ElConstDistMatrix_i p, ElDistMatrix_d B)
ElError ElSolveAfterLDLPivDist_c(ElConstDistMatrix_c A, ElConstDistMatrix_c dSub, ElConstDistMatrix_i p, ElDistMatrix_c B, bool conjugate)
ElError ElSolveAfterLDLPivDist_z(ElConstDistMatrix_z A, ElConstDistMatrix_z dSub, ElConstDistMatrix_i p, ElDistMatrix_z B, bool conjugate)

Python API

TODO: Document

Multiply vectors after factorization

C++ API

void ldl::MultiplyAfter(const Matrix<F> &A, const Matrix<F> &dSub, const Matrix<int> &p, Matrix<F> &B, bool conjugated = false)
void ldl::MultiplyAfter(const AbstractDistMatrix<F> &A, const AbstractDistMatrix<F> &dSub, const AbstractDistMatrix<int> &p, AbstractDistMatrix<F> &B, bool conjugated = false)

C API

ElError ElMultiplyAfterLDLPiv_s(ElConstMatrix_s A, ElConstMatrix_s dSub, ElConstMatrix_i p, ElMatrix_s B)
ElError ElMultiplyAfterLDLPiv_d(ElConstMatrix_d A, ElConstMatrix_d dSub, ElConstMatrix_i p, ElMatrix_d B)
ElError ElMultiplyAfterLDLPiv_c(ElConstMatrix_c A, ElConstMatrix_c dSub, ElConstMatrix_i p, ElMatrix_c B, bool conjugate)
ElError ElMultiplyAfterLDLPiv_z(ElConstMatrix_z A, ElConstMatrix_z dSub, ElConstMatrix_i p, ElMatrix_z B, bool conjugate)
ElError ElMultiplyAfterLDLPivDist_s(ElConstDistMatrix_s A, ElConstDistMatrix_s dSub, ElConstDistMatrix_i p, ElDistMatrix_s B)
ElError ElMultiplyAfterLDLPivDist_d(ElConstDistMatrix_d A, ElConstDistMatrix_d dSub, ElConstDistMatrix_i p, ElDistMatrix_d B)
ElError ElMultiplyAfterLDLPivDist_c(ElConstDistMatrix_c A, ElConstDistMatrix_c dSub, ElConstDistMatrix_i p, ElDistMatrix_c B, bool conjugate)
ElError ElMultiplyAfterLDLPivDist_z(ElConstDistMatrix_z A, ElConstDistMatrix_z dSub, ElConstDistMatrix_i p, ElDistMatrix_z B, bool conjugate)

Computing inertia after factorization

C++ API

InertiaType ldl::Inertia(const Matrix<Base<F>> &d, const Matrix<F> &dSub)
InertiaType ldl::Inertia(const AbstractDistMatrix<Base<F>> &d, const AbstractDistMatrix<F> &dSub)

C API

ElError ElInertiaAfterLDL_s(ElConstMatrix_s d, ElConstMatrix_s dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDL_d(ElConstMatrix_d d, ElConstMatrix_d dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDL_c(ElConstMatrix_s d, ElConstMatrix_c dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDL_z(ElConstMatrix_d d, ElConstMatrix_z dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDLDist_s(ElConstDistMatrix_s d, ElConstDistMatrix_s dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDLDist_d(ElConstDistMatrix_d d, ElConstDistMatrix_d dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDLDist_c(ElConstDistMatrix_s d, ElConstDistMatrix_c dSub, ElInertiaType* inertia)
ElError ElInertiaAfterLDLDist_z(ElConstDistMatrix_d d, ElConstDistMatrix_z dSub, ElInertiaType* inertia)

No pivoting

Though Cholesky factorization is ideal for HPD matrices, unpivoted LDL factorization naturally applies to a slightly larger, but harder to define, class of matrices. Upon successful completion of the factorization, a lower-triangular (with unit diagonal) \(L\) and diagonal matrix \(D\), such that \(A = L D L^H\) or \(A = L D L^T\), will be returned in the lower triangle of \(A\). If a zero pivot is attempted, then a ZeroPivotException will be thrown.

Factorization

Warning

Please use the following routines with caution, as pivoting should be employed in most cases.

C++ API

void LDLT(Matrix<F> &A)
void LDLT(AbstractDistMatrix<F> &A)
void LDLH(Matrix<F> &A)
void LDLH(AbstractDistMatrix<F> &A)

C API

ElError ElLDLT_s(ElMatrix_s A)
ElError ElLDLT_d(ElMatrix_d A)
ElError ElLDLT_c(ElMatrix_c A)
ElError ElLDLT_z(ElMatrix_z A)
ElError ElLDLTDist_s(ElDistMatrix_s A)
ElError ElLDLTDist_d(ElDistMatrix_d A)
ElError ElLDLTDist_c(ElDistMatrix_c A)
ElError ElLDLTDist_z(ElDistMatrix_z A)
ElError ElLDLH_c(ElMatrix_c A)
ElError ElLDLH_z(ElMatrix_z A)
ElError ElLDLHDist_c(ElDistMatrix_c A)
ElError ElLDLHDist_z(ElDistMatrix_z A)

Solve linear systems after factorization

C++ API

void ldl::SolveAfter(const Matrix<F> &A, Matrix<F> &B, bool conjugated = false)
void ldl::SolveAfter(const AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &B, bool conjugated = false)

C API

ElError ElSolveAfterLDL_s(ElConstMatrix_s A, ElMatrix_s B)
ElError ElSolveAfterLDL_d(ElConstMatrix_d A, ElMatrix_d B)
ElError ElSolveAfterLDL_c(ElConstMatrix_c A, ElMatrix_c B, bool conjugate)
ElError ElSolveAfterLDL_z(ElConstMatrix_z A, ElMatrix_z B, bool conjugate)
ElError ElSolveAfterLDLDist_s(ElConstDistMatrix_s A, ElDistMatrix_s B)
ElError ElSolveAfterLDLDist_d(ElConstDistMatrix_d A, ElDistMatrix_d B)
ElError ElSolveAfterLDLDist_c(ElConstDistMatrix_c A, ElDistMatrix_c B, bool conjugate)
ElError ElSolveAfterLDLDist_z(ElConstDistMatrix_z A, ElDistMatrix_z B, bool conjugate)

Multiply vectors after factorization

C++ API

void ldl::MultiplyAfter(const Matrix<F> &A, Matrix<F> &B, bool conjugated = false)
void ldl::MultiplyAfter(const AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &B, bool conjugated = false)

C API

ElError ElMultiplyAfterLDL_s(ElConstMatrix_s A, ElMatrix_s B)
ElError ElMultiplyAfterLDL_d(ElConstMatrix_d A, ElMatrix_d B)
ElError ElMultiplyAfterLDL_c(ElConstMatrix_c A, ElMatrix_c B, bool conjugate)
ElError ElMultiplyAfterLDL_z(ElConstMatrix_z A, ElMatrix_z B, bool conjugate)
ElError ElMultiplyAfterLDLDist_s(ElConstDistMatrix_s A, ElDistMatrix_s B)
ElError ElMultiplyAfterLDLDist_d(ElConstDistMatrix_d A, ElDistMatrix_d B)
ElError ElMultiplyAfterLDLDist_c(ElConstDistMatrix_c A, ElDistMatrix_c B, bool conjugate)
ElError ElMultiplyAfterLDLDist_z(ElConstDistMatrix_z A, ElDistMatrix_z B, bool conjugate)