Every square matrix <math>A</math> can be decomposed into a product of a lower triangular matrix <math>L</math> and a upper triangular matrix <math>U</math>,
as described in [[wp:LU decomposition|LU decomposition]].

:<math>A = LU</math>

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:

:<math>
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
</math>

We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of <math>L</math> are set to 1

:<math>l_{11}=1</math>
:<math>l_{22}=1</math>
:<math>l_{33}=1</math>

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

:<math>
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
</math>

Solving for the other <math>l</math> and <math>u</math>, we get the following equations:

:<math>u_{11}=a_{11}</math>
:<math>u_{12}=a_{12}</math>
:<math>u_{13}=a_{13}</math>

:<math>u_{22}=a_{22} - u_{12}l_{21}</math>
:<math>u_{23}=a_{23} - u_{13}l_{21}</math>

:<math>u_{33}=a_{33} - (u_{13}l_{31} + u_{23}l_{32})</math>

and for <math>l</math>:

:<math>l_{21}=\frac{1}{u_{11}} a_{21}</math>
:<math>l_{31}=\frac{1}{u_{11}} a_{31}</math>

:<math>l_{32}=\frac{1}{u_{22}} (a_{32} - u_{12}l_{31})</math>

We see that there is a calculation pattern, which can be expressed as the following formulas, first for <math>U</math>

:<math>u_{ij} = a_{ij} - \sum_{k=1}^{i-1} u_{kj}l_{ik}</math>

and then for <math>L</math>

:<math>l_{ij} = \frac{1}{u_{jj}} (a_{ij} - \sum_{k=1}^{j-1} u_{kj}l_{ik})</math>

We see in the second formula that to get the <math>l_{ij}</math> below the diagonal, we have to divide by the diagonal element (pivot) <math>u_{jj}</math>, so we get problems when <math>u_{jj}</math> is either 0 or very small, which leads to numerical instability.

The solution to this problem is ''pivoting'' <math>A</math>, which means rearranging the rows of <math>A</math>, prior to the <math>LU</math> decomposition, in a way that the largest element of each column gets onto the diagonal of <math>A</math>. Rearranging the rows means to multiply <math>A</math> by a permutation matrix <math>P</math>:

:<math>PA \Rightarrow A'</math>

Example:

:<math>
\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}
</math>

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

:<math>PA = LU</math>


'''Task description'''

The task is to implement a routine which will take a square nxn matrix <math>A</math> and return a lower triangular matrix <math>L</math>, a upper triangular matrix <math>U</math> and a permutation matrix <math>P</math>,
so that the above equation is fullfilled.
You should then test it on the following two examples and include your output.

Example 1:
<pre>
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
</pre>

Example 2:
<pre>
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
</pre>