luDecomposition.h

template<typename T>
Decomposition<T> Decompose(const Matrix<T> &matrix) {
    const size_t nbRows = matrix.rows();
    const size_t nbColumns = matrix.columns();
    if(nbRows != nbColumns) {
        throw std::domain_error("Matrix is not square.");
    }

    Decomposition<T> decomposition(matrix);

    for(size_t column = 0; column < nbColumns; ++column) {
        decomposition.L(column, column) = 1; //(1)

        for(size_t row = column + 1; row < nbRows; ++row) {
            const T & divisor = decomposition.U(column, column);
            if(std::fabs(divisor) < std::numeric_limits<T>::min()) {
                throw std::domain_error("Division by 0."); //a_ii != 0 is necessary because of pivoting with diaognal strategy
            }

            decomposition.L(row, column) = decomposition.U(row, column) / divisor; //(2)

            for(size_t col = column; col < nbColumns; ++col) {
                decomposition.U(row, col) -= decomposition.L(row, column) * decomposition.U(column, col); //(3)
            }
        }
    }

    return decomposition;
}