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