11 #ifndef EIGEN_INCOMPLETE_LUT_H
12 #define EIGEN_INCOMPLETE_LUT_H
28 template <
typename VectorV,
typename VectorI>
29 Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
31 typedef typename VectorV::RealScalar RealScalar;
41 if (ncut < first || ncut > last )
return 0;
45 RealScalar abskey = abs(row(mid));
46 for (Index j = first + 1; j <= last; j++) {
47 if ( abs(row(j)) > abskey) {
49 swap(row(mid), row(j));
50 swap(ind(mid), ind(j));
54 swap(row(mid), row(first));
55 swap(ind(mid), ind(first));
57 if (mid > ncut) last = mid - 1;
58 else if (mid < ncut ) first = mid + 1;
59 }
while (mid != ncut );
96 template <
typename _Scalar,
typename _StorageIndex =
int>
101 using Base::m_isInitialized;
103 typedef _Scalar Scalar;
104 typedef _StorageIndex StorageIndex;
117 m_analysisIsOk(
false), m_factorizationIsOk(
false)
120 template<
typename MatrixType>
122 : m_droptol(droptol),m_fillfactor(fillfactor),
123 m_analysisIsOk(
false),m_factorizationIsOk(
false)
125 eigen_assert(fillfactor != 0);
129 Index rows()
const {
return m_lu.
rows(); }
131 Index cols()
const {
return m_lu.
cols(); }
140 eigen_assert(m_isInitialized &&
"IncompleteLUT is not initialized.");
144 template<
typename MatrixType>
145 void analyzePattern(
const MatrixType& amat);
147 template<
typename MatrixType>
148 void factorize(
const MatrixType& amat);
155 template<
typename MatrixType>
158 analyzePattern(amat);
166 template<
typename Rhs,
typename Dest>
167 void _solve_impl(
const Rhs& b, Dest& x)
const
170 x = m_lu.template triangularView<UnitLower>().
solve(x);
171 x = m_lu.template triangularView<Upper>().
solve(x);
179 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
188 RealScalar m_droptol;
191 bool m_factorizationIsOk;
201 template<
typename Scalar,
typename StorageIndex>
204 this->m_droptol = droptol;
211 template<
typename Scalar,
typename StorageIndex>
214 this->m_fillfactor = fillfactor;
217 template <
typename Scalar,
typename StorageIndex>
218 template<
typename _MatrixType>
228 AtA.
prune(keep_diag());
229 internal::minimum_degree_ordering<Scalar, StorageIndex>(AtA, m_P);
231 m_Pinv = m_P.inverse();
233 m_analysisIsOk =
true;
234 m_factorizationIsOk =
false;
235 m_isInitialized =
true;
238 template <
typename Scalar,
typename StorageIndex>
239 template<
typename _MatrixType>
240 void IncompleteLUT<Scalar,StorageIndex>::factorize(
const _MatrixType& amat)
245 using internal::convert_index;
247 eigen_assert((amat.rows() == amat.cols()) &&
"The factorization should be done on a square matrix");
248 Index n = amat.cols();
256 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
257 SparseMatrix<Scalar,RowMajor, StorageIndex> mat;
258 mat = amat.twistedBy(m_Pinv);
266 Index fill_in = (amat.nonZeros()*m_fillfactor)/n + 1;
267 if (fill_in > n) fill_in = n;
270 Index nnzL = fill_in/2;
272 m_lu.reserve(n * (nnzL + nnzU + 1));
275 for (Index ii = 0; ii < n; ii++)
281 ju(ii) = convert_index<StorageIndex>(ii);
283 jr(ii) = convert_index<StorageIndex>(ii);
284 RealScalar rownorm = 0;
286 typename FactorType::InnerIterator j_it(mat, ii);
289 Index k = j_it.index();
293 ju(sizel) = convert_index<StorageIndex>(k);
294 u(sizel) = j_it.value();
295 jr(k) = convert_index<StorageIndex>(sizel);
300 u(ii) = j_it.value();
305 Index jpos = ii + sizeu;
306 ju(jpos) = convert_index<StorageIndex>(k);
307 u(jpos) = j_it.value();
308 jr(k) = convert_index<StorageIndex>(jpos);
311 rownorm += numext::abs2(j_it.value());
321 rownorm = sqrt(rownorm);
331 Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k);
333 if (minrow != ju(jj))
338 jr(minrow) = convert_index<StorageIndex>(jj);
339 jr(j) = convert_index<StorageIndex>(k);
346 typename FactorType::InnerIterator ki_it(m_lu, minrow);
347 while (ki_it && ki_it.index() < minrow) ++ki_it;
348 eigen_internal_assert(ki_it && ki_it.col()==minrow);
349 Scalar fact = u(jj) / ki_it.value();
352 if(abs(fact) <= m_droptol)
360 for (; ki_it; ++ki_it)
362 Scalar prod = fact * ki_it.value();
363 Index j = ki_it.index();
372 eigen_internal_assert(sizeu<=n);
378 eigen_internal_assert(sizel<=ii);
380 ju(newpos) = convert_index<StorageIndex>(j);
382 jr(j) = convert_index<StorageIndex>(newpos);
389 ju(len) = convert_index<StorageIndex>(minrow);
396 for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
402 len = (std::min)(sizel, nnzL);
403 typename Vector::SegmentReturnType ul(u.segment(0, sizel));
404 typename VectorI::SegmentReturnType jul(ju.segment(0, sizel));
405 internal::QuickSplit(ul, jul, len);
409 for(Index k = 0; k < len; k++)
410 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
414 if (u(ii) == Scalar(0))
415 u(ii) = sqrt(m_droptol) * rownorm;
416 m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
421 for(Index k = 1; k < sizeu; k++)
423 if(abs(u(ii+k)) > m_droptol * rownorm )
426 u(ii + len) = u(ii + k);
427 ju(ii + len) = ju(ii + k);
431 len = (std::min)(sizeu, nnzU);
432 typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
433 typename VectorI::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
434 internal::QuickSplit(uu, juu, len);
437 for(Index k = ii + 1; k < ii + len; k++)
438 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
441 m_lu.makeCompressed();
443 m_factorizationIsOk =
true;
449 #endif // EIGEN_INCOMPLETE_LUT_H
void setDroptol(const RealScalar &droptol)
Definition: IncompleteLUT.h:202
Index cols() const
Definition: SparseMatrix.h:132
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteLUT.h:138
const Solve< IncompleteLUT< _Scalar, _StorageIndex >, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SparseSolverBase.h:74
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Definition: IncompleteLUT.h:178
void prune(const Scalar &reference, const RealScalar &epsilon=NumTraits< RealScalar >::dummy_precision())
Definition: SparseMatrix.h:495
Definition: Constants.h:426
Incomplete LU factorization with dual-threshold strategy.
Definition: IncompleteLUT.h:97
Definition: Constants.h:424
Definition: Eigen_Colamd.h:54
IncompleteLUT & compute(const MatrixType &amat)
Definition: IncompleteLUT.h:156
ComputationInfo
Definition: Constants.h:422
Index rows() const
Definition: SparseMatrix.h:130
void setFillfactor(int fillfactor)
Definition: IncompleteLUT.h:212