25 #ifndef EIGEN_INCOMPLETE_LUT_H
26 #define EIGEN_INCOMPLETE_LUT_H
53 template <
typename _Scalar>
56 typedef _Scalar Scalar;
71 template<
typename MatrixType>
95 template<
typename MatrixType>
98 template<
typename MatrixType>
106 template<
typename MatrixType>
119 template<
typename Rhs,
typename Dest>
123 x =
m_lu.template triangularView<UnitLower>().
solve(x);
124 x =
m_lu.template triangularView<Upper>().
solve(x);
128 template<
typename Rhs>
inline const internal::solve_retval<IncompleteLUT, Rhs>
133 &&
"IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
134 return internal::solve_retval<IncompleteLUT, Rhs>(*
this, b.derived());
139 template <
typename VectorV,
typename VectorI>
168 template<
typename Scalar>
171 this->m_droptol = droptol;
178 template<
typename Scalar>
181 this->m_fillfactor = fillfactor;
194 template <
typename Scalar>
195 template <
typename VectorV,
typename VectorI>
206 if (ncut < first || ncut > last )
return 0;
211 for (
int j = first + 1; j <= last; j++) {
215 swap(ind(mid), ind(j));
219 swap(
row(mid),
row(first));
220 swap(ind(mid), ind(first));
222 if (mid > ncut) last = mid - 1;
223 else if (mid < ncut ) first = mid + 1;
224 }
while (mid != ncut );
229 template <
typename Scalar>
230 template<
typename _MatrixType>
241 internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P);
243 m_Pinv = m_P.inverse();
245 m_analysisIsOk =
true;
248 template <
typename Scalar>
249 template<
typename _MatrixType>
256 eigen_assert((amat.rows() == amat.cols()) &&
"The factorization should be done on a square matrix");
265 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
275 int fill_in =
static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1;
276 if (fill_in > n) fill_in = n;
279 int nnzL = fill_in/2;
281 m_lu.reserve(n * (nnzL + nnzU + 1));
284 for (
int ii = 0; ii < n; ii++)
293 RealScalar rownorm = 0;
295 typename FactorType::InnerIterator j_it(mat, ii);
298 int k = j_it.index();
303 u(sizel) = j_it.value();
309 u(ii) = j_it.value();
314 int jpos = ii + sizeu;
316 u(jpos) = j_it.value();
330 rownorm =
sqrt(rownorm);
340 int minrow = ju.segment(jj,sizel-jj).minCoeff(&k);
342 if (minrow != ju(jj))
347 jr(minrow) = jj; jr(j) = k;
354 typename FactorType::InnerIterator ki_it(m_lu, minrow);
355 while (ki_it && ki_it.index() < minrow) ++ki_it;
357 Scalar fact = u(jj) / ki_it.value();
360 if(
abs(fact) <= m_droptol)
368 for (; ki_it; ++ki_it)
370 Scalar prod = fact * ki_it.value();
371 int j = ki_it.index();
404 for(
int k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
410 len = (std::min)(sizel, nnzL);
411 typename Vector::SegmentReturnType ul(u.segment(0, sizel));
412 typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
413 QuickSplit(ul, jul, len);
417 for(
int k = 0; k < len; k++)
418 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
422 if (u(ii) == Scalar(0))
423 u(ii) =
sqrt(m_droptol) * rownorm;
424 m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
429 for(
int k = 1; k < sizeu; k++)
431 if(
abs(u(ii+k)) > m_droptol * rownorm )
434 u(ii + len) = u(ii + k);
435 ju(ii + len) = ju(ii + k);
439 len = (std::min)(sizeu, nnzU);
440 typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
441 typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
442 QuickSplit(uu, juu, len);
445 for(
int k = ii + 1; k < ii + len; k++)
446 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
450 m_lu.makeCompressed();
452 m_factorizationIsOk =
true;
458 template<
typename _MatrixType,
typename Rhs>
460 : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
465 template<typename Dest>
void evalTo(Dest& dst)
const
467 dec()._solve(rhs(),dst);
475 #endif // EIGEN_INCOMPLETE_LUT_H