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