diff options
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers')
8 files changed, 205 insertions, 218 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h index facdaf890..a117fc155 100644 --- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h +++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h @@ -10,7 +10,7 @@ #ifndef EIGEN_BASIC_PRECONDITIONERS_H #define EIGEN_BASIC_PRECONDITIONERS_H -namespace Eigen { +namespace Eigen { /** \ingroup IterativeLinearSolvers_Module * \brief A preconditioner based on the digonal entries @@ -52,15 +52,15 @@ class DiagonalPreconditioner compute(mat); } - Index rows() const { return m_invdiag.size(); } - Index cols() const { return m_invdiag.size(); } - + EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_invdiag.size(); } + EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_invdiag.size(); } + template<typename MatType> DiagonalPreconditioner& analyzePattern(const MatType& ) { return *this; } - + template<typename MatType> DiagonalPreconditioner& factorize(const MatType& mat) { @@ -77,7 +77,7 @@ class DiagonalPreconditioner m_isInitialized = true; return *this; } - + template<typename MatType> DiagonalPreconditioner& compute(const MatType& mat) { @@ -99,7 +99,7 @@ class DiagonalPreconditioner && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived()); } - + ComputationInfo info() { return Success; } protected: @@ -121,7 +121,7 @@ class DiagonalPreconditioner * \implsparsesolverconcept * * The diagonal entries are pre-inverted and stored into a dense vector. - * + * * \sa class LeastSquaresConjugateGradient, class DiagonalPreconditioner */ template <typename _Scalar> @@ -146,7 +146,7 @@ class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar> { return *this; } - + template<typename MatType> LeastSquareDiagonalPreconditioner& factorize(const MatType& mat) { @@ -168,7 +168,7 @@ class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar> { for(Index j=0; j<mat.outerSize(); ++j) { - RealScalar sum = mat.innerVector(j).squaredNorm(); + RealScalar sum = mat.col(j).squaredNorm(); if(sum>RealScalar(0)) m_invdiag(j) = RealScalar(1)/sum; else @@ -178,13 +178,13 @@ class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar> Base::m_isInitialized = true; return *this; } - + template<typename MatType> LeastSquareDiagonalPreconditioner& compute(const MatType& mat) { return factorize(mat); } - + ComputationInfo info() { return Success; } protected: @@ -205,19 +205,19 @@ class IdentityPreconditioner template<typename MatrixType> explicit IdentityPreconditioner(const MatrixType& ) {} - + template<typename MatrixType> IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; } - + template<typename MatrixType> IdentityPreconditioner& factorize(const MatrixType& ) { return *this; } template<typename MatrixType> IdentityPreconditioner& compute(const MatrixType& ) { return *this; } - + template<typename Rhs> inline const Rhs& solve(const Rhs& b) const { return b; } - + ComputationInfo info() { return Success; } }; diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 454f46814..153acef65 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -191,32 +191,16 @@ public: /** \internal */ template<typename Rhs,typename Dest> - void _solve_with_guess_impl(const Rhs& b, Dest& x) const + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const { - bool failed = false; - for(Index j=0; j<b.cols(); ++j) - { - m_iterations = Base::maxIterations(); - m_error = Base::m_tolerance; - - typename Dest::ColXpr xj(x,j); - if(!internal::bicgstab(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error)) - failed = true; - } - m_info = failed ? NumericalIssue + m_iterations = Base::maxIterations(); + m_error = Base::m_tolerance; + + bool ret = internal::bicgstab(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error); + + m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence; - m_isInitialized = true; - } - - /** \internal */ - using Base::_solve_impl; - template<typename Rhs,typename Dest> - void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const - { - x.resize(this->rows(),b.cols()); - x.setZero(); - _solve_with_guess_impl(b,x); } protected: diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index 395daa8e4..5d8c6b433 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -50,7 +50,8 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, tol_error = 0; return; } - RealScalar threshold = tol*tol*rhsNorm2; + const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); + RealScalar threshold = numext::maxi(RealScalar(tol*tol*rhsNorm2),considerAsZero); RealScalar residualNorm2 = residual.squaredNorm(); if (residualNorm2 < threshold) { @@ -58,7 +59,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x, tol_error = sqrt(residualNorm2 / rhsNorm2); return; } - + VectorType p(n); p = precond.solve(residual); // initial search direction @@ -194,7 +195,7 @@ public: /** \internal */ template<typename Rhs,typename Dest> - void _solve_with_guess_impl(const Rhs& b, Dest& x) const + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const { typedef typename Base::MatrixWrapper MatrixWrapper; typedef typename Base::ActualMatrixType ActualMatrixType; @@ -210,31 +211,14 @@ public: RowMajorWrapper, typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type >::type SelfAdjointWrapper; + m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - for(Index j=0; j<b.cols(); ++j) - { - m_iterations = Base::maxIterations(); - m_error = Base::m_tolerance; - - typename Dest::ColXpr xj(x,j); - RowMajorWrapper row_mat(matrix()); - internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); - } - - m_isInitialized = true; + RowMajorWrapper row_mat(matrix()); + internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error); m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; } - - /** \internal */ - using Base::_solve_impl; - template<typename Rhs,typename Dest> - void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const - { - x.setZero(); - _solve_with_guess_impl(b.derived(),x); - } protected: diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h index e45c272b4..7803fd817 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h @@ -14,8 +14,8 @@ #include <vector> #include <list> -namespace Eigen { -/** +namespace Eigen { +/** * \brief Modified Incomplete Cholesky with dual threshold * * References : C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with @@ -41,28 +41,22 @@ namespace Eigen { * the info() method, then you can either increase the initial shift, or better use another preconditioning technique. * */ -template <typename Scalar, int _UpLo = Lower, typename _OrderingType = -#ifndef EIGEN_MPL2_ONLY -AMDOrdering<int> -#else -NaturalOrdering<int> -#endif -> +template <typename Scalar, int _UpLo = Lower, typename _OrderingType = AMDOrdering<int> > class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > { protected: typedef SparseSolverBase<IncompleteCholesky<Scalar,_UpLo,_OrderingType> > Base; using Base::m_isInitialized; public: - typedef typename NumTraits<Scalar>::Real RealScalar; + typedef typename NumTraits<Scalar>::Real RealScalar; typedef _OrderingType OrderingType; typedef typename OrderingType::PermutationType PermutationType; - typedef typename PermutationType::StorageIndex StorageIndex; + typedef typename PermutationType::StorageIndex StorageIndex; typedef SparseMatrix<Scalar,ColMajor,StorageIndex> FactorType; typedef Matrix<Scalar,Dynamic,1> VectorSx; typedef Matrix<RealScalar,Dynamic,1> VectorRx; typedef Matrix<StorageIndex,Dynamic, 1> VectorIx; - typedef std::vector<std::list<StorageIndex> > VectorList; + typedef std::vector<std::list<StorageIndex> > VectorList; enum { UpLo = _UpLo }; enum { ColsAtCompileTime = Dynamic, @@ -76,22 +70,22 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up * * \sa IncompleteCholesky(const MatrixType&) */ - IncompleteCholesky() : m_initialShift(1e-3),m_factorizationIsOk(false) {} - + IncompleteCholesky() : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) {} + /** Constructor computing the incomplete factorization for the given matrix \a matrix. */ template<typename MatrixType> - IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_factorizationIsOk(false) + IncompleteCholesky(const MatrixType& matrix) : m_initialShift(1e-3),m_analysisIsOk(false),m_factorizationIsOk(false) { compute(matrix); } - + /** \returns number of rows of the factored matrix */ - Index rows() const { return m_L.rows(); } - + EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_L.rows(); } + /** \returns number of columns of the factored matrix */ - Index cols() const { return m_L.cols(); } - + EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_L.cols(); } + /** \brief Reports whether previous computation was successful. * @@ -106,19 +100,19 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up eigen_assert(m_isInitialized && "IncompleteCholesky is not initialized."); return m_info; } - + /** \brief Set the initial shift parameter \f$ \sigma \f$. */ void setInitialShift(RealScalar shift) { m_initialShift = shift; } - + /** \brief Computes the fill reducing permutation vector using the sparsity pattern of \a mat */ template<typename MatrixType> void analyzePattern(const MatrixType& mat) { - OrderingType ord; + OrderingType ord; PermutationType pinv; - ord(mat.template selfadjointView<UpLo>(), pinv); + ord(mat.template selfadjointView<UpLo>(), pinv); if(pinv.size()>0) m_perm = pinv.inverse(); else m_perm.resize(0); m_L.resize(mat.rows(), mat.cols()); @@ -126,7 +120,7 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up m_isInitialized = true; m_info = Success; } - + /** \brief Performs the numerical factorization of the input matrix \a mat * * The method analyzePattern() or compute() must have been called beforehand @@ -136,7 +130,7 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up */ template<typename MatrixType> void factorize(const MatrixType& mat); - + /** Computes or re-computes the incomplete Cholesky factorization of the input matrix \a mat * * It is a shortcut for a sequential call to the analyzePattern() and factorize() methods. @@ -149,7 +143,7 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up analyzePattern(mat); factorize(mat); } - + // internal template<typename Rhs, typename Dest> void _solve_impl(const Rhs& b, Dest& x) const @@ -176,16 +170,16 @@ class IncompleteCholesky : public SparseSolverBase<IncompleteCholesky<Scalar,_Up protected: FactorType m_L; // The lower part stored in CSC - VectorRx m_scale; // The vector for scaling the matrix + VectorRx m_scale; // The vector for scaling the matrix RealScalar m_initialShift; // The initial shift parameter - bool m_analysisIsOk; - bool m_factorizationIsOk; + bool m_analysisIsOk; + bool m_factorizationIsOk; ComputationInfo m_info; - PermutationType m_perm; + PermutationType m_perm; private: - inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); -}; + inline void updateList(Ref<const VectorIx> colPtr, Ref<VectorIx> rowIdx, Ref<VectorSx> vals, const Index& col, const Index& jk, VectorIx& firstElt, VectorList& listCol); +}; // Based on the following paper: // C-J. Lin and J. J. Moré, Incomplete Cholesky Factorizations with @@ -196,10 +190,10 @@ template<typename _MatrixType> void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType& mat) { using std::sqrt; - eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); - + eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); + // Dropping strategy : Keep only the p largest elements per column, where p is the number of elements in the column of the original matrix. Other strategies will be added - + // Apply the fill-reducing permutation computed in analyzePattern() if (m_perm.rows() == mat.rows() ) // To detect the null permutation { @@ -212,8 +206,8 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType { m_L.template selfadjointView<Lower>() = mat.template selfadjointView<_UpLo>(); } - - Index n = m_L.cols(); + + Index n = m_L.cols(); Index nnz = m_L.nonZeros(); Map<VectorSx> vals(m_L.valuePtr(), nnz); //values Map<VectorIx> rowIdx(m_L.innerIndexPtr(), nnz); //Row indices @@ -225,9 +219,9 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType VectorIx col_pattern(n); col_pattern.fill(-1); StorageIndex col_nnz; - - - // Computes the scaling factors + + + // Computes the scaling factors m_scale.resize(n); m_scale.setZero(); for (Index j = 0; j < n; j++) @@ -237,7 +231,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType if(rowIdx[k]!=j) m_scale(rowIdx[k]) += numext::abs2(vals(k)); } - + m_scale = m_scale.cwiseSqrt().cwiseSqrt(); for (Index j = 0; j < n; ++j) @@ -247,8 +241,8 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType m_scale(j) = 1; // TODO disable scaling if not needed, i.e., if it is roughly uniform? (this will make solve() faster) - - // Scale and compute the shift for the matrix + + // Scale and compute the shift for the matrix RealScalar mindiag = NumTraits<RealScalar>::highest(); for (Index j = 0; j < n; j++) { @@ -259,7 +253,7 @@ void IncompleteCholesky<Scalar,_UpLo, OrderingType>::factorize(const _MatrixType } FactorType L_save = m_L; - + RealScalar shift = 0; if(mindiag <= RealScalar(0.)) shift = m_initialShift - mindiag; @@ -381,7 +375,7 @@ inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const if (jk < colPtr(col+1) ) { Index p = colPtr(col+1) - jk; - Index minpos; + Index minpos; rowIdx.segment(jk,p).minCoeff(&minpos); minpos += jk; if (rowIdx(minpos) != rowIdx(jk)) @@ -395,6 +389,6 @@ inline void IncompleteCholesky<Scalar,_UpLo, OrderingType>::updateList(Ref<const } } -} // end namespace Eigen +} // end namespace Eigen #endif diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index 338e6f10a..cdcf709eb 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -12,19 +12,19 @@ #define EIGEN_INCOMPLETE_LUT_H -namespace Eigen { +namespace Eigen { namespace internal { - + /** \internal - * Compute a quick-sort split of a vector + * 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 + * 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> Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) { @@ -34,15 +34,15 @@ Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) 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; + first = 0; + last = n-1; if (ncut < first || ncut > last ) return 0; - + do { - mid = first; - RealScalar abskey = abs(row(mid)); + mid = first; + RealScalar abskey = abs(row(mid)); for (Index j = first + 1; j <= last; j++) { if ( abs(row(j)) > abskey) { ++mid; @@ -53,12 +53,12 @@ Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) /* 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; + else if (mid < ncut ) first = mid + 1; } while (mid != ncut ); - - return 0; /* mid is equal to ncut */ + + return 0; /* mid is equal to ncut */ } }// end namespace internal @@ -71,23 +71,23 @@ Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) * * 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 + * 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 + * 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, + * 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: + * 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: @@ -115,28 +115,28 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageInd }; public: - + IncompleteLUT() : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10), m_analysisIsOk(false), m_factorizationIsOk(false) {} - + template<typename MatrixType> explicit 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) { eigen_assert(fillfactor != 0); - compute(mat); + compute(mat); } - - Index rows() const { return m_lu.rows(); } - - Index cols() const { return m_lu.cols(); } + + EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); } + + EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); } /** \brief Reports whether previous computation was successful. * - * \returns \c Success if computation was succesful, + * \returns \c Success if computation was successful, * \c NumericalIssue if the matrix.appears to be negative. */ ComputationInfo info() const @@ -144,36 +144,36 @@ class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar, _StorageInd eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); return m_info; } - + template<typename MatrixType> void analyzePattern(const MatrixType& amat); - + template<typename MatrixType> void factorize(const MatrixType& amat); - + /** * Compute an incomplete LU factorization with dual threshold on the matrix mat * No pivoting is done in this version - * + * **/ template<typename MatrixType> IncompleteLUT& compute(const MatrixType& amat) { - analyzePattern(amat); + analyzePattern(amat); factorize(amat); return *this; } - void setDroptol(const RealScalar& droptol); - void setFillfactor(int fillfactor); - + void setDroptol(const RealScalar& droptol); + void setFillfactor(int fillfactor); + template<typename Rhs, typename Dest> void _solve_impl(const Rhs& b, Dest& x) const { x = m_Pinv * b; x = m_lu.template triangularView<UnitLower>().solve(x); x = m_lu.template triangularView<Upper>().solve(x); - x = m_P * x; + x = m_P * x; } protected: @@ -200,22 +200,22 @@ protected: /** * Set control parameter droptol - * \param droptol Drop any element whose magnitude is less than this tolerance - **/ + * \param droptol Drop any element whose magnitude is less than this tolerance + **/ template<typename Scalar, typename StorageIndex> void IncompleteLUT<Scalar,StorageIndex>::setDroptol(const RealScalar& droptol) { - this->m_droptol = droptol; + this->m_droptol = droptol; } /** * Set control parameter fillfactor - * \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row. - **/ + * \param fillfactor This is used to compute the number @p fill_in of largest elements to keep on each row. + **/ template<typename Scalar, typename StorageIndex> void IncompleteLUT<Scalar,StorageIndex>::setFillfactor(int fillfactor) { - this->m_fillfactor = fillfactor; + this->m_fillfactor = fillfactor; } template <typename Scalar, typename StorageIndex> @@ -225,24 +225,15 @@ void IncompleteLUT<Scalar,StorageIndex>::analyzePattern(const _MatrixType& amat) // Compute the Fill-reducing permutation // Since ILUT does not perform any numerical pivoting, // it is highly preferable to keep the diagonal through symmetric permutations. -#ifndef EIGEN_MPL2_ONLY // To this end, let's symmetrize the pattern and perform AMD on it. SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat; SparseMatrix<Scalar,ColMajor, StorageIndex> mat2 = amat.transpose(); // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice. - // on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered... + // on the other hand for a really non-symmetric pattern, mat2*mat1 should be preferred... SparseMatrix<Scalar,ColMajor, StorageIndex> AtA = mat2 + mat1; AMDOrdering<StorageIndex> ordering; ordering(AtA,m_P); m_Pinv = m_P.inverse(); // cache the inverse permutation -#else - // If AMD is not available, (MPL2-only), then let's use the slower COLAMD routine. - SparseMatrix<Scalar,ColMajor, StorageIndex> mat1 = amat; - COLAMDOrdering<StorageIndex> ordering; - ordering(mat1,m_Pinv); - m_P = m_Pinv.inverse(); -#endif - m_analysisIsOk = true; m_factorizationIsOk = false; m_isInitialized = true; diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 7c2326eb7..28a0c5109 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -10,7 +10,7 @@ #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H #define EIGEN_ITERATIVE_SOLVER_BASE_H -namespace Eigen { +namespace Eigen { namespace internal { @@ -145,7 +145,7 @@ class IterativeSolverBase : public SparseSolverBase<Derived> protected: typedef SparseSolverBase<Derived> Base; using Base::m_isInitialized; - + public: typedef typename internal::traits<Derived>::MatrixType MatrixType; typedef typename internal::traits<Derived>::Preconditioner Preconditioner; @@ -169,10 +169,10 @@ public: } /** Initialize the solver with matrix \a A for further \c Ax=b solving. - * + * * This constructor is a shortcut for the default constructor followed * by a call to compute(). - * + * * \warning this class stores a reference to the matrix A as well as some * precomputed values that depend on it. Therefore, if \a A is changed * this class becomes invalid. Call compute() to update it with the new @@ -187,7 +187,7 @@ public: } ~IterativeSolverBase() {} - + /** Initializes the iterative solver for the sparsity pattern of the matrix \a A for further solving \c Ax=b problems. * * Currently, this function mostly calls analyzePattern on the preconditioner. In the future @@ -203,7 +203,7 @@ public: m_info = m_preconditioner.info(); return derived(); } - + /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems. * * Currently, this function mostly calls factorize on the preconditioner. @@ -216,7 +216,7 @@ public: template<typename MatrixDerived> Derived& factorize(const EigenBase<MatrixDerived>& A) { - eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); + eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); grab(A.derived()); m_preconditioner.factorize(matrix()); m_factorizationIsOk = true; @@ -247,16 +247,16 @@ public: } /** \internal */ - Index rows() const { return matrix().rows(); } + EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return matrix().rows(); } /** \internal */ - Index cols() const { return matrix().cols(); } + EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return matrix().cols(); } /** \returns the tolerance threshold used by the stopping criteria. * \sa setTolerance() */ RealScalar tolerance() const { return m_tolerance; } - + /** Sets the tolerance threshold used by the stopping criteria. * * This value is used as an upper bound to the relative residual error: |Ax-b|/|b|. @@ -270,19 +270,19 @@ public: /** \returns a read-write reference to the preconditioner for custom configuration. */ Preconditioner& preconditioner() { return m_preconditioner; } - + /** \returns a read-only reference to the preconditioner. */ const Preconditioner& preconditioner() const { return m_preconditioner; } /** \returns the max number of iterations. - * It is either the value setted by setMaxIterations or, by default, + * It is either the value set by setMaxIterations or, by default, * twice the number of columns of the matrix. */ Index maxIterations() const { return (m_maxIterations<0) ? 2*matrix().cols() : m_maxIterations; } - + /** Sets the max number of iterations. * Default is twice the number of columns of the matrix. */ @@ -328,13 +328,13 @@ public: eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); return m_info; } - + /** \internal */ template<typename Rhs, typename DestDerived> - void _solve_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const + void _solve_with_guess_impl(const Rhs& b, SparseMatrixBase<DestDerived> &aDest) const { eigen_assert(rows()==b.rows()); - + Index rhsCols = b.cols(); Index size = b.rows(); DestDerived& dest(aDest.derived()); @@ -344,15 +344,65 @@ public: // We do not directly fill dest because sparse expressions have to be free of aliasing issue. // For non square least-square problems, b and dest might not have the same size whereas they might alias each-other. typename DestDerived::PlainObject tmp(cols(),rhsCols); + ComputationInfo global_info = Success; for(Index k=0; k<rhsCols; ++k) { tb = b.col(k); - tx = derived().solve(tb); + tx = dest.col(k); + derived()._solve_vector_with_guess_impl(tb,tx); tmp.col(k) = tx.sparseView(0); + + // The call to _solve_vector_with_guess_impl updates m_info, so if it failed for a previous column + // we need to restore it to the worst value. + if(m_info==NumericalIssue) + global_info = NumericalIssue; + else if(m_info==NoConvergence) + global_info = NoConvergence; } + m_info = global_info; dest.swap(tmp); } + template<typename Rhs, typename DestDerived> + typename internal::enable_if<Rhs::ColsAtCompileTime!=1 && DestDerived::ColsAtCompileTime!=1>::type + _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &aDest) const + { + eigen_assert(rows()==b.rows()); + + Index rhsCols = b.cols(); + DestDerived& dest(aDest.derived()); + ComputationInfo global_info = Success; + for(Index k=0; k<rhsCols; ++k) + { + typename DestDerived::ColXpr xk(dest,k); + typename Rhs::ConstColXpr bk(b,k); + derived()._solve_vector_with_guess_impl(bk,xk); + + // The call to _solve_vector_with_guess updates m_info, so if it failed for a previous column + // we need to restore it to the worst value. + if(m_info==NumericalIssue) + global_info = NumericalIssue; + else if(m_info==NoConvergence) + global_info = NoConvergence; + } + m_info = global_info; + } + + template<typename Rhs, typename DestDerived> + typename internal::enable_if<Rhs::ColsAtCompileTime==1 || DestDerived::ColsAtCompileTime==1>::type + _solve_with_guess_impl(const Rhs& b, MatrixBase<DestDerived> &dest) const + { + derived()._solve_vector_with_guess_impl(b,dest.derived()); + } + + /** \internal default initial guess = 0 */ + template<typename Rhs,typename Dest> + void _solve_impl(const Rhs& b, Dest& x) const + { + x.setZero(); + derived()._solve_with_guess_impl(b,x); + } + protected: void init() { @@ -370,19 +420,19 @@ protected: { return m_matrixWrapper.matrix(); } - + template<typename InputType> void grab(const InputType &A) { m_matrixWrapper.grab(A); } - + MatrixWrapper m_matrixWrapper; Preconditioner m_preconditioner; Index m_maxIterations; RealScalar m_tolerance; - + mutable RealScalar m_error; mutable Index m_iterations; mutable ComputationInfo m_info; diff --git a/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h index 0aea0e099..203fd0ec6 100644 --- a/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/LeastSquareConjugateGradient.h @@ -182,32 +182,14 @@ public: /** \internal */ template<typename Rhs,typename Dest> - void _solve_with_guess_impl(const Rhs& b, Dest& x) const + void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; - for(Index j=0; j<b.cols(); ++j) - { - m_iterations = Base::maxIterations(); - m_error = Base::m_tolerance; - - typename Dest::ColXpr xj(x,j); - internal::least_square_conjugate_gradient(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error); - } - - m_isInitialized = true; + internal::least_square_conjugate_gradient(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error); m_info = m_error <= Base::m_tolerance ? Success : NoConvergence; } - - /** \internal */ - using Base::_solve_impl; - template<typename Rhs,typename Dest> - void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const - { - x.setZero(); - _solve_with_guess_impl(b.derived(),x); - } }; diff --git a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h index 0ace45177..7b8965754 100644 --- a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h +++ b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h @@ -13,7 +13,7 @@ namespace Eigen { template<typename Decomposition, typename RhsType, typename GuessType> class SolveWithGuess; - + /** \class SolveWithGuess * \ingroup IterativeLinearSolvers_Module * @@ -45,13 +45,15 @@ public: typedef typename internal::traits<SolveWithGuess>::PlainObject PlainObject; typedef typename internal::generic_xpr_base<SolveWithGuess<Decomposition,RhsType,GuessType>, MatrixXpr, typename internal::traits<RhsType>::StorageKind>::type Base; typedef typename internal::ref_selector<SolveWithGuess>::type Nested; - + SolveWithGuess(const Decomposition &dec, const RhsType &rhs, const GuessType &guess) : m_dec(dec), m_rhs(rhs), m_guess(guess) {} - - EIGEN_DEVICE_FUNC Index rows() const { return m_dec.cols(); } - EIGEN_DEVICE_FUNC Index cols() const { return m_rhs.cols(); } + + EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR + Index rows() const EIGEN_NOEXCEPT { return m_dec.cols(); } + EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR + Index cols() const EIGEN_NOEXCEPT { return m_rhs.cols(); } EIGEN_DEVICE_FUNC const Decomposition& dec() const { return m_dec; } EIGEN_DEVICE_FUNC const RhsType& rhs() const { return m_rhs; } @@ -61,7 +63,7 @@ protected: const Decomposition &m_dec; const RhsType &m_rhs; const GuessType &m_guess; - + private: Scalar coeff(Index row, Index col) const; Scalar coeff(Index i) const; @@ -85,8 +87,8 @@ struct evaluator<SolveWithGuess<Decomposition,RhsType, GuessType> > m_result = solve.guess(); solve.dec()._solve_with_guess_impl(solve.rhs(), m_result); } - -protected: + +protected: PlainObject m_result; }; @@ -108,7 +110,7 @@ struct Assignment<DstXprType, SolveWithGuess<DecType,RhsType,GuessType>, interna } }; -} // end namepsace internal +} // end namespace internal } // end namespace Eigen |