LU factorization is a popular method to decompose a matrix as the product of a Lower triangle and an Upper triangle matrices. Lot of simulations perform the solving of a given system of linear equations. In general, the resolution of these systems is quite hard so matrix decomposition method is frequently used to subdivide the problem. As it's easier to solve a triangular matrix than a standard matrix, this is usually the first step on many simulations.

How it is done on the code side

Requirements:

• Should be a square matrix

This method performs a PLU because we are using a P pivot vector to use the most convenient rows during the process. Following the Gaussian elimination we are going to make zeros on the upper and lower areas of the matrix.

To reduce the error we are going to look for the absolute highest value on the current row:

```    // Iterate through each column
for ( unsigned int col=0; col<this->getColumnsCount()-1; col++ )
{
permute(col);```

permute function will find that value and interchange its row if it was necessary. Next, we start to iterate through the rows and on each row we compute the pivot for that row. This is done dividing the current row,col value by current step col,col value.

```        // Iterate through each row to do zero
for ( unsigned int row=col+1; row<this->getRowsCount(); row++ )
{
// Compute the pivot
T p = static_cast<T>(-static_cast<T>(this->get(row, col))/static_cast<T>(this->get(col, col)));```

Finally, we update that row adding the to that pair row,column the product of the pivot by the pair row,column step:

```            // Update row
for ( unsigned int k=col; k<this->getColumnsCount(); k++ )
{
this->set(row, k, this->get(row, k)+p*this->get(col, k));
}
// Store the pivot for L
this->set(row, col, -p);
}```

Putting all together:

LU method

```template <typename T>
void SquareMatrix<T>::lu()
{
// Iterate through each column
for ( unsigned int col=0; col<this->getColumnsCount()-1; col++ )
{
permute(col);
// Iterate through each row to do zero
for ( unsigned int row=col+1; row<this->getRowsCount(); row++ )
{
// Compute the pivot
T p = static_cast<T>(-static_cast<T>(this->get(row, col))/static_cast<T>(this->get(col, col)));
// Update row
for ( unsigned int k=col; k<this->getColumnsCount(); k++ )
{
this->set(row, k, this->get(row, k)+p*this->get(col, k));
}
// Store the pivot for L
this->set(row, col, -p);
}
}
}```

Permute method

```template <typename T>
void SquareMatrix<T>::permute(const unsigned int startRow)
{
T maxValue = static_cast<double>(std::abs(this->get(startRow, startRow)));
unsigned int maxValueRow = startRow;
// Find the absolute max value of a column in a rows range(startRow:nRows)
for ( unsigned int currRow=startRow+1; currRow<this->getColumnsCount(); currRow++ )
{
if ( std::abs(this->get(currRow, startRow)) > maxValue ) {
maxValue = std::abs(this->get(currRow, startRow));
maxValueRow = currRow;
}
}
// Interchange row (startRow <-> maxValueRow)
if ( maxValueRow != startRow ) {
T tmp;
for ( unsigned int k=0; k<this->getColumnsCount(); k++ )
{
tmp = this->get(startRow, k);
this->set(startRow, k, this->get(maxValueRow, k));
this->set(maxValueRow, k, tmp);
}
}
}```

The application

The application has been developed using Qt and it is able to factorize matrices with ranges between three and eight. It's very easy to use:

1. Set the proper range
2. Fill the values by yourself or using automatic generator
3. Factorize the matrix

All the code is available here.