diff options
author | 2014-09-21 14:02:51 +0300 | |
---|---|---|
committer | 2014-09-21 14:02:51 +0300 | |
commit | 60e093a9dce2f8d4c0f3b2ea3e0386d5f01bff8d (patch) | |
tree | 05442eeff0bcfe7fe85ce59cf5fa72aa06ee2a07 /Eigen/src/IterativeLinearSolvers | |
parent | 56408504e4e3fa5f9c59d9edac14ca1ba1255e5a (diff) | |
parent | 03dd4dd91a5d8963f56eebe3b9d2eb924bc06e02 (diff) |
Merged eigen/eigen into default
Diffstat (limited to 'Eigen/src/IterativeLinearSolvers')
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h | 25 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 49 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/ConjugateGradient.h | 43 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IncompleteLUT.h | 40 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h | 63 | ||||
-rw-r--r-- | Eigen/src/IterativeLinearSolvers/SolveWithGuess.h | 114 |
6 files changed, 159 insertions, 175 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h index 1f3c060d0..98b169868 100644 --- a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h +++ b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -80,19 +80,20 @@ class DiagonalPreconditioner return factorize(mat); } + /** \internal */ template<typename Rhs, typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, Dest& x) const { x = m_invdiag.array() * b.array() ; } - template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs> + template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs> solve(const MatrixBase<Rhs>& b) const { eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized."); eigen_assert(m_invdiag.size()==b.rows() && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived()); + return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived()); } protected: @@ -100,22 +101,6 @@ class DiagonalPreconditioner bool m_isInitialized; }; -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs> - : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs> -{ - typedef DiagonalPreconditioner<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} /** \ingroup IterativeLinearSolvers_Module * \brief A naive preconditioner which approximates any matrix as the identity matrix diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 27824b9d5..051940dc7 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla @@ -181,26 +181,10 @@ public: BiCGSTAB(const MatrixType& A) : Base(A) {} ~BiCGSTAB() {} - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A - * \a x0 as an initial solution. - * - * \sa compute() - */ - template<typename Rhs,typename Guess> - inline const internal::solve_retval_with_guess<BiCGSTAB, Rhs, Guess> - solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const - { - eigen_assert(m_isInitialized && "BiCGSTAB is not initialized."); - eigen_assert(Base::rows()==b.rows() - && "BiCGSTAB::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval_with_guess - <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0); - } - + /** \internal */ template<typename Rhs,typename Dest> - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { bool failed = false; for(int j=0; j<b.cols(); ++j) @@ -219,36 +203,19 @@ public: } /** \internal */ + using Base::_solve_impl; template<typename Rhs,typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const { -// x.setZero(); - x = b; - _solveWithGuess(b,x); + // x.setZero(); + x = b; + _solve_with_guess_impl(b,x); } protected: }; - -namespace internal { - - template<typename _MatrixType, typename _Preconditioner, typename Rhs> -struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> - : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> -{ - typedef BiCGSTAB<_MatrixType, _Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_BICGSTAB_H diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index 3ce517940..f72cf86a5 100644 --- a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h +++ b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -192,26 +192,10 @@ public: ConjugateGradient(const MatrixType& A) : Base(A) {} ~ConjugateGradient() {} - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A - * \a x0 as an initial solution. - * - * \sa compute() - */ - template<typename Rhs,typename Guess> - inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess> - solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const - { - eigen_assert(m_isInitialized && "ConjugateGradient is not initialized."); - eigen_assert(Base::rows()==b.rows() - && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval_with_guess - <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0); - } /** \internal */ template<typename Rhs,typename Dest> - void _solveWithGuess(const Rhs& b, Dest& x) const + void _solve_with_guess_impl(const Rhs& b, Dest& x) const { m_iterations = Base::maxIterations(); m_error = Base::m_tolerance; @@ -231,35 +215,18 @@ public: } /** \internal */ + using Base::_solve_impl; template<typename Rhs,typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const { x.setOnes(); - _solveWithGuess(b,x); + _solve_with_guess_impl(b.derived(),x); } protected: }; - -namespace internal { - -template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs> -struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs> - : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs> -{ - typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_CONJUGATE_GRADIENT_H diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index b55afc136..7adbbc489 100644 --- a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h +++ b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> +// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -93,8 +94,12 @@ Index QuickSplit(VectorV &row, VectorI &ind, Index ncut) * http://comments.gmane.org/gmane.comp.lib.eigen/3302 */ template <typename _Scalar> -class IncompleteLUT : internal::noncopyable +class IncompleteLUT : public SparseSolverBase<IncompleteLUT<_Scalar> > { + protected: + typedef SparseSolverBase<IncompleteLUT<_Scalar> > Base; + using Base::m_isInitialized; + public: typedef _Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar,Dynamic,1> Vector; @@ -107,13 +112,13 @@ class IncompleteLUT : internal::noncopyable IncompleteLUT() : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10), - m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false) + m_analysisIsOk(false), m_factorizationIsOk(false) {} template<typename MatrixType> 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) + m_analysisIsOk(false),m_factorizationIsOk(false) { eigen_assert(fillfactor != 0); compute(mat); @@ -158,7 +163,7 @@ class IncompleteLUT : internal::noncopyable void setFillfactor(int fillfactor); template<typename Rhs, typename Dest> - void _solve(const Rhs& b, Dest& x) const + void _solve_impl(const Rhs& b, Dest& x) const { x = m_Pinv * b; x = m_lu.template triangularView<UnitLower>().solve(x); @@ -166,15 +171,6 @@ class IncompleteLUT : internal::noncopyable x = m_P * x; } - template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "IncompleteLUT is not initialized."); - eigen_assert(cols()==b.rows() - && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived()); - } - protected: /** keeps off-diagonal entries; drops diagonal entries */ @@ -192,7 +188,6 @@ protected: int m_fillfactor; bool m_analysisIsOk; bool m_factorizationIsOk; - bool m_isInitialized; ComputationInfo m_info; PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation @@ -445,23 +440,6 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) m_info = Success; } -namespace internal { - -template<typename _MatrixType, typename Rhs> -struct solve_retval<IncompleteLUT<_MatrixType>, Rhs> - : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs> -{ - typedef IncompleteLUT<_MatrixType> Dec; - EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec()._solve(rhs(),dst); - } -}; - -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_INCOMPLETE_LUT_H diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 2036922d6..fd9285087 100644 --- a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h +++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -18,8 +18,12 @@ namespace Eigen { * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner */ template< typename Derived> -class IterativeSolverBase : internal::noncopyable +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; @@ -29,8 +33,7 @@ public: public: - Derived& derived() { return *static_cast<Derived*>(this); } - const Derived& derived() const { return *static_cast<const Derived*>(this); } + using Base::derived; /** Default constructor. */ IterativeSolverBase() @@ -159,31 +162,18 @@ public: return m_error; } - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. - * - * \sa compute() - */ - template<typename Rhs> inline const internal::solve_retval<Derived, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); - eigen_assert(rows()==b.rows() - && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<Derived, Rhs>(derived(), b.derived()); - } - - /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. + /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A + * and \a x0 as an initial solution. * - * \sa compute() + * \sa solve(), compute() */ - template<typename Rhs> - inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs> - solve(const SparseMatrixBase<Rhs>& b) const + template<typename Rhs,typename Guess> + inline const SolveWithGuess<Derived, Rhs, Guess> + solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const { - eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized."); - eigen_assert(rows()==b.rows() - && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); - return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived()); + eigen_assert(m_isInitialized && "Solver is not initialized."); + eigen_assert(derived().rows()==b.rows() && "solve(): invalid number of rows of the right hand side matrix b"); + return SolveWithGuess<Derived, Rhs, Guess>(derived(), b.derived(), x0); } /** \returns Success if the iterations converged, and NoConvergence otherwise. */ @@ -195,7 +185,7 @@ public: /** \internal */ template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> - void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + void _solve_impl(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const { eigen_assert(rows()==b.rows()); @@ -229,26 +219,9 @@ protected: mutable RealScalar m_error; mutable int m_iterations; mutable ComputationInfo m_info; - mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk; -}; - -namespace internal { - -template<typename Derived, typename Rhs> -struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs> - : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs> -{ - typedef IterativeSolverBase<Derived> Dec; - EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) - - template<typename Dest> void evalTo(Dest& dst) const - { - dec().derived()._solve_sparse(rhs(),dst); - } + mutable bool m_analysisIsOk, m_factorizationIsOk; }; -} // end namespace internal - } // end namespace Eigen #endif // EIGEN_ITERATIVE_SOLVER_BASE_H diff --git a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h new file mode 100644 index 000000000..46dd5babe --- /dev/null +++ b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h @@ -0,0 +1,114 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_SOLVEWITHGUESS_H +#define EIGEN_SOLVEWITHGUESS_H + +namespace Eigen { + +template<typename Decomposition, typename RhsType, typename GuessType> class SolveWithGuess; + +/** \class SolveWithGuess + * \ingroup IterativeLinearSolvers_Module + * + * \brief Pseudo expression representing a solving operation + * + * \tparam Decomposition the type of the matrix or decomposion object + * \tparam Rhstype the type of the right-hand side + * + * This class represents an expression of A.solve(B) + * and most of the time this is the only way it is used. + * + */ +namespace internal { + + +template<typename Decomposition, typename RhsType, typename GuessType> +struct traits<SolveWithGuess<Decomposition, RhsType, GuessType> > + : traits<Solve<Decomposition,RhsType> > +{}; + +} + + +template<typename Decomposition, typename RhsType, typename GuessType> +class SolveWithGuess : public internal::generic_xpr_base<SolveWithGuess<Decomposition,RhsType,GuessType>, MatrixXpr, typename internal::traits<RhsType>::StorageKind>::type +{ +public: + typedef typename RhsType::Index Index; + 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; + + 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 const Decomposition& dec() const { return m_dec; } + EIGEN_DEVICE_FUNC const RhsType& rhs() const { return m_rhs; } + EIGEN_DEVICE_FUNC const GuessType& guess() const { return m_guess; } + +protected: + const Decomposition &m_dec; + const RhsType &m_rhs; + const GuessType &m_guess; + + typedef typename internal::traits<SolveWithGuess>::Scalar Scalar; + +private: + Scalar coeff(Index row, Index col) const; + Scalar coeff(Index i) const; +}; + +namespace internal { + +// Evaluator of SolveWithGuess -> eval into a temporary +template<typename Decomposition, typename RhsType, typename GuessType> +struct evaluator<SolveWithGuess<Decomposition,RhsType, GuessType> > + : public evaluator<typename SolveWithGuess<Decomposition,RhsType,GuessType>::PlainObject>::type +{ + typedef SolveWithGuess<Decomposition,RhsType,GuessType> SolveType; + typedef typename SolveType::PlainObject PlainObject; + typedef typename evaluator<PlainObject>::type Base; + + typedef evaluator type; + typedef evaluator nestedType; + + evaluator(const SolveType& solve) + : m_result(solve.rows(), solve.cols()) + { + ::new (static_cast<Base*>(this)) Base(m_result); + solve.dec()._solve_with_guess_impl(solve.rhs(), m_result, solve().guess()); + } + +protected: + PlainObject m_result; +}; + +// Specialization for "dst = dec.solve(rhs)" +// NOTE we need to specialize it for Dense2Dense to avoid ambiguous specialization error and a Sparse2Sparse specialization must exist somewhere +template<typename DstXprType, typename DecType, typename RhsType, typename GuessType, typename Scalar> +struct Assignment<DstXprType, SolveWithGuess<DecType,RhsType,GuessType>, internal::assign_op<Scalar>, Dense2Dense, Scalar> +{ + typedef Solve<DecType,RhsType> SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &) + { + // FIXME shall we resize dst here? + dst = src.guess(); + src.dec()._solve_with_guess_impl(src.rhs(), dst/*, src.guess()*/); + } +}; + +} // end namepsace internal + +} // end namespace Eigen + +#endif // EIGEN_SOLVEWITHGUESS_H |