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);
153 m_isInitialized = m_factorizationIsOk;
160 template<
typename Rhs,
typename Dest>
161 void _solve(
const Rhs& b, Dest& x)
const
164 x = m_lu.template triangularView<UnitLower>().solve(x);
165 x = m_lu.template triangularView<Upper>().solve(x);
169 template<
typename Rhs>
inline const internal::solve_retval<IncompleteLUT, Rhs>
170 solve(
const MatrixBase<Rhs>& b)
const
172 eigen_assert(m_isInitialized &&
"IncompleteLUT is not initialized.");
173 eigen_assert(cols()==b.rows()
174 &&
"IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
175 return internal::solve_retval<IncompleteLUT, Rhs>(*
this, b.derived());
182 inline bool operator() (
const Index& row,
const Index& col,
const Scalar&)
const
191 RealScalar m_droptol;
194 bool m_factorizationIsOk;
195 bool m_isInitialized;
205 template<
typename Scalar>
208 this->m_droptol = droptol;
215 template<
typename Scalar>
218 this->m_fillfactor = fillfactor;
221 template <
typename Scalar>
222 template<
typename _MatrixType>
232 AtA.
prune(keep_diag());
233 internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P);
235 m_Pinv = m_P.inverse();
237 m_analysisIsOk =
true;
240 template <
typename Scalar>
241 template<
typename _MatrixType>
242 void IncompleteLUT<Scalar>::factorize(
const _MatrixType& amat)
248 eigen_assert((amat.rows() == amat.cols()) &&
"The factorization should be done on a square matrix");
249 Index n = amat.cols();
257 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
258 SparseMatrix<Scalar,RowMajor, Index> mat;
259 mat = amat.twistedBy(m_Pinv);
267 Index fill_in =
static_cast<Index
> (amat.nonZeros()*m_fillfactor)/n+1;
268 if (fill_in > n) fill_in = n;
271 Index nnzL = fill_in/2;
273 m_lu.reserve(n * (nnzL + nnzU + 1));
276 for (Index ii = 0; ii < n; ii++)
285 RealScalar rownorm = 0;
287 typename FactorType::InnerIterator j_it(mat, ii);
290 Index k = j_it.index();
295 u(sizel) = j_it.value();
301 u(ii) = j_it.value();
306 Index jpos = ii + sizeu;
308 u(jpos) = j_it.value();
312 rownorm += numext::abs2(j_it.value());
322 rownorm = sqrt(rownorm);
332 Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k);
334 if (minrow != ju(jj))
339 jr(minrow) = jj; jr(j) = 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);
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 VectorXi::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 VectorXi::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);
442 m_lu.makeCompressed();
444 m_factorizationIsOk =
true;
450 template<
typename _MatrixType,
typename Rhs>
451 struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
452 : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
454 typedef IncompleteLUT<_MatrixType> Dec;
455 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
457 template<typename Dest>
void evalTo(Dest& dst)
const
459 dec()._solve(rhs(),dst);
467 #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:181
void setFillfactor(int fillfactor)
Definition: IncompleteLUT.h:216
Incomplete LU factorization with dual-threshold strategy.
Definition: IncompleteLUT.h:96
IncompleteLUT< Scalar > & compute(const MatrixType &amat)
Definition: IncompleteLUT.h:149
Definition: Eigen_Colamd.h:54
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:206
Matrix< int, Dynamic, 1 > VectorXi
Definition: Matrix.h:393
ComputationInfo
Definition: Constants.h:374