LU factorization

Implementation

Subroutines

Test driver

Example driver

Partial pivoting

Since LU factorization without pivoting is known to be unstable for general matrices, it is standard practice to pivot the rows of \(A\) during the factorization (this is called partial pivoting since the columns are not also pivoted). An LU factorization with partial pivoting therefore computes \(P\), \(L\), and \(U\) such that \(PA=LU\), where \(L\) and \(U\) are as described above and \(P\) is a permutation matrix.

Overwrites the matrix \(A\) with the LU decomposition of \(PA\), where \(P\) is represented by the permutation vector p, which consists of the column indices of the nonzero entry in each row of \(P\).

Factorization

C++ API

void LU(Matrix<F> &A, Matrix<int> &p)
void LU(AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &p)

C API

ElError ElLUPartialPiv_s(ElMatrix_s A, ElMatrix_i p)
ElError ElLUPartialPiv_d(ElMatrix_d A, ElMatrix_i p)
ElError ElLUPartialPiv_c(ElMatrix_c A, ElMatrix_i p)
ElError ElLUPartialPiv_z(ElMatrix_z A, ElMatrix_i p)
ElError ElLUPartialPivDist_s(ElDistMatrix_s A, ElDistMatrix_i p)
ElError ElLUPartialPivDist_d(ElDistMatrix_d A, ElDistMatrix_i p)
ElError ElLUPartialPivDist_c(ElDistMatrix_c A, ElDistMatrix_i p)
ElError ElLUPartialPivDist_z(ElDistMatrix_z A, ElDistMatrix_i p)

Solving linear systems with the factorization

C++ API

void lu::SolveAfter(Orientation orientation, const Matrix<F> &A, const Matrix<int> &p, Matrix<F> &B)
void lu::SolveAfter(Orientation orientation, const AbstractDistMatrix<F> &A, const AbstractDistMatrix<int> &p, AbstractDistMatrix<F> &B)

C API

ElError ElSolveAfterLUPartialPiv_s(ElOrientation orientation, ElConstMatrix_s A, ElConstMatrix_i p, ElMatrix_s B)
ElError ElSolveAfterLUPartialPiv_d(ElOrientation orientation, ElConstMatrix_d A, ElConstMatrix_i p, ElMatrix_d B)
ElError ElSolveAfterLUPartialPiv_c(ElOrientation orientation, ElConstMatrix_c A, ElConstMatrix_i p, ElMatrix_c B)
ElError ElSolveAfterLUPartialPiv_z(ElOrientation orientation, ElConstMatrix_z A, ElConstMatrix_i p, ElMatrix_z B)
ElError ElSolveAfterLUPartialPivDist_s(ElOrientation orientation, ElConstDistMatrix_s A, ElConstDistMatrix_i p, ElDistMatrix_s B)
ElError ElSolveAfterLUPartialPivDist_d(ElOrientation orientation, ElConstDistMatrix_d A, ElConstDistMatrix_i p, ElDistMatrix_d B)
ElError ElSolveAfterLUPartialPivDist_c(ElOrientation orientation, ElConstDistMatrix_c A, ElConstDistMatrix_i p, ElDistMatrix_c B)
ElError ElSolveAfterLUPartialPivDist_z(ElOrientation orientation, ElConstDistMatrix_z A, ElConstDistMatrix_i p, ElDistMatrix_z B)

Full pivoting

Overwrites the matrix \(A\) with the LU decomposition of \(PAQ^T\), where \(P\) and \(Q\) are represented by the permutation vectors p and q, which consist of the column indices of the nonzero entry in each row of \(P\) and \(Q\), respectively.

Factorization

C++ API

void LU(Matrix<F> &A, Matrix<int> &p, Matrix<int> &q)
void LU(AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &p, AbstractDistMatrix<F> &q)

C API

ElError ElLUFullPiv_s(ElMatrix_s A, ElMatrix_i p, ElMatrix_i q)
ElError ElLUFullPiv_d(ElMatrix_d A, ElMatrix_i p, ElMatrix_i q)
ElError ElLUFullPiv_c(ElMatrix_c A, ElMatrix_i p, ElMatrix_i q)
ElError ElLUFullPiv_z(ElMatrix_z A, ElMatrix_i p, ElMatrix_i q)
ElError ElLUFullPivDist_s(ElDistMatrix_s A, ElDistMatrix_i p, ElDistMatrix_i q)
ElError ElLUFullPivDist_d(ElDistMatrix_d A, ElDistMatrix_i p, ElDistMatrix_i q)
ElError ElLUFullPivDist_c(ElDistMatrix_c A, ElDistMatrix_i p, ElDistMatrix_i q)
ElError ElLUFullPivDist_z(ElDistMatrix_z A, ElDistMatrix_i p, ElDistMatrix_i q)

Solving linear systems with the factorization

C++ API

void lu::SolveAfter(Orientation orientation, const Matrix<F> &A, const Matrix<int> &p, const Matrix<int> &q, Matrix<F> &B)
void lu::SolveAfter(Orientation orientation, const AbstractDistMatrix<F> &A, const AbstractDistMatrix<int> &p, const AbstractDistMatrix<int> &q, AbstractDistMatrix<F> &B)

C API

ElError ElSolveAfterLUFullPiv_s(ElOrientation orientation, ElConstMatrix_s A, ElConstMatrix_i p, ElConstMatrix_i q, ElMatrix_s B)
ElError ElSolveAfterLUFullPiv_d(ElOrientation orientation, ElConstMatrix_d A, ElConstMatrix_i p, ElConstMatrix_i q, ElMatrix_d B)
ElError ElSolveAfterLUFullPiv_c(ElOrientation orientation, ElConstMatrix_c A, ElConstMatrix_i p, ElConstMatrix_i q, ElMatrix_c B)
ElError ElSolveAfterLUFullPiv_z(ElOrientation orientation, ElConstMatrix_z A, ElConstMatrix_i p, ElConstMatrix_i q, ElMatrix_z B)
ElError ElSolveAfterLUFullPivDist_s(ElOrientation orientation, ElConstDistMatrix_s A, ElConstDistMatrix_i p, ElConstDistMatrix_i q, ElDistMatrix_s B)
ElError ElSolveAfterLUFullPivDist_d(ElOrientation orientation, ElConstDistMatrix_d A, ElConstDistMatrix_i p, ElConstDistMatrix_i q, ElDistMatrix_d B)
ElError ElSolveAfterLUFullPivDist_c(ElOrientation orientation, ElConstDistMatrix_c A, ElConstDistMatrix_i p, ElConstDistMatrix_i q, ElDistMatrix_c B)
ElError ElSolveAfterLUFullPivDist_z(ElOrientation orientation, ElConstDistMatrix_z A, ElConstDistMatrix_i p, ElConstDistMatrix_i q, ElDistMatrix_z B)

No pivoting

Given \(A \in \mathbb{F}^{m \times n}\), an LU factorization (without pivoting) attempts to find a unit lower-trapezoidal \(L \in \mathbb{F}^{m \times \mbox{min}(m,n)}\) and upper-trapezoidal \(U \in \mathbb{F}^{\mbox{min}(m,n) \times n}\) such that \(A=LU\). Since \(L\) is required to have its diaganal entries set to one: the upper portion of \(A\) can be overwritten with U, and the strictly lower portion of \(A\) can be overwritten with the strictly lower portion of \(L\). If a numerically zero diagonal entry of \(U\) is created, then a SingularMatrixException will be thrown.

Note

It might be appropriate to switch this routine to a ZeroPivotException, as it is strange to refer to non-square matrices as singular.

Factorization

The following routines overwrite \(A\) with its LU decomposition.

C++ API

void LU(Matrix<F> &A)
void LU(AbstractDistMatrix<F> &A)

C API

ElError ElLU_s(ElMatrix_s A)
ElError ElLU_d(ElMatrix_d A)
ElError ElLU_c(ElMatrix_c A)
ElError ElLU_z(ElMatrix_z A)
ElError ElLUDist_s(ElDistMatrix_s A)
ElError ElLUDist_d(ElDistMatrix_d A)
ElError ElLUDist_c(ElDistMatrix_c A)
ElError ElLUDist_z(ElDistMatrix_z A)

Solving linear systems with the factorization

C++ API

void lu::SolveAfter(Orientation orientation, const Matrix<F> &A, Matrix<F> &B)
void lu::SolveAfter(Orientation orientation, const AbstractDistMatrix<F> &A, AbstractDistMatrix<F> &B)

C API

ElError ElSolveAfterLU_s(ElOrientation orientation, ElConstMatrix_s A, ElMatrix_s B)
ElError ElSolveAfterLU_d(ElOrientation orientation, ElConstMatrix_d A, ElMatrix_d B)
ElError ElSolveAfterLU_c(ElOrientation orientation, ElConstMatrix_c A, ElMatrix_c B)
ElError ElSolveAfterLU_z(ElOrientation orientation, ElConstMatrix_z A, ElMatrix_z B)
ElError ElSolveAfterLUDist_s(ElOrientation orientation, ElConstDistMatrix_s A, ElDistMatrix_s B)
ElError ElSolveAfterLUDist_d(ElOrientation orientation, ElConstDistMatrix_d A, ElDistMatrix_d B)
ElError ElSolveAfterLUDist_c(ElOrientation orientation, ElConstDistMatrix_c A, ElDistMatrix_c B)
ElError ElSolveAfterLUDist_z(ElOrientation orientation, ElConstDistMatrix_z A, ElDistMatrix_z B)

Rank-one modification to a factorization

Modify an existing LU factorization, \(A = P^T L U\), to incorporate the rank-one update \(A + u v^T\) or \(A + u v^H\). This algorithm only requires a quadratic number of operations.

Note

The current implementation has only been tested for square matrices.

C++ API

void LUMod(Matrix<F> &A, Matrix<int> &p, const Matrix<F> &u, const Matrix<F> &v, bool conjugate = true, Base<F> tau = 0.1)
void LUMod(AbstractDistMatrix<F> &A, AbstractDistMatrix<int> &p, const AbstractDistMatrix<F> &u, const AbstractDistMatrix<F> &v, bool conjugate = true, Base<F> tau = 0.1)

C API

ElError ElLUMod_s(ElMatrix_s A, ElMatrix_i p, ElConstMatrix_s u, ElConstMatrix_s v, float tau)
ElError ElLUMod_d(ElMatrix_d A, ElMatrix_i p, ElConstMatrix_d u, ElConstMatrix_d v, double tau)
ElError ElLUMod_c(ElMatrix_c A, ElMatrix_i p, ElConstMatrix_c u, ElConstMatrix_c v, float tau)
ElError ElLUMod_z(ElMatrix_z A, ElMatrix_i p, ElConstMatrix_z u, ElConstMatrix_z v, double tau)
ElError ElLUModDist_s(ElDistMatrix_s A, ElDistMatrix_i p, ElConstDistMatrix_s u, ElConstDistMatrix_s v, float tau)
ElError ElLUModDist_d(ElDistMatrix_d A, ElDistMatrix_i p, ElConstDistMatrix_d u, ElConstDistMatrix_d v, double tau)
ElError ElLUModDist_c(ElDistMatrix_c A, ElDistMatrix_i p, ElConstDistMatrix_c u, ElConstDistMatrix_c v, float tau)
ElError ElLUModDist_z(ElDistMatrix_z A, ElDistMatrix_i p, ElConstDistMatrix_z u, ElConstDistMatrix_z v, double tau)