choleskyDecomposition.h

namespace CholeskyDecomposition {
    template<typename T>
    Matrix<T> SimplifySymmetricMatrix(Matrix<T> matrix) {
        const size_t nbRows = matrix.rows();
        const size_t nbColumns = matrix.columns();

        for(size_t row = 0; row < nbRows; ++row) {
            for(size_t column = row + 1; column < nbColumns; ++column) {
                matrix(row, column) = 0;
            }
        }

        return matrix;
    }

    template<typename T>
    Matrix<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.");
        }

        Matrix<T> L = SimplifySymmetricMatrix(matrix);
        
        for(size_t k = 0; k < nbColumns; ++k) {
            const T & a_kk = L(k, k);
            
            if(a_kk > 0) { //(1)
                L(k, k) = std::sqrt(a_kk);

                for(size_t i = k + 1; i < nbRows; ++i) {
                    L(i, k) /= L(k, k); //(2)
                    
                    for(size_t j = k + 1; j <= i; ++j) { 
                        L(i, j) -= L(i, k) * L(j, k); //(3)
                    }
                }
            }
            else {
                throw std::domain_error("Matrix is not positive definit.");
            }
        }
        
        return L;
    }
}