25 #ifndef EIGEN_SKYLINEINPLACELU_H
26 #define EIGEN_SKYLINEINPLACELU_H
39 template<
typename MatrixType>
42 typedef typename MatrixType::Scalar Scalar;
43 typedef typename MatrixType::Index Index;
45 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
52 : m_flags(
flags), m_status(0), m_lu(matrix) {
53 m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar > ();
54 m_lu.IsRowMajor ? computeRowMajor() :
compute();
95 void setOrderingMethod(
int m) {
99 int orderingMethod()
const {
105 void computeRowMajor();
113 template<
typename BDerived,
typename XDerived>
114 bool solve(
const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x,
115 const int transposed = 0)
const;
123 RealScalar m_precision;
125 mutable int m_status;
133 template<
typename MatrixType>
136 const size_t rows = m_lu.rows();
137 const size_t cols = m_lu.cols();
139 eigen_assert(rows == cols &&
"We do not (yet) support rectangular LU.");
140 eigen_assert(!m_lu.IsRowMajor &&
"LU decomposition does not work with rowMajor Storage");
143 const double pivot = m_lu.coeffDiag(
row);
147 for (
typename MatrixType::InnerLowerIterator lIt(m_lu, col); lIt; ++lIt) {
148 lIt.valueRef() /= pivot;
152 typename MatrixType::InnerLowerIterator lIt(m_lu, col);
153 for (Index rrow =
row + 1; rrow < m_lu.rows(); rrow++) {
154 typename MatrixType::InnerUpperIterator uItPivot(m_lu,
row);
155 typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
156 const double coef = lIt.value();
158 uItPivot += (rrow -
row - 1);
161 for (++uItPivot; uIt && uItPivot;) {
162 uIt.valueRef() -= uItPivot.value() * coef;
171 typename MatrixType::InnerLowerIterator lIt3(m_lu, col);
172 for (Index rrow =
row + 1; rrow < m_lu.rows(); rrow++) {
173 typename MatrixType::InnerUpperIterator uItPivot(m_lu,
row);
174 const double coef = lIt3.value();
177 for (Index i = 0; i < rrow -
row - 1; i++) {
178 m_lu.coeffRefLower(rrow, row + i + 1) -= uItPivot.value() * coef;
184 typename MatrixType::InnerLowerIterator lIt2(m_lu, col);
185 for (Index rrow =
row + 1; rrow < m_lu.rows(); rrow++) {
187 typename MatrixType::InnerUpperIterator uItPivot(m_lu,
row);
188 typename MatrixType::InnerUpperIterator uIt(m_lu, rrow);
189 const double coef = lIt2.value();
191 uItPivot += (rrow -
row - 1);
192 m_lu.coeffRefDiag(rrow) -= uItPivot.value() * coef;
198 template<
typename MatrixType>
200 const size_t rows = m_lu.rows();
201 const size_t cols = m_lu.cols();
203 eigen_assert(rows == cols &&
"We do not (yet) support rectangular LU.");
204 eigen_assert(m_lu.IsRowMajor &&
"You're trying to apply rowMajor decomposition on a ColMajor matrix !");
206 for (Index row = 0;
row < rows;
row++) {
207 typename MatrixType::InnerLowerIterator llIt(m_lu, row);
210 for (Index col = llIt.col();
col <
row;
col++) {
211 if (m_lu.coeffExistLower(row, col)) {
212 const double diag = m_lu.coeffDiag(col);
214 typename MatrixType::InnerLowerIterator lIt(m_lu, row);
215 typename MatrixType::InnerUpperIterator uIt(m_lu, col);
218 const Index offset = lIt.col() - uIt.row();
221 Index stop = offset > 0 ?
col - lIt.col() :
col - uIt.row();
225 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
226 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
229 Scalar newCoeff = m_lu.coeffLower(row, col) - rowVal.dot(colVal);
235 Scalar newCoeff = m_lu.coeffLower(row, col);
237 for (Index k = 0; k < stop; ++k) {
238 const Scalar tmp = newCoeff;
239 newCoeff = tmp - lIt.value() * uIt.value();
245 m_lu.coeffRefLower(row, col) = newCoeff / diag;
250 const Index
col = row;
251 typename MatrixType::InnerUpperIterator uuIt(m_lu, col);
252 for (Index rrow = uuIt.row(); rrow < col; rrow++) {
254 typename MatrixType::InnerLowerIterator lIt(m_lu, rrow);
255 typename MatrixType::InnerUpperIterator uIt(m_lu, col);
256 const Index offset = lIt.col() - uIt.row();
258 Index stop = offset > 0 ? rrow - lIt.col() : rrow - uIt.row();
261 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
262 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
264 Scalar newCoeff = m_lu.coeffUpper(rrow, col) - rowVal.dot(colVal);
270 Scalar newCoeff = m_lu.coeffUpper(rrow, col);
271 for (Index k = 0; k < stop; ++k) {
272 const Scalar tmp = newCoeff;
273 newCoeff = tmp - lIt.value() * uIt.value();
279 m_lu.coeffRefUpper(rrow, col) = newCoeff;
284 typename MatrixType::InnerLowerIterator lIt(m_lu, row);
285 typename MatrixType::InnerUpperIterator uIt(m_lu, row);
287 const Index offset = lIt.col() - uIt.row();
290 Index stop = offset > 0 ? lIt.size() : uIt.size();
292 Map<VectorXd > rowVal(lIt.valuePtr() + (offset > 0 ? 0 : -offset), stop);
293 Map<VectorXd > colVal(uIt.valuePtr() + (offset > 0 ? offset : 0), stop);
294 Scalar newCoeff = m_lu.coeffDiag(row) - rowVal.dot(colVal);
300 Scalar newCoeff = m_lu.coeffDiag(row);
301 for (Index k = 0; k < stop; ++k) {
302 const Scalar tmp = newCoeff;
303 newCoeff = tmp - lIt.value() * uIt.value();
308 m_lu.coeffRefDiag(row) = newCoeff;
320 template<
typename MatrixType>
321 template<
typename BDerived,
typename XDerived>
323 const size_t rows = m_lu.rows();
324 const size_t cols = m_lu.cols();
327 for (Index row = 0; row < rows; row++) {
328 x->coeffRef(row) = b.coeff(row);
329 Scalar newVal = x->coeff(row);
330 typename MatrixType::InnerLowerIterator lIt(m_lu, row);
332 Index col = lIt.col();
333 while (lIt.col() < row) {
335 newVal -= x->coeff(col++) * lIt.value();
339 x->coeffRef(row) = newVal;
343 for (Index col = rows - 1; col > 0; col--) {
344 x->coeffRef(col) = x->coeff(col) / m_lu.coeffDiag(col);
346 const Scalar x_col = x->coeff(col);
348 typename MatrixType::InnerUpperIterator uIt(m_lu, col);
353 x->coeffRef(uIt.row()) -= x_col * uIt.value();
360 x->coeffRef(0) = x->coeff(0) / m_lu.coeffDiag(0);
367 #endif // EIGEN_SKYLINELU_H