diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers/IncompleteLUT.h')
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 215 |
1 files changed, 108 insertions, 107 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 224304f0e..b55afc136 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -10,36 +10,88 @@ #ifndef EIGEN_INCOMPLETE_LUT_H #define EIGEN_INCOMPLETE_LUT_H + namespace Eigen { -/** - * \brief Incomplete LU factorization with dual-threshold strategy - * During the numerical factorization, two dropping rules are used : - * 1) any element whose magnitude is less than some tolerance is dropped. - * This tolerance is obtained by multiplying the input tolerance @p droptol - * by the average magnitude of all the original elements in the current row. - * 2) After the elimination of the row, only the @p fill largest elements in - * the L part and the @p fill largest elements in the U part are kept - * (in addition to the diagonal element ). Note that @p fill is computed from - * the input parameter @p fillfactor which is used the ratio to control the fill_in - * relatively to the initial number of nonzero elements. - * - * The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements) - * and when @p fill=n/2 with @p droptol being different to zero. - * - * References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization, - * Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994. - * - * NOTE : The following implementation is derived from the ILUT implementation - * in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota - * released under the terms of the GNU LGPL: - * http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README - * However, Yousef Saad gave us permission to relicense his ILUT code to MPL2. - * See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012: - * http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html - * alternatively, on GMANE: - * http://comments.gmane.org/gmane.comp.lib.eigen/3302 - */ +namespace internal { + +/** \internal + * Compute a quick-sort split of a vector + * On output, the vector row is permuted such that its elements satisfy + * abs(row(i)) >= abs(row(ncut)) if i<ncut + * abs(row(i)) <= abs(row(ncut)) if i>ncut + * \param row The vector of values + * \param ind The array of index for the elements in @p row + * \param ncut The number of largest elements to keep + **/ +template <typename VectorV, typename VectorI, typename Index> +Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) +{ + typedef typename VectorV::RealScalar RealScalar; + using std::swap; + using std::abs; + Index mid; + Index n = row.size(); /* length of the vector */ + Index first, last ; + + ncut--; /* to fit the zero-based indices */ + first = 0; + last = n-1; + if (ncut < first || ncut > last ) return 0; + + do { + mid = first; + RealScalar abskey = abs(row(mid)); + for (Index j = first + 1; j <= last; j++) { + if ( abs(row(j)) > abskey) { + ++mid; + swap(row(mid), row(j)); + swap(ind(mid), ind(j)); + } + } + /* Interchange for the pivot element */ + swap(row(mid), row(first)); + swap(ind(mid), ind(first)); + + if (mid > ncut) last = mid - 1; + else if (mid < ncut ) first = mid + 1; + } while (mid != ncut ); + + return 0; /* mid is equal to ncut */ +} + +}// end namespace internal + +/** \ingroup IterativeLinearSolvers_Module + * \class IncompleteLUT + * \brief Incomplete LU factorization with dual-threshold strategy + * + * During the numerical factorization, two dropping rules are used : + * 1) any element whose magnitude is less than some tolerance is dropped. + * This tolerance is obtained by multiplying the input tolerance @p droptol + * by the average magnitude of all the original elements in the current row. + * 2) After the elimination of the row, only the @p fill largest elements in + * the L part and the @p fill largest elements in the U part are kept + * (in addition to the diagonal element ). Note that @p fill is computed from + * the input parameter @p fillfactor which is used the ratio to control the fill_in + * relatively to the initial number of nonzero elements. + * + * The two extreme cases are when @p droptol=0 (to keep all the @p fill*2 largest elements) + * and when @p fill=n/2 with @p droptol being different to zero. + * + * References : Yousef Saad, ILUT: A dual threshold incomplete LU factorization, + * Numerical Linear Algebra with Applications, 1(4), pp 387-402, 1994. + * + * NOTE : The following implementation is derived from the ILUT implementation + * in the SPARSKIT package, Copyright (C) 2005, the Regents of the University of Minnesota + * released under the terms of the GNU LGPL: + * http://www-users.cs.umn.edu/~saad/software/SPARSKIT/README + * However, Yousef Saad gave us permission to relicense his ILUT code to MPL2. + * See the Eigen mailing list archive, thread: ILUT, date: July 8, 2012: + * http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/07/msg00064.html + * alternatively, on GMANE: + * http://comments.gmane.org/gmane.comp.lib.eigen/3302 + */ template <typename _Scalar> class IncompleteLUT : internal::noncopyable { @@ -59,7 +111,7 @@ class IncompleteLUT : internal::noncopyable {} template<typename MatrixType> - IncompleteLUT(const MatrixType& mat, RealScalar droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10) + IncompleteLUT(const MatrixType& mat, const RealScalar& droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10) : m_droptol(droptol),m_fillfactor(fillfactor), m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false) { @@ -98,12 +150,11 @@ class IncompleteLUT : internal::noncopyable { analyzePattern(amat); factorize(amat); - eigen_assert(m_factorizationIsOk == true); - m_isInitialized = true; + m_isInitialized = m_factorizationIsOk; return *this; } - void setDroptol(RealScalar droptol); + void setDroptol(const RealScalar& droptol); void setFillfactor(int fillfactor); template<typename Rhs, typename Dest> @@ -126,10 +177,6 @@ class IncompleteLUT : internal::noncopyable protected: - template <typename VectorV, typename VectorI> - int QuickSplit(VectorV &row, VectorI &ind, int ncut); - - /** keeps off-diagonal entries; drops diagonal entries */ struct keep_diag { inline bool operator() (const Index& row, const Index& col, const Scalar&) const @@ -156,7 +203,7 @@ protected: * \param droptol Drop any element whose magnitude is less than this tolerance **/ template<typename Scalar> -void IncompleteLUT<Scalar>::setDroptol(RealScalar droptol) +void IncompleteLUT<Scalar>::setDroptol(const RealScalar& droptol) { this->m_droptol = droptol; } @@ -171,51 +218,6 @@ void IncompleteLUT<Scalar>::setFillfactor(int fillfactor) this->m_fillfactor = fillfactor; } - -/** - * Compute a quick-sort split of a vector - * On output, the vector row is permuted such that its elements satisfy - * abs(row(i)) >= abs(row(ncut)) if i<ncut - * abs(row(i)) <= abs(row(ncut)) if i>ncut - * \param row The vector of values - * \param ind The array of index for the elements in @p row - * \param ncut The number of largest elements to keep - **/ -template <typename Scalar> -template <typename VectorV, typename VectorI> -int IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut) -{ - using std::swap; - int mid; - int n = row.size(); /* length of the vector */ - int first, last ; - - ncut--; /* to fit the zero-based indices */ - first = 0; - last = n-1; - if (ncut < first || ncut > last ) return 0; - - do { - mid = first; - RealScalar abskey = std::abs(row(mid)); - for (int j = first + 1; j <= last; j++) { - if ( std::abs(row(j)) > abskey) { - ++mid; - swap(row(mid), row(j)); - swap(ind(mid), ind(j)); - } - } - /* Interchange for the pivot element */ - swap(row(mid), row(first)); - swap(ind(mid), ind(first)); - - if (mid > ncut) last = mid - 1; - else if (mid < ncut ) first = mid + 1; - } while (mid != ncut ); - - return 0; /* mid is equal to ncut */ -} - template <typename Scalar> template<typename _MatrixType> void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat) @@ -244,7 +246,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) using std::abs; eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix"); - int n = amat.cols(); // Size of the matrix + Index n = amat.cols(); // Size of the matrix m_lu.resize(n,n); // Declare Working vectors and variables Vector u(n) ; // real values of the row -- maximum size is n -- @@ -262,21 +264,21 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) u.fill(0); // number of largest elements to keep in each row: - int fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1; + Index fill_in = static_cast<Index> (amat.nonZeros()*m_fillfactor)/n+1; if (fill_in > n) fill_in = n; // number of largest nonzero elements to keep in the L and the U part of the current row: - int nnzL = fill_in/2; - int nnzU = nnzL; + Index nnzL = fill_in/2; + Index nnzU = nnzL; m_lu.reserve(n * (nnzL + nnzU + 1)); // global loop over the rows of the sparse matrix - for (int ii = 0; ii < n; ii++) + for (Index ii = 0; ii < n; ii++) { // 1 - copy the lower and the upper part of the row i of mat in the working vector u - int sizeu = 1; // number of nonzero elements in the upper part of the current row - int sizel = 0; // number of nonzero elements in the lower part of the current row + Index sizeu = 1; // number of nonzero elements in the upper part of the current row + Index sizel = 0; // number of nonzero elements in the lower part of the current row ju(ii) = ii; u(ii) = 0; jr(ii) = ii; @@ -285,7 +287,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii for (; j_it; ++j_it) { - int k = j_it.index(); + Index k = j_it.index(); if (k < ii) { // copy the lower part @@ -301,13 +303,13 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) else { // copy the upper part - int jpos = ii + sizeu; + Index jpos = ii + sizeu; ju(jpos) = k; u(jpos) = j_it.value(); jr(k) = jpos; ++sizeu; } - rownorm += internal::abs2(j_it.value()); + rownorm += numext::abs2(j_it.value()); } // 2 - detect possible zero row @@ -320,19 +322,19 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) rownorm = sqrt(rownorm); // 3 - eliminate the previous nonzero rows - int jj = 0; - int len = 0; + Index jj = 0; + Index len = 0; while (jj < sizel) { // In order to eliminate in the correct order, // we must select first the smallest column index among ju(jj:sizel) - int k; - int minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment + Index k; + Index minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment k += jj; if (minrow != ju(jj)) { // swap the two locations - int j = ju(jj); + Index j = ju(jj); swap(ju(jj), ju(k)); jr(minrow) = jj; jr(j) = k; swap(u(jj), u(k)); @@ -358,11 +360,11 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) for (; ki_it; ++ki_it) { Scalar prod = fact * ki_it.value(); - int j = ki_it.index(); - int jpos = jr(j); + Index j = ki_it.index(); + Index jpos = jr(j); if (jpos == -1) // fill-in element { - int newpos; + Index newpos; if (j >= ii) // dealing with the upper part { newpos = ii + sizeu; @@ -391,7 +393,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) } // end of the elimination on the row ii // reset the upper part of the pointer jr to zero - for(int k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; + for(Index k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1; // 4 - partially sort and insert the elements in the m_lu matrix @@ -400,11 +402,11 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) len = (std::min)(sizel, nnzL); typename Vector::SegmentReturnType ul(u.segment(0, sizel)); typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel)); - QuickSplit(ul, jul, len); + internal::QuickSplit(ul, jul, len); // store the largest m_fill elements of the L part m_lu.startVec(ii); - for(int k = 0; k < len; k++) + for(Index k = 0; k < len; k++) m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); // store the diagonal element @@ -416,7 +418,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) // sort the U-part of the row // apply the dropping rule first len = 0; - for(int k = 1; k < sizeu; k++) + for(Index k = 1; k < sizeu; k++) { if(abs(u(ii+k)) > m_droptol * rownorm ) { @@ -429,10 +431,10 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) len = (std::min)(sizeu, nnzU); typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1)); typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1)); - QuickSplit(uu, juu, len); + internal::QuickSplit(uu, juu, len); // store the largest elements of the U part - for(int k = ii + 1; k < ii + len; k++) + for(Index k = ii + 1; k < ii + len; k++) m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k); } @@ -463,4 +465,3 @@ struct solve_retval<IncompleteLUT<_MatrixType>, Rhs> } // end namespace Eigen #endif // EIGEN_INCOMPLETE_LUT_H - |