#include <linalg/RowPivotLU.hpp>

Detailed Description

LU factorization with row pivoting.

Factorizes a square matrix into a lower triangular and an upper-triangular factor.

This version uses row pivoting, but it is not rank-revealing.

Definition at line 16 of file RowPivotLU.hpp.

Collaboration diagram for RowPivotLU:

Constructors
	RowPivotLU ()=default
	Default constructor. More...

	RowPivotLU (const SquareMatrix &matrix)
	Factorize the given matrix. More...

	RowPivotLU (SquareMatrix &&matrix)
	Factorize the given matrix. More...

Factorization
void	compute (SquareMatrix &&matrix)
	Perform the LU factorization of the given matrix. More...

void	compute (const SquareMatrix &matrix)
	Perform the LU factorization of the given matrix. More...

Retrieving the L factor
SquareMatrix &&	steal_L ()
	Get the lower-triangular matrix L, reusing the internal storage. More...

void	get_L_inplace (Matrix &L) const
	Copy the lower-triangular matrix L to the given matrix. More...

SquareMatrix	get_L () const &
	Get a copy of the lower-triangular matrix L. More...

SquareMatrix &&	get_L () &&
	Get the lower-triangular matrix L. More...

Retrieving the U factor
SquareMatrix &&	steal_U ()
	Get the upper-triangular matrix U, reusing the internal storage. More...

void	get_U_inplace (Matrix &U) const
	Copy the upper-triangular matrix U to the given matrix. More...

SquareMatrix	get_U () const &
	Get a copy of the upper-triangular matrix U. More...

SquareMatrix &&	get_U () &&
	Get the upper-triangular matrix U. More...

Retrieving the P factor
PermutationMatrix &&	steal_P ()
	Get the permutation matrix P, reusing the internal storage. More...

PermutationMatrix	get_P () const &
	Get a copy of the permutation matrix P. More...

PermutationMatrix &&	get_P () &&
	Get the permutation matrix P. More...

Solving systems of equations problems
void	solve_inplace (Matrix &B) const
	Solve the system AX = B or LUX = B. More...

Matrix	solve (const Matrix &B) const
	Solve the system AX = B or LUX = B. More...

Matrix &&	solve (Matrix &&B) const
	Solve the system AX = B or LUX = B. More...

Vector	solve (const Vector &B) const
	Solve the system Ax = b or LUx = b. More...

Vector &&	solve (Vector &&B) const
	Solve the system Ax = b or LUx = b. More...

Access to internal representation
bool	is_factored () const
	Check if this object contains a factorization. More...

bool	has_LU () const
	Check if this object contains valid L and U factors. More...

bool	has_P () const
	Check if this object contains a valid permutation matrix P. More...

SquareMatrix &&	steal_LU ()
	Get the internal storage of the upper-triangular matrix U and the strict lower-triangular part of matrix L. More...

const SquareMatrix &	get_LU () const &
	Get a copy of the internal storage of the upper-triangular matrix U and the strict lower-triangular part of matrix L. More...

SquareMatrix &&	get_LU () &&
	Get the internal storage of the upper-triangular matrix U and the strict lower-triangular part of matrix L. More...

Private Types
enum	State { NotFactored = 0 , Factored = 1 }

Private Member Functions
void	compute_factorization ()
	The actual LU factorization algorithm. More...

void	back_subs (const Matrix &B, Matrix &X) const
	Back substitution algorithm for solving upper-triangular systems UX = B. More...

void	forward_subs (const Matrix &B, Matrix &X) const
	Forward substitution algorithm for solving lower-triangular systems LX = B. More...

Private Attributes
SquareMatrix	LU
	Result of a LU factorization: stores the upper-triangular matrix U and the strict lower-triangular part of matrix L. More...

PermutationMatrix	P = PermutationMatrix::RowPermutation
	The permutation of A that maximizes pivot size. More...

enum RowPivotLU::State	state = NotFactored

bool	valid_LU = false

bool	valid_P = false

Related Functions
(Note that these are not member functions.)
std::ostream &	operator<< (std::ostream &os, const RowPivotLU &lu)
	Print the L and U matrices of an LU object. More...

Member Enumeration Documentation

◆ State

enum State

private

Enumerator
NotFactored
Factored

Definition at line 163 of file RowPivotLU.hpp.

Constructor & Destructor Documentation

◆ RowPivotLU() [1/3]

RowPivotLU ( )

default

Default constructor.

◆ RowPivotLU() [2/3]

RowPivotLU ( const SquareMatrix & matrix )

inline

Factorize the given matrix.

Definition at line 24 of file RowPivotLU.hpp.

◆ RowPivotLU() [3/3]

RowPivotLU ( SquareMatrix && matrix )

inline

Factorize the given matrix.

Definition at line 26 of file RowPivotLU.hpp.

Member Function Documentation

◆ compute() [1/2]

void compute ( SquareMatrix && matrix )

Perform the LU factorization of the given matrix.

Definition at line 7 of file RowPivotLU.ipp.

◆ compute() [2/2]

void compute ( const SquareMatrix & matrix )

Perform the LU factorization of the given matrix.

Definition at line 14 of file RowPivotLU.ipp.

◆ steal_L()

SquareMatrix && steal_L ( )

Get the lower-triangular matrix L, reusing the internal storage.

Warning: After calling this function, the LU object is no longer valid, because this function steals its storage. Stealing both L and P is allowed (if you do not steal U, because it shares storage with L).

Definition at line 21 of file RowPivotLU.ipp.

◆ get_L_inplace()

void get_L_inplace ( Matrix & L ) const

Copy the lower-triangular matrix L to the given matrix.

Definition at line 36 of file RowPivotLU.ipp.

◆ get_L() [1/2]

SquareMatrix get_L ( ) const &

Get a copy of the lower-triangular matrix L.

Definition at line 52 of file RowPivotLU.ipp.

◆ get_L() [2/2]

SquareMatrix&& get_L ( ) &&

inline

Get the lower-triangular matrix L.

Definition at line 57 of file RowPivotLU.hpp.

◆ steal_U()

SquareMatrix && steal_U ( )

Get the upper-triangular matrix U, reusing the internal storage.

Warning: After calling this function, the LU object is no longer valid, because this function steals its storage. Stealing both U and P is allowed (if you do not steal L, because it shares storage with U).

Definition at line 58 of file RowPivotLU.ipp.

◆ get_U_inplace()

void get_U_inplace ( Matrix & U ) const

Copy the upper-triangular matrix U to the given matrix.

Definition at line 71 of file RowPivotLU.ipp.

◆ get_U() [1/2]

SquareMatrix get_U ( ) const &

Get a copy of the upper-triangular matrix U.

Definition at line 85 of file RowPivotLU.ipp.

◆ get_U() [2/2]

SquareMatrix&& get_U ( ) &&

inline

Get the upper-triangular matrix U.

Definition at line 77 of file RowPivotLU.hpp.

◆ steal_P()

PermutationMatrix && steal_P ( )

Get the permutation matrix P, reusing the internal storage.

Warning: After calling this function, the LU object is no longer valid, because this function steals its storage. Stealing P and either L or U (not both) is allowed.

Definition at line 91 of file RowPivotLU.ipp.

◆ get_P() [1/2]

PermutationMatrix get_P ( ) const &

inline

Get a copy of the permutation matrix P.

Definition at line 92 of file RowPivotLU.hpp.

◆ get_P() [2/2]

PermutationMatrix&& get_P ( ) &&

inline

Get the permutation matrix P.

Definition at line 94 of file RowPivotLU.hpp.

◆ solve_inplace()

void solve_inplace ( Matrix & B ) const

Solve the system AX = B or LUX = B.

Matrix B is overwritten with the result X.

Implementation

void RowPivotLU::solve_inplace(Matrix &B) const {
    // Solve the system AX = B, PAX = PB or LUX = PB.
    //
    // Let UX = Z, and first solve LZ = PB, which is a simple lower-triangular
    // system of equations.
    // Now that Z is known, solve UX = Z, which is a simple upper-triangular
    // system of equations.
    assert(is_factored());
 
    P.permute_rows(B);
    forward_subs(B, B); // overwrite B with Z
    back_subs(B, B);    // overwrite B (Z) with X
}

Definition at line 214 of file RowPivotLU.cpp.

◆ solve() [1/4]

Matrix solve ( const Matrix & B ) const

Solve the system AX = B or LUX = B.

Definition at line 103 of file RowPivotLU.ipp.

◆ solve() [2/4]

Matrix && solve ( Matrix && B ) const

Solve the system AX = B or LUX = B.

Definition at line 109 of file RowPivotLU.ipp.

◆ solve() [3/4]

Vector solve ( const Vector & B ) const

Solve the system Ax = b or LUx = b.

Definition at line 114 of file RowPivotLU.ipp.

◆ solve() [4/4]

Vector && solve ( Vector && B ) const

Solve the system Ax = b or LUx = b.

Definition at line 118 of file RowPivotLU.ipp.

◆ is_factored()

bool is_factored ( ) const

inline

Check if this object contains a factorization.

Definition at line 121 of file RowPivotLU.hpp.

◆ has_LU()

bool has_LU ( ) const

inline

Check if this object contains valid L and U factors.

Definition at line 124 of file RowPivotLU.hpp.

◆ has_P()

bool has_P ( ) const

inline

Check if this object contains a valid permutation matrix P.

Definition at line 127 of file RowPivotLU.hpp.

◆ steal_LU()

SquareMatrix && steal_LU ( )

Get the internal storage of the upper-triangular matrix U and the strict lower-triangular part of matrix L.

Warning: After calling this function, the LU object is no longer valid, because this function steals its storage. Stealing both LU and P is allowed (but not L or U individually).

Definition at line 97 of file RowPivotLU.ipp.

◆ get_LU() [1/2]

const SquareMatrix& get_LU ( ) const &

inline

Get a copy of the internal storage of the upper-triangular matrix U and the strict lower-triangular part of matrix L.

Definition at line 139 of file RowPivotLU.hpp.

◆ get_LU() [2/2]

SquareMatrix&& get_LU ( ) &&

inline

Get the internal storage of the upper-triangular matrix U and the strict lower-triangular part of matrix L.

Definition at line 142 of file RowPivotLU.hpp.

◆ compute_factorization()

void compute_factorization ( )

private

The actual LU factorization algorithm.

Precondition: LU contains the matrix A to be factorized; P contains the identity matrix (no permutations); LU.rows() == LU.cols(); P.size() == LU.rows()

Postcondition: The complete upper-triangular part of LU contains the full upper-triangular matrix U and the strict lower-triangular part of matrix L. The diagonal elements of L are implicitly 1.; get_L() * get_U() == get_P() * A (up to rounding errors)

See also: NoPivotLU::compute_factorization() for an introduction to the LU factorization. The basics presented there will not be reiterated here.

Implementation

void RowPivotLU::compute_factorization() {
    // For the intermediate calculations, we'll be working with LU.
    // It is initialized to the square n×n matrix to be factored.
 
    assert(LU.rows() == LU.cols());
    assert(P.size() == LU.rows());
 
    // The goal of the LU factorization algorithm is to repeatedly apply
    // transformations Lₖ to the matrix A to eventually end up with an upper-
    // triangular matrix U. When row pivoting is used, the rows of A are
    // permuted using a permutation matrix P:
    //
    //     Lₙ'⋯L₂'L₁'PA = U
    //
    // The main steps of the algorithm are exactly the same as the original
    // LU algorithm explained in LU.cpp, and will not be repeated here.
    // The only difference is that instead of using the diagonal element as the
    // pivot, rows are swapped so that the element with the largest magnitude
    // ends up on the diagonal and can be used as the pivot.
 
    // Loop over all columns of A:
    for (size_t k = 0; k < LU.cols(); ++k) {
        // In the following comments, k = [1, n], because this is more intuitive
        // and it follows the usual mathematical convention.
        // In the code, however, array indices start at zero, so k = [0, n-1].
 
        // On each iteration, the largest element on or below the diagonal in
        // the current (k-th) column will be used as the pivot.
        // To this end, the k-th row and the row that contains the largest
        // element are swapped, and the swapping is stored in the permutation
        // matrix, so that it can later be undone, when solving systems of
        // equations for example.
 
        // Find the largest element (in absolute value)
        double max_elem = std::abs(LU(k, k));
        size_t max_index = k;
        for (size_t i = k + 1; i < LU.rows(); ++i) {
            double abs_elem = std::abs(LU(i, k));
            if (abs_elem > max_elem) {
                max_elem = abs_elem;
                max_index = i;
            }
        }
 
        // Select the index of the element that is largest in absolute value as
        // the new pivot index.
        // If this index is not the diagonal element, rows have to be swapped:
        if (max_index != k) {
            P(k) = max_index;           // save the permutation
            LU.swap_rows(k, max_index); // actually perfrom the permutation
        }
 
        // Note how all columns of the two rows are permuted, not just the
        // columns greater than k. You might wonder how that'll ever work out
        // correctly in the end.
        // Recall that for the LU factorization without pivoting, the result was
        //     Lₙ⋯L₂L₁A = U.
        // When pivoting is used, however, the matrix is permuted before each
        // elimination step:
        //
        //     LₙPₙ⋯L₂P₂L₁P₁A = U
        //
        // Luckily, the product can be reordered, and all permutation matrices
        // Pₖ can be grouped together.
        // Without loss of generality, consider the pivoted LU factorization of
        // a 3×3 matrix:
        //
        //     L₃P₃L₂P₂L₁P₁A = U
        //
        // Now introduce the following permuted matrices Lₖ':
        //
        //     L₃' = L₃
        //     L₂' = P₃L₂P₃⁻¹
        //     L₁' = P₃P₂L₁P₂⁻¹P₃⁻¹
        //
        // You can then easily see that the following equation holds:
        //
        //     L₃'L₂'L₁'P₂P₃P₁A = U
        //   ⇔ L₃(P₃L₂P₃⁻¹)(P₃P₂L₁P₂⁻¹P₃⁻¹)P₃P₂P₁A = U
        //   ⇔ L₃P₃L₂P₂L₁P₁A = U
        //
        // Furthermore, the matrices Lₖ' have the same structure as Lₖ, because
        // only rows below the k-th pivot are permuted.
        //
        // All of this allows us to group all row permutations into a single
        // permutation matrix P, and use the same algorithm and storage format
        // as for the unpivoted case.
        //
        // The factors Lₖ' are computed implicitly by applying the row
        // permutations to the entire matrix that stores both the U and L
        // factors, rather than just to the elements of the trailing submatrix.
 
        // The rest of the algorithm is identical to the one explained in
        // NoPivotLU.cpp.
 
        double pivot = LU(k, k);
 
        // Compute the k-th column of L, the coefficients lᵢₖ:
        for (size_t i = k + 1; i < LU.rows(); ++i)
            LU(i, k) /= pivot;
 
        // Update the trailing submatrix A'(k+1:n,k+1:n) = LₖA(k+1:n,k+1:n):
        for (size_t c = k + 1; c < LU.cols(); ++c)
            // Subtract lᵢₖ times the current pivot row A(k,:):
            for (size_t i = k + 1; i < LU.rows(); ++i)
                // A'(i,c) = 1·A(i,c) - lᵢₖ·A(k,c)
                LU(i, c) -= LU(i, k) * LU(k, c);
 
        // Because of the row pivoting, zero pivots are no longer an issue,
        // since the pivot is always chosen to be the largest possible element.
        // When the matrix is singular, the algorithm will still fail, of
        // course.
    }
    state = Factored;
    valid_LU = true;
    valid_P = true;
}

Definition at line 25 of file RowPivotLU.cpp.

◆ back_subs()

void back_subs	(	const Matrix &	B,
		Matrix &	X
	)		const

private

Back substitution algorithm for solving upper-triangular systems UX = B.

Implementation

void RowPivotLU::back_subs(const Matrix &B, Matrix &X) const {
    // Solve upper triangular system UX = B by solving each column of B as a
    // vector system Uxᵢ = bᵢ
    //
    //     ┌                 ┐┌     ┐   ┌     ┐
    //     │ u₁₁ u₁₂ u₁₃ u₁₄ ││ x₁ᵢ │   │ b₁ᵢ │
    //     │     u₂₂ u₂₃ u₂₄ ││ x₂ᵢ │ = │ b₂ᵢ │
    //     │         u₃₃ u₃₄ ││ x₃ᵢ │   │ b₃ᵢ │
    //     │             u₄₄ ││ x₄ᵢ │   │ b₄ᵢ │
    //     └                 ┘└     ┘   └     ┘
    //
    // b₄ᵢ = u₄₄·x₄ᵢ                     ⟺ x₄ᵢ = b₄ᵢ/u₄₄
    // b₃ᵢ = u₃₃·x₃ᵢ + u₃₄·x₄ᵢ           ⟺ x₃ᵢ = (b₃ᵢ - u₃₄·x₄ᵢ)/u₃₃
    // b₂ᵢ = u₂₂·x₂ᵢ + u₂₃·x₃ᵢ + u₂₄·x₄ᵢ ⟺ x₂ᵢ = (b₂ᵢ - u₂₃·x₃ᵢ + u₂₄·x₄ᵢ)/u₂₂
    // ...
 
    for (size_t i = 0; i < B.cols(); ++i) {
        for (size_t r = LU.rows(); r-- > 0;) {
            X(r, i) = B(r, i);
            for (size_t c = r + 1; c < LU.cols(); ++c)
                X(r, i) -= LU(r, c) * X(c, i);
            X(r, i) /= LU(r, r);
        }
    }
}

Definition at line 149 of file RowPivotLU.cpp.

◆ forward_subs()

void forward_subs	(	const Matrix &	B,
		Matrix &	X
	)		const

private

Forward substitution algorithm for solving lower-triangular systems LX = B.

Implementation

void RowPivotLU::forward_subs(const Matrix &B, Matrix &X) const {
    // Solve lower triangular system LX = B by solving each column of B as a
    // vector system Lxᵢ = bᵢ.
    // The diagonal is always 1, due to the construction of the L matrix in the
    // LU algorithm.
    //
    //     ┌               ┐┌     ┐   ┌     ┐
    //     │ 1             ││ x₁ᵢ │   │ b₁ᵢ │
    //     │ l₂₁ 1         ││ x₂ᵢ │ = │ b₂ᵢ │
    //     │ l₃₁ l₃₂ 1     ││ x₃ᵢ │   │ b₃ᵢ │
    //     │ l₄₁ l₄₂ l₄₃ 1 ││ x₄ᵢ │   │ b₄ᵢ │
    //     └               ┘└     ┘   └     ┘
    //
    // b₁ᵢ =   1·x₁ᵢ                   ⟺ x₁ᵢ = b₁ᵢ
    // b₂ᵢ = l₂₁·x₁ᵢ +   1·x₂ᵢ         ⟺ x₂ᵢ = b₂ᵢ - l₂₁·x₁ᵢ
    // b₃ᵢ = l₃₁·x₁ᵢ + l₃₂·x₂ᵢ + 1·x₃ᵢ ⟺ x₃ᵢ = b₃ᵢ - l₃₂·x₂ᵢ - l₃₁·x₁ᵢ
    // ...
 
    for (size_t i = 0; i < B.cols(); ++i) {
        for (size_t r = 0; r < LU.rows(); ++r) {
            X(r, i) = B(r, i);
            for (size_t c = 0; c < r; ++c)
                X(r, i) -= LU(r, c) * X(c, i);
        }
    }
}

Definition at line 181 of file RowPivotLU.cpp.

Friends And Related Function Documentation

◆ operator<<()

std::ostream & operator<<	(	std::ostream &	os,
		const RowPivotLU &	lu
	)

Definition at line 125 of file RowPivotLU.ipp.

Member Data Documentation

◆ LU

SquareMatrix LU

private

Result of a LU factorization: stores the upper-triangular matrix U and the strict lower-triangular part of matrix L.

The diagonal elements of L are implicitly 1.

Definition at line 159 of file RowPivotLU.hpp.

◆ P

PermutationMatrix P = PermutationMatrix::RowPermutation

private

The permutation of A that maximizes pivot size.

Definition at line 161 of file RowPivotLU.hpp.

◆ state

enum RowPivotLU::State state = NotFactored

private

◆ valid_LU

bool valid_LU = false

private

Definition at line 168 of file RowPivotLU.hpp.

◆ valid_P

bool valid_P = false

private

Definition at line 169 of file RowPivotLU.hpp.

The documentation for this class was generated from the following files:

Detailed Description

Constructors

Factorization

Retrieving the L factor

Retrieving the U factor

Retrieving the P factor

Solving systems of equations problems

Access to internal representation

Private Types

Private Member Functions

Private Attributes

Related Functions

Member Enumeration Documentation

◆ State

Constructor & Destructor Documentation

◆ RowPivotLU() [1/3]

◆ RowPivotLU() [2/3]

◆ RowPivotLU() [3/3]

Member Function Documentation

◆ compute() [1/2]

◆ compute() [2/2]

◆ steal_L()

◆ get_L_inplace()

◆ get_L() [1/2]

◆ get_L() [2/2]

◆ steal_U()

◆ get_U_inplace()

◆ get_U() [1/2]

◆ get_U() [2/2]

◆ steal_P()

◆ get_P() [1/2]

◆ get_P() [2/2]

◆ solve_inplace()

Implementation

◆ solve() [1/4]

◆ solve() [2/4]

◆ solve() [3/4]

◆ solve() [4/4]

◆ is_factored()

◆ has_LU()

◆ has_P()

◆ steal_LU()

◆ get_LU() [1/2]

◆ get_LU() [2/2]

◆ compute_factorization()

Implementation

◆ back_subs()

Implementation

◆ forward_subs()

Implementation

Friends And Related Function Documentation

◆ operator<<()

Member Data Documentation

◆ LU

◆ P

◆ state

◆ valid_LU

◆ valid_P