Skip to content

Instantly share code, notes, and snippets.

@Ben1980

Ben1980/luDecomposition.h

Last active Nov 14, 2019
Embed
What would you like to do?
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);
for(size_t column = 0; column < nbColumns; ++column) {
decomposition.L(column, column) = 1; //(1)
for(size_t row = column + 1; row < nbRows; ++row) {
const T & divisor = decomposition.U(column, column);
if(std::fabs(divisor) < std::numeric_limits<T>::min()) {
throw std::domain_error("Division by 0."); //a_ii != 0 is necessary because of pivoting with diaognal strategy
}
decomposition.L(row, column) = decomposition.U(row, column) / divisor; //(2)
for(size_t col = column; col < nbColumns; ++col) {
decomposition.U(row, col) -= decomposition.L(row, column) * decomposition.U(column, col); //(3)
}
}
}
return decomposition;
}
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.