10 #ifndef EIGEN_SKYLINEINPLACELU_H 11 #define EIGEN_SKYLINEINPLACELU_H 24 template<
typename MatrixType>
27 typedef typename MatrixType::Scalar
Scalar;
28 typedef typename MatrixType::Index
Index;
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();
124 eigen_assert(rows == cols &&
"We do not (yet) support rectangular LU.");
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();
162 for (
Index i = 0; i < rrow -
row - 1; i++) {
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();
188 eigen_assert(rows == cols &&
"We do not (yet) support rectangular LU.");
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);
203 const Index offset = lIt.col() - uIt.row();
206 Index stop = offset > 0 ?
col - lIt.col() :
col - uIt.row();
210 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
211 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
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);
241 const Index offset = lIt.col() - uIt.row();
243 Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
246 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
247 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
249 Scalar newCoeff =
m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
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);
272 const Index offset = lIt.col() - uIt.row();
275 Index stop = offset > 0 ? lIt.size() : uIt.size();
277 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
278 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
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();
313 x->coeffRef(
row) = b.coeff(
row);
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;
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_SKYLINELU_H
A matrix or vector expression mapping an existing array of data.
void setOrderingMethod(int m)
NumTraits< typename MatrixType::Scalar >::Real RealScalar
SkylineInplaceLU(MatrixType &matrix, int flags=0)
Inplace LU decomposition of a skyline matrix and associated features.
iterative scaling algorithm to equilibrate rows and column norms in matrices
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
bool succeeded(void) const
bool solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > *x, const int transposed=0) const
CoeffReturnType value() const
RealScalar precision() const
MatrixType::Scalar Scalar
int orderingMethod() const
Base class for all dense matrices, vectors, and expressions.
void setPrecision(RealScalar v)