Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
namespace PivotLUDecomposition {
template<typename T>
struct Decomposition {
Matrix<T> P;
Matrix<T> L;
Matrix<T> U;
Decomposition(const Matrix<T> &matrix) : P(matrix.rows(), matrix.columns()), L(matrix.rows(), matrix.columns()), U(matrix) {}
};
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);
decomposition.P = MatrixFactory::IdentityMatrix<T>(nbRows);
for(size_t k = 0; k < nbRows; ++k) {
T max = 0;
size_t pk = 0;
for(size_t i = k; i < nbRows; ++i) {
T s = 0;
for(size_t j = k; j < nbColumns; ++j) {
s += std::fabs(decomposition.U(i, j)); //(1)
}
T q = std::fabs(decomposition.U(i, k)) / s; //(2)
if(q > max) {
max = q;
pk = i;
}
}
if(std::fabs(max) < std::numeric_limits<T>::min()) {
throw std::domain_error("Pivot has 0 value.");
}
if(pk != k) {
for(size_t j = 0; j < nbColumns; ++j) { //(3)
std::swap(decomposition.P(k, j), decomposition.P(pk, j));
std::swap(decomposition.L(k, j), decomposition.L(pk, j));
std::swap(decomposition.U(k, j), decomposition.U(pk, j));
}
}
for(size_t i = k+1; i < nbRows; ++i) {
decomposition.L(i, k) = decomposition.U(i, k)/decomposition.U(k, k); //(4)
for(size_t j = k; j < nbColumns; ++j) {
decomposition.U(i, j) = decomposition.U(i, j) - decomposition.L(i, k) * decomposition.U(k, j); //(5)
}
}
}
for(size_t k = 0; k < nbRows; ++k) {
decomposition.L(k,k) = 1;
}
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.