10 #ifndef EIGEN_INCOMPLETE_LUT_H
11 #define EIGEN_INCOMPLETE_LUT_H
27 template <
typename VectorV,
typename VectorI,
typename Index>
28 Index QuickSplit(VectorV &row, VectorI &ind, Index ncut)
30 typedef typename VectorV::RealScalar RealScalar;
40 if (ncut < first || ncut > last )
return 0;
44 RealScalar abskey = abs(row(mid));
45 for (Index j = first + 1; j <= last; j++) {
46 if ( abs(row(j)) > abskey) {
48 swap(row(mid), row(j));
49 swap(ind(mid), ind(j));
53 swap(row(mid), row(first));
54 swap(ind(mid), ind(first));
56 if (mid > ncut) last = mid - 1;
57 else if (mid < ncut ) first = mid + 1;
58 }
while (mid != ncut );
95 template <
typename _Scalar>
98 typedef _Scalar Scalar;
103 typedef typename FactorType::Index Index;
110 m_analysisIsOk(
false), m_factorizationIsOk(
false), m_isInitialized(
false)
113 template<
typename MatrixType>
115 : m_droptol(droptol),m_fillfactor(fillfactor),
116 m_analysisIsOk(
false),m_factorizationIsOk(
false),m_isInitialized(
false)
118 eigen_assert(fillfactor != 0);
122 Index rows()
const {
return m_lu.
rows(); }
124 Index cols()
const {
return m_lu.
cols(); }
133 eigen_assert(m_isInitialized &&
"IncompleteLUT is not initialized.");
137 template<
typename MatrixType>
138 void analyzePattern(
const MatrixType& amat);
140 template<
typename MatrixType>
141 void factorize(
const MatrixType& amat);
148 template<
typename MatrixType>
151 analyzePattern(amat);
159 template<
typename Rhs,
typename Dest>
160 void _solve(
const Rhs& b, Dest& x)
const
163 x = m_lu.template triangularView<UnitLower>().solve(x);
164 x = m_lu.template triangularView<Upper>().solve(x);
168 template<
typename Rhs>
inline const internal::solve_retval<IncompleteLUT, Rhs>
169 solve(
const MatrixBase<Rhs>& b)
const
171 eigen_assert(m_isInitialized &&
"IncompleteLUT is not initialized.");
172 eigen_assert(cols()==b.rows()
173 &&
"IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
174 return internal::solve_retval<IncompleteLUT, Rhs>(*
this, b.derived());
181 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
190 RealScalar m_droptol;
193 bool m_factorizationIsOk;
194 bool m_isInitialized;
204 template<
typename Scalar>
207 this->m_droptol = droptol;
214 template<
typename Scalar>
217 this->m_fillfactor = fillfactor;
220 template <
typename Scalar>
221 template<
typename _MatrixType>
231 AtA.
prune(keep_diag());
232 internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P);
234 m_Pinv = m_P.inverse();
236 m_analysisIsOk =
true;
237 m_factorizationIsOk =
false;
238 m_isInitialized =
false;
241 template <
typename Scalar>
242 template<
typename _MatrixType>
243 void IncompleteLUT<Scalar>::factorize(
const _MatrixType& amat)
249 eigen_assert((amat.rows() == amat.cols()) &&
"The factorization should be done on a square matrix");
250 Index n = amat.cols();
258 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
259 SparseMatrix<Scalar,RowMajor, Index> mat;
260 mat = amat.twistedBy(m_Pinv);
268 Index fill_in =
static_cast<Index
> (amat.nonZeros()*m_fillfactor)/n+1;
269 if (fill_in > n) fill_in = n;
272 Index nnzL = fill_in/2;
274 m_lu.reserve(n * (nnzL + nnzU + 1));
277 for (Index ii = 0; ii < n; ii++)
286 RealScalar rownorm = 0;
288 typename FactorType::InnerIterator j_it(mat, ii);
291 Index k = j_it.index();
296 u(sizel) = j_it.value();
302 u(ii) = j_it.value();
307 Index jpos = ii + sizeu;
309 u(jpos) = j_it.value();
313 rownorm += numext::abs2(j_it.value());
323 rownorm = sqrt(rownorm);
333 Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k);
335 if (minrow != ju(jj))
340 jr(minrow) = jj; jr(j) = k;
347 typename FactorType::InnerIterator ki_it(m_lu, minrow);
348 while (ki_it && ki_it.index() < minrow) ++ki_it;
349 eigen_internal_assert(ki_it && ki_it.col()==minrow);
350 Scalar fact = u(jj) / ki_it.value();
353 if(abs(fact) <= m_droptol)
361 for (; ki_it; ++ki_it)
363 Scalar prod = fact * ki_it.value();
364 Index j = ki_it.index();
373 eigen_internal_assert(sizeu<=n);
379 eigen_internal_assert(sizel<=ii);
397 for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
403 len = (std::min)(sizel, nnzL);
404 typename Vector::SegmentReturnType ul(u.segment(0, sizel));
405 typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
406 internal::QuickSplit(ul, jul, len);
410 for(Index k = 0; k < len; k++)
411 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
415 if (u(ii) == Scalar(0))
416 u(ii) = sqrt(m_droptol) * rownorm;
417 m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
422 for(Index k = 1; k < sizeu; k++)
424 if(abs(u(ii+k)) > m_droptol * rownorm )
427 u(ii + len) = u(ii + k);
428 ju(ii + len) = ju(ii + k);
432 len = (std::min)(sizeu, nnzU);
433 typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
434 typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
435 internal::QuickSplit(uu, juu, len);
438 for(Index k = ii + 1; k < ii + len; k++)
439 m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
443 m_lu.makeCompressed();
445 m_factorizationIsOk =
true;
446 m_isInitialized = m_factorizationIsOk;
452 template<
typename _MatrixType,
typename Rhs>
453 struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
454 : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
456 typedef IncompleteLUT<_MatrixType> Dec;
457 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
459 template<typename Dest>
void evalTo(Dest& dst)
const
461 dec()._solve(rhs(),dst);
469 #endif // EIGEN_INCOMPLETE_LUT_H
Index rows() const
Definition: SparseMatrix.h:119
Index cols() const
Definition: SparseMatrix.h:121
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: IncompleteLUT.h:131
void fill(const Scalar &value)
Definition: CwiseNullaryOp.h:322
Definition: Constants.h:378
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Definition: IncompleteLUT.h:180
void setFillfactor(int fillfactor)
Definition: IncompleteLUT.h:215
Incomplete LU factorization with dual-threshold strategy.
Definition: IncompleteLUT.h:96
IncompleteLUT< Scalar > & compute(const MatrixType &amat)
Definition: IncompleteLUT.h:149
void prune(const Scalar &reference, const RealScalar &epsilon=NumTraits< RealScalar >::dummy_precision())
Definition: SparseMatrix.h:490
Definition: Constants.h:376
void setDroptol(const RealScalar &droptol)
Definition: IncompleteLUT.h:205
Matrix< int, Dynamic, 1 > VectorXi
Definition: Matrix.h:393
ComputationInfo
Definition: Constants.h:374