Skip to content
{{ message }}

Instantly share code, notes, and snippets.

Ben1980/choleskyDecomposition.h

Last active Nov 27, 2019
 namespace CholeskyDecomposition { template Matrix SimplifySymmetricMatrix(Matrix 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 Matrix Decompose(const Matrix &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 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; } }
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.