pivotLUDecomposition.h

namespace PivotLUDecomposition {
    template<typename T>
    struct Decomposition {
        Matrix<T> P;
        Matrix<T> L;
        Matrix<T> U;

        Decomposition(const Matrix<T> &matrix) : P(matrix.rows(), matrix.columns()), L(matrix.rows(), matrix.columns()), U(matrix) {}
    };

    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);

        decomposition.P = MatrixFactory::IdentityMatrix<T>(nbRows);

        for(size_t k = 0; k < nbRows; ++k) {
            T max = 0;
            size_t pk = 0;
            for(size_t i = k; i < nbRows; ++i) {
                T s = 0;
                for(size_t j = k; j < nbColumns; ++j) {
                    s += std::fabs(decomposition.U(i, j)); //(1)
                }
                T q = std::fabs(decomposition.U(i, k)) / s; //(2)
                if(q > max) {
                    max = q;
                    pk = i;
                }
            }

            if(std::fabs(max) < std::numeric_limits<T>::min()) {
                throw std::domain_error("Pivot has 0 value.");
            }

            if(pk != k) {
                for(size_t j = 0; j < nbColumns; ++j) { //(3)
                    std::swap(decomposition.P(k, j), decomposition.P(pk, j));
                    std::swap(decomposition.L(k, j), decomposition.L(pk, j));
                    std::swap(decomposition.U(k, j), decomposition.U(pk, j));
                }
            }

            for(size_t i = k+1; i < nbRows; ++i) {
                decomposition.L(i, k) = decomposition.U(i, k)/decomposition.U(k, k); //(4)

                for(size_t j = k; j < nbColumns; ++j) {
                    decomposition.U(i, j) = decomposition.U(i, j) - decomposition.L(i, k) * decomposition.U(k, j); //(5)
                }
            }
        }

        for(size_t k = 0; k < nbRows; ++k) {
            decomposition.L(k,k) = 1;
        }

        return decomposition;
    }
}