diff options
Diffstat (limited to 'extern/Eigen3/Eigen/src/IterativeLinearSolvers')
4 files changed, 155 insertions, 119 deletions
diff --git a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 126341be8d8..6fc6ab85225 100644 --- a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -39,10 +39,17 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, int maxIters = iters; int n = mat.cols(); + x = precond.solve(x); VectorType r = rhs - mat * x; VectorType r0 = r; RealScalar r0_sqnorm = r0.squaredNorm(); + RealScalar rhs_sqnorm = rhs.squaredNorm(); + if(rhs_sqnorm == 0) + { + x.setZero(); + return true; + } Scalar rho = 1; Scalar alpha = 1; Scalar w = 1; @@ -55,13 +62,22 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, RealScalar tol2 = tol*tol; int i = 0; + int restarts = 0; - while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters ) + while ( r.squaredNorm()/rhs_sqnorm > tol2 && i<maxIters ) { Scalar rho_old = rho; rho = r0.dot(r); - if (rho == Scalar(0)) return false; /* New search directions cannot be found */ + if (internal::isMuchSmallerThan(rho,r0_sqnorm)) + { + // The new residual vector became too orthogonal to the arbitrarily choosen direction r0 + // Let's restart with a new r0: + r0 = r; + rho = r0_sqnorm = r.squaredNorm(); + if(restarts++ == 0) + i = 0; + } Scalar beta = (rho/rho_old) * (alpha / w); p = r + beta * (p - w * v); @@ -75,12 +91,16 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, z = precond.solve(s); t.noalias() = mat * z; - w = t.dot(s) / t.squaredNorm(); + RealScalar tmp = t.squaredNorm(); + if(tmp>RealScalar(0)) + w = t.dot(s) / tmp; + else + w = Scalar(0); x += alpha * y + w * z; r = s - w * t; ++i; } - tol_error = sqrt(r.squaredNorm()/r0_sqnorm); + tol_error = sqrt(r.squaredNorm()/rhs_sqnorm); iters = i; return true; } @@ -223,7 +243,8 @@ public: template<typename Rhs,typename Dest> void _solve(const Rhs& b, Dest& x) const { - x.setZero(); +// x.setZero(); + x = b; _solveWithGuess(b,x); } diff --git a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index f64f2534d28..a74a8155e68 100644 --- a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -41,15 +41,29 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, int n = mat.cols(); VectorType residual = rhs - mat * x; //initial residual - VectorType p(n); + RealScalar rhsNorm2 = rhs.squaredNorm(); + if(rhsNorm2 == 0) + { + x.setZero(); + iters = 0; + tol_error = 0; + return; + } + RealScalar threshold = tol*tol*rhsNorm2; + RealScalar residualNorm2 = residual.squaredNorm(); + if (residualNorm2 < threshold) + { + iters = 0; + tol_error = sqrt(residualNorm2 / rhsNorm2); + return; + } + + VectorType p(n); p = precond.solve(residual); //initial search direction VectorType z(n), tmp(n); - RealScalar absNew = internal::real(residual.dot(p)); // the square of the absolute value of r scaled by invM - RealScalar rhsNorm2 = rhs.squaredNorm(); - RealScalar residualNorm2 = 0; - RealScalar threshold = tol*tol*rhsNorm2; + RealScalar absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM int i = 0; while(i < maxIters) { @@ -66,7 +80,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, z = precond.solve(residual); // approximately solve for "A z = residual" RealScalar absOld = absNew; - absNew = internal::real(residual.dot(z)); // update the absolute value of r + absNew = numext::real(residual.dot(z)); // update the absolute value of r RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction p = z + beta * p; // update search direction i++; diff --git a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 224304f0eb8..b55afc13636 100644 --- a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/extern/Eigen3/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 - diff --git a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 11706cebabd..2036922d69c 100644 --- a/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/extern/Eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -120,7 +120,7 @@ public: RealScalar tolerance() const { return m_tolerance; } /** Sets the tolerance threshold used by the stopping criteria */ - Derived& setTolerance(RealScalar tolerance) + Derived& setTolerance(const RealScalar& tolerance) { m_tolerance = tolerance; return derived(); |