10 #ifndef EIGEN_SKYLINEINPLACELU_H
11 #define EIGEN_SKYLINEINPLACELU_H
24 template<
typename MatrixType>
98 template<
typename BDerived,
typename XDerived>
100 const int transposed = 0)
const;
118 template<
typename MatrixType>
121 const size_t rows = m_lu.rows();
122 const size_t cols = m_lu.cols();
125 eigen_assert(!m_lu.IsRowMajor &&
"LU decomposition does not work with rowMajor Storage");
128 const double pivot = m_lu.coeffDiag(
row);
132 for (
typename MatrixType::InnerLowerIterator lIt(m_lu,
col); lIt; ++lIt) {
133 lIt.valueRef() /= pivot;
137 typename MatrixType::InnerLowerIterator lIt(m_lu,
col);
138 for (
Index rrow =
row + 1; rrow < m_lu.rows(); rrow++) {
139 typename MatrixType::InnerUpperIterator uItPivot(m_lu,
row);
140 typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
141 const double coef = lIt.value();
143 uItPivot += (rrow -
row - 1);
146 for (++uItPivot; uIt && uItPivot;) {
147 uIt.valueRef() -= uItPivot.value() * coef;
156 typename MatrixType::InnerLowerIterator lIt3(m_lu,
col);
157 for (
Index rrow =
row + 1; rrow < m_lu.rows(); rrow++) {
158 typename MatrixType::InnerUpperIterator uItPivot(m_lu,
row);
159 const double coef = lIt3.value();
163 m_lu.coeffRefLower(rrow,
row +
i + 1) -= uItPivot.value() * coef;
169 typename MatrixType::InnerLowerIterator lIt2(m_lu,
col);
170 for (
Index rrow =
row + 1; rrow < m_lu.rows(); rrow++) {
172 typename MatrixType::InnerUpperIterator uItPivot(m_lu,
row);
173 typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
174 const double coef = lIt2.value();
176 uItPivot += (rrow -
row - 1);
177 m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
183 template<
typename MatrixType>
185 const size_t rows = m_lu.rows();
186 const size_t cols = m_lu.cols();
189 eigen_assert(m_lu.IsRowMajor &&
"You're trying to apply rowMajor decomposition on a ColMajor matrix !");
192 typename MatrixType::InnerLowerIterator llIt(m_lu,
row);
196 if (m_lu.coeffExistLower(
row,
col)) {
197 const double diag = m_lu.coeffDiag(
col);
199 typename MatrixType::InnerLowerIterator lIt(m_lu,
row);
200 typename MatrixType::InnerUpperIterator uIt(m_lu,
col);
214 Scalar newCoeff = m_lu.coeffLower(
row,
col) - rowVal.dot(colVal);
222 for (
Index k = 0; k < stop; ++k) {
223 const Scalar tmp = newCoeff;
224 newCoeff = tmp - lIt.value() * uIt.value();
230 m_lu.coeffRefLower(
row,
col) = newCoeff /
diag;
236 typename MatrixType::InnerUpperIterator uuIt(m_lu,
col);
237 for (
Index rrow = uuIt.row(); rrow <
col; rrow++) {
239 typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
240 typename MatrixType::InnerUpperIterator uIt(m_lu,
col);
243 Index stop =
offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
249 Scalar newCoeff = m_lu.coeffUpper(rrow,
col) - rowVal.dot(colVal);
255 Scalar newCoeff = m_lu.coeffUpper(rrow,
col);
256 for (
Index k = 0; k < stop; ++k) {
257 const Scalar tmp = newCoeff;
258 newCoeff = tmp - lIt.value() * uIt.value();
264 m_lu.coeffRefUpper(rrow,
col) = newCoeff;
269 typename MatrixType::InnerLowerIterator lIt(m_lu,
row);
270 typename MatrixType::InnerUpperIterator uIt(m_lu,
row);
279 Scalar newCoeff = m_lu.coeffDiag(
row) - rowVal.dot(colVal);
286 for (
Index k = 0; k < stop; ++k) {
287 const Scalar tmp = newCoeff;
288 newCoeff = tmp - lIt.value() * uIt.value();
293 m_lu.coeffRefDiag(
row) = newCoeff;
305 template<
typename MatrixType>
306 template<
typename BDerived,
typename XDerived>
308 const size_t rows = m_lu.rows();
309 const size_t cols = m_lu.cols();
315 typename MatrixType::InnerLowerIterator lIt(m_lu,
row);
318 while (lIt.col() <
row) {
320 newVal -=
x->coeff(
col++) * lIt.value();
324 x->coeffRef(
row) = newVal;
329 x->coeffRef(
col) =
x->coeff(
col) / m_lu.coeffDiag(
col);
333 typename MatrixType::InnerUpperIterator uIt(m_lu,
col);
338 x->coeffRef(uIt.row()) -= x_col * uIt.value();
345 x->coeffRef(0) =
x->coeff(0) / m_lu.coeffDiag(0);
352 #endif // EIGEN_SKYLINEINPLACELU_H