Every square matrix 



A


{\displaystyle A}

 can be decomposed into a product of a lower triangular matrix 



L


{\displaystyle L}

 and a upper triangular matrix 



U


{\displaystyle U}

,
as described in LU decomposition.





A
=
L
U


{\displaystyle A=LU}


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:





A
=


(




a

11





a

12





a

13







a

21





a

22





a

23







a

31





a

32





a

33





)


=


(




l

11




0


0





l

21





l

22




0





l

31





l

32





l

33





)




(




u

11





u

12





u

13






0



u

22





u

23






0


0



u

33





)


=
L
U


{\displaystyle A={\begin{pmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}l_{11}&0&0\\l_{21}&l_{22}&0\\l_{31}&l_{32}&l_{33}\\\end{pmatrix}}{\begin{pmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}}=LU}


We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of 



L


{\displaystyle L}

 are set to 1






l

11


=
1


{\displaystyle l_{11}=1}







l

22


=
1


{\displaystyle l_{22}=1}







l

33


=
1


{\displaystyle l_{33}=1}


so we get a solvable system of 9 unknowns and 9 equations.





A
=


(




a

11





a

12





a

13







a

21





a

22





a

23







a

31





a

32





a

33





)


=


(



1


0


0





l

21




1


0





l

31





l

32




1



)




(




u

11





u

12





u

13






0



u

22





u

23






0


0



u

33





)


=


(




u

11





u

12





u

13







u

11



l

21





u

12



l

21


+

u

22





u

13



l

21


+

u

23







u

11



l

31





u

12



l

31


+

u

22



l

32





u

13



l

31


+

u

23



l

32


+

u

33





)


=
L
U


{\displaystyle A={\begin{pmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\\\end{pmatrix}}={\begin{pmatrix}1&0&0\\l_{21}&1&0\\l_{31}&l_{32}&1\\\end{pmatrix}}{\begin{pmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\end{pmatrix}}={\begin{pmatrix}u_{11}&u_{12}&u_{13}\\u_{11}l_{21}&u_{12}l_{21}+u_{22}&u_{13}l_{21}+u_{23}\\u_{11}l_{31}&u_{12}l_{31}+u_{22}l_{32}&u_{13}l_{31}+u_{23}l_{32}+u_{33}\end{pmatrix}}=LU}


Solving for the other 



l


{\displaystyle l}

 and 



u


{\displaystyle u}

, we get the following equations:






u

11


=

a

11




{\displaystyle u_{11}=a_{11}}







u

12


=

a

12




{\displaystyle u_{12}=a_{12}}







u

13


=

a

13




{\displaystyle u_{13}=a_{13}}







u

22


=

a

22


−

u

12



l

21




{\displaystyle u_{22}=a_{22}-u_{12}l_{21}}







u

23


=

a

23


−

u

13



l

21




{\displaystyle u_{23}=a_{23}-u_{13}l_{21}}







u

33


=

a

33


−
(

u

13



l

31


+

u

23



l

32


)


{\displaystyle u_{33}=a_{33}-(u_{13}l_{31}+u_{23}l_{32})}


and for 



l


{\displaystyle l}

:






l

21


=


1

u

11





a

21




{\displaystyle l_{21}={\frac {1}{u_{11}}}a_{21}}







l

31


=


1

u

11





a

31




{\displaystyle l_{31}={\frac {1}{u_{11}}}a_{31}}







l

32


=


1

u

22




(

a

32


−

u

12



l

31


)


{\displaystyle l_{32}={\frac {1}{u_{22}}}(a_{32}-u_{12}l_{31})}


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



U


{\displaystyle U}








u

i
j


=

a

i
j


−

∑

k
=
1


i
−
1



u

k
j



l

i
k




{\displaystyle u_{ij}=a_{ij}-\sum _{k=1}^{i-1}u_{kj}l_{ik}}


and then for 



L


{\displaystyle L}








l

i
j


=


1

u

j
j




(

a

i
j


−

∑

k
=
1


j
−
1



u

k
j



l

i
k


)


{\displaystyle l_{ij}={\frac {1}{u_{jj}}}(a_{ij}-\sum _{k=1}^{j-1}u_{kj}l_{ik})}


We see in the second formula that to get the 




l

i
j




{\displaystyle l_{ij}}

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




u

j
j




{\displaystyle u_{jj}}

, so we get problems when 




u

j
j




{\displaystyle u_{jj}}

 is either 0 or very small, which leads to numerical instability.

The solution to this problem is pivoting 



A


{\displaystyle A}

, which means rearranging the rows of 



A


{\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 



A


{\displaystyle A}

. Rearranging the rows means to multiply 



A


{\displaystyle A}

 by a permutation matrix 



P


{\displaystyle P}

:





P
A
⇒

A
′



{\displaystyle PA\Rightarrow A'}


Example:







(



0


1




1


0



)




(



1


4




2


3



)


⇒


(



2


3




1


4



)




{\displaystyle {\begin{pmatrix}0&1\\1&0\end{pmatrix}}{\begin{pmatrix}1&4\\2&3\end{pmatrix}}\Rightarrow {\begin{pmatrix}2&3\\1&4\end{pmatrix}}}


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





P
A
=
L
U


{\displaystyle PA=LU}




Task description

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



A


{\displaystyle A}

 and return a lower triangular matrix 



L


{\displaystyle L}

, a upper triangular matrix 



U


{\displaystyle U}

 and a permutation matrix 



P


{\displaystyle P}

,
so that the above equation is fulfilled.

You should then test it on the following two examples and include your output.



Example 1

A

1   3   5
2   4   7
1   1   0

L

1.00000   0.00000   0.00000
0.50000   1.00000   0.00000
0.50000  -1.00000   1.00000

U

2.00000   4.00000   7.00000
0.00000   1.00000   1.50000
0.00000   0.00000  -2.00000

P

0   1   0
1   0   0
0   0   1

Example 2

A

11    9   24    2
 1    5    2    6
 3   17   18    1
 2    5    7    1

L

1.00000   0.00000   0.00000   0.00000
0.27273   1.00000   0.00000   0.00000
0.09091   0.28750   1.00000   0.00000
0.18182   0.23125   0.00360   1.00000

U

11.00000    9.00000   24.00000    2.00000
 0.00000   14.54545   11.45455    0.45455
 0.00000    0.00000   -3.47500    5.68750
 0.00000    0.00000    0.00000    0.51079

P

1   0   0   0
0   0   1   0
0   1   0   0
0   0   0   1



#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;
}