How to resolve the algorithm LU decomposition step by step in the C++ programming language

Problem Statement
Step by Step Solution
Sourcecode

Problem Statement

Every square matrix

{\displaystyle A}

can be decomposed into a product of a lower triangular matrix

{\displaystyle L}

and a upper triangular matrix

{\displaystyle U}

, as described in LU decomposition. It is a modified form of Gaussian elimination. While the Cholesky decomposition only works for symmetric, positive definite matrices, the more general LU decomposition works for any square matrix. There are several algorithms for calculating L and U. To derive Crout's algorithm for a 3x3 example, we have to solve the following system: We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of

{\displaystyle L}

are set to 1 so we get a solvable system of 9 unknowns and 9 equations. Solving for the other

{\displaystyle l}

and

{\displaystyle u}

, we get the following equations: and for

{\displaystyle l}

: We see that there is a calculation pattern, which can be expressed as the following formulas, first for

{\displaystyle U}

and then for

{\displaystyle L}

We see in the second formula that to get the

i j

{\displaystyle l_{ij}}

below the diagonal, we have to divide by the diagonal element (pivot)

j j

{\displaystyle u_{jj}}

, so we get problems when

j j

{\displaystyle u_{jj}}

is either 0 or very small, which leads to numerical instability. The solution to this problem is pivoting

{\displaystyle A}

, which means rearranging the rows of

{\displaystyle A}

, prior to the

L U

{\displaystyle LU}

decomposition, in a way that the largest element of each column gets onto the diagonal of

{\displaystyle A}

. Rearranging the rows means to multiply

{\displaystyle A}

by a permutation matrix

{\displaystyle P}

: Example: The decomposition algorithm is then applied on the rearranged matrix so that

The task is to implement a routine which will take a square nxn matrix

{\displaystyle A}

and return a lower triangular matrix

{\displaystyle L}

, a upper triangular matrix

{\displaystyle U}

and a permutation matrix

{\displaystyle P}

, so that the above equation is fulfilled. You should then test it on the following two examples and include your output.

Let's start with the solution:

Step by Step solution about How to resolve the algorithm LU decomposition step by step in the C++ programming language

The provided C++ code implements LU decomposition for a matrix. LU decomposition is a matrix factorization technique that decomposes a matrix into a lower triangular matrix and an upper triangular matrix. The code works for matrices with both double and integer elements.

Key Features of the Code:

matrix class:
- Represents a matrix with customizable rows, columns, and element type.
- Initializes matrices using various constructors, including from initializer lists.
- Provides access to elements using () operator with range checks.
print function:
- Prints a matrix in a visually appealing format using Unicode box drawing characters.
- Includes alignment and precision control.
lu_decompose function:
- Decomposes a square matrix (where the number of rows equals the number of columns) into lower (L), upper (U), and pivot (P) matrices using Gaussian elimination with partial pivoting.
- Throws a runtime error if the matrix is singular (non-invertible).
show_lu_decomposition function:
- Displays the original matrix (A), its LU decomposition results (L, U, and P), and handles exceptions during decomposition.

How the Code Works:

Define a matrix input using an initializer list to represent the matrix to be decomposed.
Call lu_decompose(input) to perform LU decomposition and obtain the L, U, and P matrices.
Use the print function to display the original matrix, L, U, and P matrices in a formatted way.
Handle any exceptions thrown during decomposition, such as singular matrices.

Example Usage:

The provided examples in the main function illustrate the LU decomposition of matrices with various sizes and values. The output is a formatted representation of the original matrix A, followed by its decomposition into L, U, and P matrices.

Exception Handling:

The lu_decompose function throws a runtime error if the input matrix is singular (non-invertible). This check prevents division by zero errors during decomposition.

Source code in the cpp programming language

#include <cassert>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <sstream>
#include <vector>

template <typename scalar_type> class matrix {
public:
    matrix(size_t rows, size_t columns)
        : rows_(rows), columns_(columns), elements_(rows * columns) {}

    matrix(size_t rows, size_t columns, scalar_type value)
        : rows_(rows), columns_(columns), elements_(rows * columns, value) {}

    matrix(size_t rows, size_t columns,
        const std::initializer_list<std::initializer_list<scalar_type>>& values)
        : rows_(rows), columns_(columns), elements_(rows * columns) {
        assert(values.size() <= rows_);
        size_t i = 0;
        for (const auto& row : values) {
            assert(row.size() <= columns_);
            std::copy(begin(row), end(row), &elements_[i]);
            i += columns_;
        }
    }

    size_t rows() const { return rows_; }
    size_t columns() const { return columns_; }

    const scalar_type& operator()(size_t row, size_t column) const {
        assert(row < rows_);
        assert(column < columns_);
        return elements_[row * columns_ + column];
    }
    scalar_type& operator()(size_t row, size_t column) {
        assert(row < rows_);
        assert(column < columns_);
        return elements_[row * columns_ + column];
    }
private:
    size_t rows_;
    size_t columns_;
    std::vector<scalar_type> elements_;
};

template <typename scalar_type>
void print(std::wostream& out, const matrix<scalar_type>& a) {
    const wchar_t* box_top_left = L"\x23a1";
    const wchar_t* box_top_right = L"\x23a4";
    const wchar_t* box_left = L"\x23a2";
    const wchar_t* box_right = L"\x23a5";
    const wchar_t* box_bottom_left = L"\x23a3";
    const wchar_t* box_bottom_right = L"\x23a6";

    const int precision = 5;
    size_t rows = a.rows(), columns = a.columns();
    std::vector<size_t> width(columns);
    for (size_t column = 0; column < columns; ++column) {
        size_t max_width = 0;
        for (size_t row = 0; row < rows; ++row) {
            std::ostringstream str;
            str << std::fixed << std::setprecision(precision) << a(row, column);
            max_width = std::max(max_width, str.str().length());
        }
        width[column] = max_width;
    }
    out << std::fixed << std::setprecision(precision);
    for (size_t row = 0; row < rows; ++row) {
        const bool top(row == 0), bottom(row + 1 == rows);
        out << (top ? box_top_left : (bottom ? box_bottom_left : box_left));
        for (size_t column = 0; column < columns; ++column) {
            if (column > 0)
                out << L' ';
            out << std::setw(width[column]) << a(row, column);
        }
        out << (top ? box_top_right : (bottom ? box_bottom_right : box_right));
        out << L'\n';
    }
}

// Return value is a tuple with elements (lower, upper, pivot)
template <typename scalar_type>
auto lu_decompose(const matrix<scalar_type>& input) {
    assert(input.rows() == input.columns());
    size_t n = input.rows();
    std::vector<size_t> perm(n);
    std::iota(perm.begin(), perm.end(), 0);
    matrix<scalar_type> lower(n, n);
    matrix<scalar_type> upper(n, n);
    matrix<scalar_type> input1(input);
    for (size_t j = 0; j < n; ++j) {
        size_t max_index = j;
        scalar_type max_value = 0;
        for (size_t i = j; i < n; ++i) {
            scalar_type value = std::abs(input1(perm[i], j));
            if (value > max_value) {
                max_index = i;
                max_value = value;
            }
        }
        if (max_value <= std::numeric_limits<scalar_type>::epsilon())
            throw std::runtime_error("matrix is singular");
        if (j != max_index)
            std::swap(perm[j], perm[max_index]);
        size_t jj = perm[j];
        for (size_t i = j + 1; i < n; ++i) {
            size_t ii = perm[i];
            input1(ii, j) /= input1(jj, j);
            for (size_t k = j + 1; k < n; ++k)
                input1(ii, k) -= input1(ii, j) * input1(jj, k);
        }
    }
    
    for (size_t j = 0; j < n; ++j) {
        lower(j, j) = 1;
        for (size_t i = j + 1; i < n; ++i)
            lower(i, j) = input1(perm[i], j);
        for (size_t i = 0; i <= j; ++i)
            upper(i, j) = input1(perm[i], j);
    }
    
    matrix<scalar_type> pivot(n, n);
    for (size_t i = 0; i < n; ++i)
        pivot(i, perm[i]) = 1;

    return std::make_tuple(lower, upper, pivot);
}

template <typename scalar_type>
void show_lu_decomposition(const matrix<scalar_type>& input) {
    try {
        std::wcout << L"A\n";
        print(std::wcout, input);
        auto result(lu_decompose(input));
        std::wcout << L"\nL\n";
        print(std::wcout, std::get<0>(result));
        std::wcout << L"\nU\n";
        print(std::wcout, std::get<1>(result));
        std::wcout << L"\nP\n";
        print(std::wcout, std::get<2>(result));
    } catch (const std::exception& ex) {
        std::cerr << ex.what() << '\n';
    }
}

int main() {
    std::wcout.imbue(std::locale(""));
    std::wcout << L"Example 1:\n";
    matrix<double> matrix1(3, 3,
       {{1, 3, 5},
        {2, 4, 7},
        {1, 1, 0}});
    show_lu_decomposition(matrix1);
    std::wcout << '\n';

    std::wcout << L"Example 2:\n";
    matrix<double> matrix2(4, 4,
      {{11, 9, 24, 2},
        {1, 5, 2, 6},
        {3, 17, 18, 1},
        {2, 5, 7, 1}});
    show_lu_decomposition(matrix2);
    std::wcout << '\n';
    
    std::wcout << L"Example 3:\n";
    matrix<double> matrix3(3, 3,
      {{-5, -6, -3},
       {-1,  0, -2},
       {-3, -4, -7}});
    show_lu_decomposition(matrix3);
    std::wcout << '\n';
    
    std::wcout << L"Example 4:\n";
    matrix<double> matrix4(3, 3,
      {{1, 2, 3},
       {4, 5, 6},
       {7, 8, 9}});
    show_lu_decomposition(matrix4);

    return 0;
}

You may also check:How to resolve the algorithm Address of a variable step by step in the Smalltalk programming language
You may also check:How to resolve the algorithm Bitcoin/public point to address step by step in the Racket programming language
You may also check:How to resolve the algorithm Abelian sandpile model/Identity step by step in the Wren programming language
You may also check:How to resolve the algorithm Array concatenation step by step in the Little programming language
You may also check:How to resolve the algorithm Pentagram step by step in the jq programming language

How to resolve the algorithm LU decomposition step by step in the C++ programming language

Table of Contents

Problem Statement

Step by Step solution about How to resolve the algorithm LU decomposition step by step in the C++ programming language

Source code in the cpp programming language