Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.