diff options
20 files changed, 565 insertions, 131 deletions
diff --git a/Eigen/IterativeLinearSolvers b/Eigen/IterativeLinearSolvers index 0f4159dc1..4df2c5a14 100644 --- a/Eigen/IterativeLinearSolvers +++ b/Eigen/IterativeLinearSolvers @@ -26,8 +26,12 @@ * \endcode */ +#ifndef EIGEN_TEST_EVALUATORS #include "src/misc/Solve.h" #include "src/misc/SparseSolve.h" +#else +#include "src/IterativeLinearSolvers/SolveWithGuess.h" +#endif #include "src/IterativeLinearSolvers/IterativeSolverBase.h" #include "src/IterativeLinearSolvers/BasicPreconditioners.h" diff --git a/Eigen/SparseCholesky b/Eigen/SparseCholesky index 9f5056aa1..2c4f66105 100644 --- a/Eigen/SparseCholesky +++ b/Eigen/SparseCholesky @@ -34,8 +34,11 @@ #error The SparseCholesky module has nothing to offer in MPL2 only mode #endif +#ifndef EIGEN_TEST_EVALUATORS #include "src/misc/Solve.h" #include "src/misc/SparseSolve.h" +#endif + #include "src/SparseCholesky/SimplicialCholesky.h" #ifndef EIGEN_MPL2_ONLY diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 41cae58b0..b68c8fa8a 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -58,6 +58,7 @@ struct Sparse {}; #include "src/SparseCore/TriangularSolver.h" #include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseFuzzy.h" +#include "src/SparseCore/SparseSolverBase.h" #include "src/Core/util/ReenableStupidWarnings.h" diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index c449960de..4416c5308 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -157,8 +157,12 @@ enum CholmodMode { * \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT */ template<typename _MatrixType, int _UpLo, typename Derived> -class CholmodBase : internal::noncopyable +class CholmodBase : public SparseSolverBase<Derived> { + protected: + typedef SparseSolverBase<Derived> Base; + using Base::derived; + using Base::m_isInitialized; public: typedef _MatrixType MatrixType; enum { UpLo = _UpLo }; @@ -170,14 +174,14 @@ class CholmodBase : internal::noncopyable public: CholmodBase() - : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) + : m_cholmodFactor(0), m_info(Success) { m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0); cholmod_start(&m_cholmod); } CholmodBase(const MatrixType& matrix) - : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) + : m_cholmodFactor(0), m_info(Success) { m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0); cholmod_start(&m_cholmod); @@ -194,9 +198,6 @@ class CholmodBase : internal::noncopyable inline Index cols() const { return m_cholmodFactor->n; } inline Index rows() const { return m_cholmodFactor->n; } - Derived& derived() { return *static_cast<Derived*>(this); } - const Derived& derived() const { return *static_cast<const Derived*>(this); } - /** \brief Reports whether previous computation was successful. * * \returns \c Success if computation was succesful, @@ -216,6 +217,7 @@ class CholmodBase : internal::noncopyable return derived(); } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -235,14 +237,15 @@ class CholmodBase : internal::noncopyable * \sa compute() */ template<typename Rhs> - inline const internal::sparse_solve_retval<CholmodBase, Rhs> + inline const internal::sparse_solve_retval<Derived, Rhs> solve(const SparseMatrixBase<Rhs>& b) const { eigen_assert(m_isInitialized && "LLT is not initialized."); eigen_assert(rows()==b.rows() && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b"); - return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived()); + return internal::sparse_solve_retval<Derived, Rhs>(derived(), b.derived()); } +#endif // EIGEN_TEST_EVALUATORS /** Performs a symbolic decomposition on the sparsity pattern of \a matrix. * @@ -290,7 +293,7 @@ class CholmodBase : internal::noncopyable #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ template<typename Rhs,typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const { eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); const Index size = m_cholmodFactor->n; @@ -312,7 +315,7 @@ class CholmodBase : internal::noncopyable /** \internal */ template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex> - void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + void _solve_impl(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const { eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); const Index size = m_cholmodFactor->n; @@ -357,7 +360,6 @@ class CholmodBase : internal::noncopyable cholmod_factor* m_cholmodFactor; RealScalar m_shiftOffset[2]; mutable ComputationInfo m_info; - bool m_isInitialized; int m_factorizationIsOk; int m_analysisIsOk; }; @@ -572,6 +574,7 @@ class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecom } }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs> @@ -583,7 +586,7 @@ struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; @@ -596,11 +599,12 @@ struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; } // end namespace internal +#endif } // end namespace Eigen diff --git a/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h b/Eigen/src/IterativeLinearSolvers/BasicPreconditioners.h index 1f3c060d0..92af28cc8 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,12 +80,14 @@ 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() ; } +#ifndef EIGEN_TEST_EVALUATORS template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs> solve(const MatrixBase<Rhs>& b) const { @@ -94,12 +96,23 @@ class DiagonalPreconditioner && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b"); return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived()); } +#else + 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 Solve<DiagonalPreconditioner, Rhs>(*this, b.derived()); + } +#endif protected: Vector m_invdiag; bool m_isInitialized; }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Rhs> @@ -111,11 +124,12 @@ struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; } +#endif // EIGEN_TEST_EVALUATORS /** \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..2cad946d3 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,7 +181,8 @@ public: BiCGSTAB(const MatrixType& A) : Base(A) {} ~BiCGSTAB() {} - + +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A * \a x0 as an initial solution. * @@ -197,10 +198,26 @@ public: return internal::solve_retval_with_guess <BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0); } - +#else + /** \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 SolveWithGuess<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 SolveWithGuess<BiCGSTAB, Rhs, Guess>(*this, b.derived(), x0); + } +#endif + /** \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,22 +236,23 @@ 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: }; - +#ifndef EIGEN_TEST_EVALUATORS namespace internal { - template<typename _MatrixType, typename _Preconditioner, typename Rhs> +template<typename _MatrixType, typename _Preconditioner, typename Rhs> struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> : solve_retval_base<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> { @@ -243,11 +261,12 @@ struct solve_retval<BiCGSTAB<_MatrixType, _Preconditioner>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; } // end namespace internal +#endif } // end namespace Eigen diff --git a/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h index 3ce517940..000780eff 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 @@ -193,6 +193,7 @@ public: ~ConjugateGradient() {} +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A * \a x0 as an initial solution. * @@ -208,10 +209,26 @@ public: return internal::solve_retval_with_guess <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0); } +#else + /** \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 SolveWithGuess<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 SolveWithGuess<ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0); + } +#endif /** \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,18 +248,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.setOnes(); - _solveWithGuess(b,x); + _solve_with_guess_impl(b.derived(),x); } protected: }; - +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs> @@ -254,11 +272,12 @@ struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; } // end namespace internal +#endif } // end namespace Eigen diff --git a/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h b/Eigen/src/IterativeLinearSolvers/IncompleteLUT.h index b55afc136..a39ed7748 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 @@ -158,7 +159,11 @@ class IncompleteLUT : internal::noncopyable void setFillfactor(int fillfactor); template<typename Rhs, typename Dest> +#ifndef EIGEN_TEST_EVALUATORS void _solve(const Rhs& b, Dest& x) const +#else + void _solve_impl(const Rhs& b, Dest& x) const +#endif { x = m_Pinv * b; x = m_lu.template triangularView<UnitLower>().solve(x); @@ -166,14 +171,25 @@ class IncompleteLUT : internal::noncopyable x = m_P * x; } +#ifndef EIGEN_TEST_EVALUATORS template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs> - solve(const MatrixBase<Rhs>& b) const + 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()); } +#else + template<typename Rhs> inline const Solve<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 Solve<IncompleteLUT, Rhs>(*this, b.derived()); + } +#endif protected: @@ -445,6 +461,7 @@ void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat) m_info = Success; } +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Rhs> @@ -461,6 +478,7 @@ struct solve_retval<IncompleteLUT<_MatrixType>, Rhs> }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h index 2036922d6..2fc1a511b 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,6 +162,7 @@ public: return m_error; } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -171,7 +175,9 @@ public: && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); return internal::solve_retval<Derived, Rhs>(derived(), b.derived()); } +#endif +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -185,6 +191,7 @@ public: && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b"); return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived()); } +#endif // EIGEN_TEST_EVALUATORS /** \returns Success if the iterations converged, and NoConvergence otherwise. */ ComputationInfo info() const @@ -193,6 +200,7 @@ public: return m_info; } +#ifndef EIGEN_TEST_EVALUATORS /** \internal */ template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const @@ -210,6 +218,25 @@ public: dest.col(k) = tx.sparseView(0); } } +#else + /** \internal */ + template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> + void _solve_impl(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const + { + eigen_assert(rows()==b.rows()); + + int rhsCols = b.cols(); + int size = b.rows(); + Eigen::Matrix<DestScalar,Dynamic,1> tb(size); + Eigen::Matrix<DestScalar,Dynamic,1> tx(size); + for(int k=0; k<rhsCols; ++k) + { + tb = b.col(k); + tx = derived().solve(tb); + dest.col(k) = tx.sparseView(0); + } + } +#endif // EIGEN_TEST_EVALUATORS protected: void init() @@ -229,9 +256,10 @@ protected: mutable RealScalar m_error; mutable int m_iterations; mutable ComputationInfo m_info; - mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk; + mutable bool m_analysisIsOk, m_factorizationIsOk; }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename Derived, typename Rhs> @@ -248,6 +276,7 @@ struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs> }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen 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 diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index 8a546dc2f..0dc5c6a6f 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -125,9 +125,15 @@ namespace internal // This is the base class to interface with PaStiX functions. // Users should not used this class directly. template <class Derived> -class PastixBase : internal::noncopyable +class PastixBase : public SparseSolverBase<Derived> { + protected: + typedef SparseSolverBase<Derived> Base; + using Base::derived; + using Base::m_isInitialized; public: + using Base::_solve_impl; + typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; @@ -138,7 +144,7 @@ class PastixBase : internal::noncopyable public: - PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0) + PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_pastixdata(0), m_size(0) { init(); } @@ -148,6 +154,7 @@ class PastixBase : internal::noncopyable clean(); } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -161,19 +168,11 @@ class PastixBase : internal::noncopyable && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); return internal::solve_retval<PastixBase, Rhs>(*this, b.derived()); } +#endif template<typename Rhs,typename Dest> - bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; + bool _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const; - Derived& derived() - { - return *static_cast<Derived*>(this); - } - const Derived& derived() const - { - return *static_cast<const Derived*>(this); - } - /** Returns a reference to the integer vector IPARM of PaStiX parameters * to modify the default parameters. * The statistics related to the different phases of factorization and solve are saved here as well @@ -228,6 +227,7 @@ class PastixBase : internal::noncopyable return m_info; } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -241,6 +241,7 @@ class PastixBase : internal::noncopyable && "PastixBase::solve(): invalid number of rows of the right hand side matrix b"); return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived()); } +#endif // EIGEN_TEST_EVALUATORS protected: @@ -268,7 +269,6 @@ class PastixBase : internal::noncopyable int m_initisOk; int m_analysisIsOk; int m_factorizationIsOk; - bool m_isInitialized; mutable ComputationInfo m_info; mutable pastix_data_t *m_pastixdata; // Data structure for pastix mutable int m_comm; // The MPI communicator identifier @@ -393,7 +393,7 @@ void PastixBase<Derived>::factorize(ColSpMatrix& mat) /* Solve the system */ template<typename Base> template<typename Rhs,typename Dest> -bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const +bool PastixBase<Base>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const { eigen_assert(m_isInitialized && "The matrix should be factorized first"); EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0, @@ -694,6 +694,7 @@ class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> > } }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Rhs> @@ -705,7 +706,7 @@ struct solve_retval<PastixBase<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; @@ -723,6 +724,7 @@ struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs> }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen diff --git a/Eigen/src/PardisoSupport/PardisoSupport.h b/Eigen/src/PardisoSupport/PardisoSupport.h index b6571069e..e3b1c5bc0 100644 --- a/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/Eigen/src/PardisoSupport/PardisoSupport.h @@ -96,10 +96,17 @@ namespace internal } template<class Derived> -class PardisoImpl : internal::noncopyable +class PardisoImpl : public SparseSolveBase<PardisoImpl<Derived> { + protected: + typedef SparseSolveBase<PardisoImpl<Derived> Base; + using Base::derived; + using Base::m_isInitialized; + typedef internal::pardiso_traits<Derived> Traits; public: + using base::_solve_impl; + typedef typename Traits::MatrixType MatrixType; typedef typename Traits::Scalar Scalar; typedef typename Traits::RealScalar RealScalar; @@ -118,7 +125,7 @@ class PardisoImpl : internal::noncopyable eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type"); m_iparm.setZero(); m_msglvl = 0; // No output - m_initialized = false; + m_isInitialized = false; } ~PardisoImpl() @@ -136,7 +143,7 @@ class PardisoImpl : internal::noncopyable */ ComputationInfo info() const { - eigen_assert(m_initialized && "Decomposition is not initialized."); + eigen_assert(m_isInitialized && "Decomposition is not initialized."); return m_info; } @@ -166,6 +173,7 @@ class PardisoImpl : internal::noncopyable Derived& compute(const MatrixType& matrix); +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -174,7 +182,7 @@ class PardisoImpl : internal::noncopyable inline const internal::solve_retval<PardisoImpl, Rhs> solve(const MatrixBase<Rhs>& b) const { - eigen_assert(m_initialized && "Pardiso solver is not initialized."); + eigen_assert(m_isInitialized && "Pardiso solver is not initialized."); eigen_assert(rows()==b.rows() && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b"); return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived()); @@ -188,28 +196,20 @@ class PardisoImpl : internal::noncopyable inline const internal::sparse_solve_retval<PardisoImpl, Rhs> solve(const SparseMatrixBase<Rhs>& b) const { - eigen_assert(m_initialized && "Pardiso solver is not initialized."); + eigen_assert(m_isInitialized && "Pardiso solver is not initialized."); eigen_assert(rows()==b.rows() && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b"); return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived()); } - - Derived& derived() - { - return *static_cast<Derived*>(this); - } - const Derived& derived() const - { - return *static_cast<const Derived*>(this); - } +#endif template<typename BDerived, typename XDerived> - bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const; + bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const; protected: void pardisoRelease() { - if(m_initialized) // Factorization ran at least once + if(m_isInitialized) // Factorization ran at least once { internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0, m_iparm.data(), m_msglvl, 0, 0); @@ -270,7 +270,7 @@ class PardisoImpl : internal::noncopyable mutable SparseMatrixType m_matrix; ComputationInfo m_info; - bool m_initialized, m_analysisIsOk, m_factorizationIsOk; + bool m_analysisIsOk, m_factorizationIsOk; Index m_type, m_msglvl; mutable void *m_pt[64]; mutable ParameterType m_iparm; @@ -298,7 +298,7 @@ Derived& PardisoImpl<Derived>::compute(const MatrixType& a) manageErrorCode(error); m_analysisIsOk = true; m_factorizationIsOk = true; - m_initialized = true; + m_isInitialized = true; return derived(); } @@ -321,7 +321,7 @@ Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a) manageErrorCode(error); m_analysisIsOk = true; m_factorizationIsOk = false; - m_initialized = true; + m_isInitialized = true; return derived(); } @@ -345,7 +345,7 @@ Derived& PardisoImpl<Derived>::factorize(const MatrixType& a) template<class Base> template<typename BDerived,typename XDerived> -bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const +bool PardisoImpl<Base>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const { if(m_iparm[0] == 0) // Factorization was not computed return false; @@ -546,6 +546,7 @@ class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> > } }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _Derived, typename Rhs> @@ -557,7 +558,7 @@ struct solve_retval<PardisoImpl<_Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; @@ -575,6 +576,7 @@ struct sparse_solve_retval<PardisoImpl<Derived>, Rhs> }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 8808e6c0d..136e80ce9 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -350,7 +350,8 @@ void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) c dst.bottomRows(cols()-rank).setZero(); } #endif - + +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Rhs> @@ -366,6 +367,7 @@ struct solve_retval<HouseholderQR<_MatrixType>, Rhs> }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS /** Performs the QR factorization of the given matrix \a matrix. The result of * the factorization is stored into \c *this, and a reference to \c *this diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h index a2cc2a9e2..483aaef07 100644 --- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h +++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2012 Desire Nuentsa <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 @@ -54,8 +55,11 @@ namespace Eigen { * */ template<typename _MatrixType> -class SPQR +class SPQR : public SparseSolverBase<SPQR<_MatrixType> > { + protected: + typedef SparseSolverBase<SPQR<_MatrixType> > Base; + using Base::m_isInitialized; public: typedef typename _MatrixType::Scalar Scalar; typedef typename _MatrixType::RealScalar RealScalar; @@ -64,19 +68,13 @@ class SPQR typedef PermutationMatrix<Dynamic, Dynamic> PermutationType; public: SPQR() - : m_isInitialized(false), - m_ordering(SPQR_ORDERING_DEFAULT), - m_allow_tol(SPQR_DEFAULT_TOL), - m_tolerance (NumTraits<Scalar>::epsilon()) + : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()) { cholmod_l_start(&m_cc); } SPQR(const _MatrixType& matrix) - : m_isInitialized(false), - m_ordering(SPQR_ORDERING_DEFAULT), - m_allow_tol(SPQR_DEFAULT_TOL), - m_tolerance (NumTraits<Scalar>::epsilon()) + : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()) { cholmod_l_start(&m_cc); compute(matrix); @@ -127,10 +125,11 @@ class SPQR */ inline Index cols() const { return m_cR->ncol; } - /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. - * - * \sa compute() - */ +#ifndef EIGEN_TEST_EVALUATORS + /** \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<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const { @@ -139,9 +138,10 @@ class SPQR && "SPQR::solve(): invalid number of rows of the right hand side matrix B"); return internal::solve_retval<SPQR, Rhs>(*this, B.derived()); } +#endif // EIGEN_TEST_EVALUATORS template<typename Rhs, typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const { eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()"); eigen_assert(b.cols()==1 && "This method is for vectors only"); @@ -214,7 +214,6 @@ class SPQR return m_info; } protected: - bool m_isInitialized; bool m_analysisIsOk; bool m_factorizationIsOk; mutable bool m_isRUpToDate; @@ -293,6 +292,7 @@ struct SPQRMatrixQTransposeReturnType{ const SPQRType& m_spqr; }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Rhs> @@ -304,11 +304,12 @@ struct solve_retval<SPQR<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; } // end namespace internal +#endif }// End namespace Eigen #endif diff --git a/Eigen/src/SparseCholesky/SimplicialCholesky.h b/Eigen/src/SparseCholesky/SimplicialCholesky.h index 1abd31304..ac8cd29b0 100644 --- a/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -33,8 +33,11 @@ enum SimplicialCholeskyMode { * */ template<typename Derived> -class SimplicialCholeskyBase : internal::noncopyable +class SimplicialCholeskyBase : public SparseSolverBase<Derived> { + typedef SparseSolverBase<Derived> Base; + using Base::m_isInitialized; + public: typedef typename internal::traits<Derived>::MatrixType MatrixType; typedef typename internal::traits<Derived>::OrderingType OrderingType; @@ -46,14 +49,16 @@ class SimplicialCholeskyBase : internal::noncopyable typedef Matrix<Scalar,Dynamic,1> VectorType; public: + + using Base::derived; /** Default constructor */ SimplicialCholeskyBase() - : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) + : m_info(Success), m_shiftOffset(0), m_shiftScale(1) {} SimplicialCholeskyBase(const MatrixType& matrix) - : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) + : m_info(Success), m_shiftOffset(0), m_shiftScale(1) { derived().compute(matrix); } @@ -79,6 +84,7 @@ class SimplicialCholeskyBase : internal::noncopyable return m_info; } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -106,6 +112,7 @@ class SimplicialCholeskyBase : internal::noncopyable && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b"); return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); } +#endif // EIGEN_TEST_EVALUATORS /** \returns the permutation P * \sa permutationPinv() */ @@ -150,7 +157,11 @@ class SimplicialCholeskyBase : internal::noncopyable /** \internal */ template<typename Rhs,typename Dest> +#ifndef EIGEN_TEST_EVALUATORS void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const +#else + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const +#endif { eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); eigen_assert(m_matrix.rows()==b.rows()); @@ -176,6 +187,14 @@ class SimplicialCholeskyBase : internal::noncopyable dest = m_Pinv * dest; } +#ifdef EIGEN_TEST_EVALUATORS + template<typename Rhs,typename Dest> + void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const + { + internal::solve_sparse_through_dense_panels(derived(), b, dest); + } +#endif + #endif // EIGEN_PARSED_BY_DOXYGEN protected: @@ -226,7 +245,6 @@ class SimplicialCholeskyBase : internal::noncopyable }; mutable ComputationInfo m_info; - bool m_isInitialized; bool m_factorizationIsOk; bool m_analysisIsOk; @@ -560,7 +578,11 @@ public: /** \internal */ template<typename Rhs,typename Dest> +#ifndef EIGEN_TEST_EVALUATORS void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const +#else + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const +#endif { eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); eigen_assert(Base::m_matrix.rows()==b.rows()); @@ -596,6 +618,15 @@ public: dest = Base::m_Pinv * dest; } +#ifdef EIGEN_TEST_EVALUATORS + /** \internal */ + template<typename Rhs,typename Dest> + void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const + { + internal::solve_sparse_through_dense_panels(*this, b, dest); + } +#endif + Scalar determinant() const { if(m_LDLT) @@ -636,6 +667,7 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixTy ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); } +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename Derived, typename Rhs> @@ -665,6 +697,7 @@ struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs> }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen diff --git a/Eigen/src/SparseCore/SparseSolverBase.h b/Eigen/src/SparseCore/SparseSolverBase.h new file mode 100644 index 000000000..9a6ed1292 --- /dev/null +++ b/Eigen/src/SparseCore/SparseSolverBase.h @@ -0,0 +1,113 @@ +// 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_SPARSESOLVERBASE_H +#define EIGEN_SPARSESOLVERBASE_H + +namespace Eigen { + +namespace internal { + + /** \internal + * Helper functions to solve with a sparse right-hand-side and result. + * The rhs is decomposed into small vertical panels which are solved through dense temporaries. + */ +template<typename Decomposition, typename Rhs, typename Dest> +void solve_sparse_through_dense_panels(const Decomposition &dec, const Rhs& rhs, Dest &dest) +{ + EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); + typedef typename Dest::Scalar DestScalar; + // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. + static const int NbColsAtOnce = 4; + int rhsCols = rhs.cols(); + int size = rhs.rows(); + // the temporary matrices do not need more columns than NbColsAtOnce: + int tmpCols = (std::min)(rhsCols, NbColsAtOnce); + Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,tmpCols); + Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,tmpCols); + for(int k=0; k<rhsCols; k+=NbColsAtOnce) + { + int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); + tmp.leftCols(actualCols) = rhs.middleCols(k,actualCols); + tmpX.leftCols(actualCols) = dec.solve(tmp.leftCols(actualCols)); + dest.middleCols(k,actualCols) = tmpX.leftCols(actualCols).sparseView(); + } +} + +} // end namespace internal + +/** \class SparseSolverBase + * \ingroup SparseCore_Module + * \brief A base class for sparse solvers + * + * \tparam Derived the actual type of the solver. + * + */ +template<typename Derived> +class SparseSolverBase : internal::noncopyable +{ + public: + + /** Default constructor */ + SparseSolverBase() + : m_isInitialized(false) + {} + + ~SparseSolverBase() + {} + + Derived& derived() { return *static_cast<Derived*>(this); } + const Derived& derived() const { return *static_cast<const Derived*>(this); } + +#ifdef EIGEN_TEST_EVALUATORS + /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * \sa compute() + */ + template<typename Rhs> + inline const Solve<Derived, Rhs> + solve(const MatrixBase<Rhs>& b) const + { + 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 Solve<Derived, Rhs>(derived(), b.derived()); + } + + /** \returns an expression of the solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * \sa compute() + */ + template<typename Rhs> + inline const Solve<Derived, Rhs> + solve(const SparseMatrixBase<Rhs>& b) const + { + 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 Solve<Derived, Rhs>(derived(), b.derived()); + } +#endif + + +#ifndef EIGEN_PARSED_BY_DOXYGEN + /** \internal default implementation of solving with a sparse rhs */ + template<typename Rhs,typename Dest> + void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const + { + internal::solve_sparse_through_dense_panels(derived(), b.derived(), dest.derived()); + } +#endif // EIGEN_PARSED_BY_DOXYGEN + + protected: + + mutable bool m_isInitialized; +}; + +} // end namespace Eigen + +#endif // EIGEN_SPARSESOLVERBASE_H diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h index 7a9aeec2d..eb61fe3d9 100644 --- a/Eigen/src/SparseLU/SparseLU.h +++ b/Eigen/src/SparseLU/SparseLU.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> -// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2012-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 @@ -70,9 +70,14 @@ template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixURetu * \sa \ref OrderingMethods_Module */ template <typename _MatrixType, typename _OrderingType> -class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index> +class SparseLU : public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >, public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index> { + protected: + typedef SparseSolverBase<SparseLU<_MatrixType,_OrderingType> > APIBase; + using APIBase::m_isInitialized; public: + using APIBase::_solve_impl; + typedef _MatrixType MatrixType; typedef _OrderingType OrderingType; typedef typename MatrixType::Scalar Scalar; @@ -86,11 +91,11 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ typedef internal::SparseLUImpl<Scalar, Index> Base; public: - SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) + SparseLU():m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) { initperfvalues(); } - SparseLU(const MatrixType& matrix):m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) + SparseLU(const MatrixType& matrix):m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1) { initperfvalues(); compute(matrix); @@ -168,6 +173,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ m_diagpivotthresh = thresh; } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. * * \warning the destination matrix X in X = this->solve(B) must be colmun-major. @@ -195,6 +201,20 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ && "SparseLU::solve(): invalid number of rows of the right hand side matrix B"); return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived()); } +#else + +#ifdef EIGEN_PARSED_BY_DOXYGEN + /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. + * + * \warning the destination matrix X in X = this->solve(B) must be colmun-major. + * + * \sa compute() + */ + template<typename Rhs> + inline const Solve<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const; +#endif // EIGEN_PARSED_BY_DOXYGEN + +#endif // EIGEN_TEST_EVALUATORS /** \brief Reports whether previous computation was successful. * @@ -219,7 +239,7 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ } template<typename Rhs, typename Dest> - bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const + bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base) const { Dest& X(X_base.derived()); eigen_assert(m_factorizationIsOk && "The matrix should be factorized first"); @@ -332,7 +352,6 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ // Variables mutable ComputationInfo m_info; - bool m_isInitialized; bool m_factorizationIsOk; bool m_analysisIsOk; std::string m_lastError; @@ -728,6 +747,7 @@ struct SparseLUMatrixUReturnType : internal::no_assignment_operator const MatrixUType& m_mapU; }; +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Derived, typename Rhs> @@ -739,7 +759,7 @@ struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; @@ -756,6 +776,7 @@ struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs> } }; } // end namespace internal +#endif } // End namespace Eigen diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index 6be569533..8e946b045 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -62,9 +62,13 @@ namespace internal { * */ template<typename _MatrixType, typename _OrderingType> -class SparseQR +class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > { + protected: + typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base; + using Base::m_isInitialized; public: + using Base::_solve_impl; typedef _MatrixType MatrixType; typedef _OrderingType OrderingType; typedef typename MatrixType::Scalar Scalar; @@ -75,7 +79,7 @@ class SparseQR typedef Matrix<Scalar, Dynamic, 1> ScalarVector; typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; public: - SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) + SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { } /** Construct a QR factorization of the matrix \a mat. @@ -84,7 +88,7 @@ class SparseQR * * \sa compute() */ - SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) + SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { compute(mat); } @@ -162,7 +166,7 @@ class SparseQR /** \internal */ template<typename Rhs, typename Dest> - bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const + bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const { eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); @@ -186,7 +190,6 @@ class SparseQR m_info = Success; return true; } - /** Sets the threshold that is used to determine linearly dependent columns during the factorization. * @@ -199,6 +202,7 @@ class SparseQR m_threshold = threshold; } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. * * \sa compute() @@ -217,6 +221,26 @@ class SparseQR eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); return internal::sparse_solve_retval<SparseQR, Rhs>(*this, B.derived()); } +#else + /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. + * + * \sa compute() + */ + template<typename Rhs> + inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const + { + eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); + eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); + return Solve<SparseQR, Rhs>(*this, B.derived()); + } + template<typename Rhs> + inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const + { + eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); + eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); + return Solve<SparseQR, Rhs>(*this, B.derived()); + } +#endif // EIGEN_TEST_EVALUATORS /** \brief Reports whether previous computation was successful. * @@ -244,7 +268,6 @@ class SparseQR protected: - bool m_isInitialized; bool m_analysisIsok; bool m_factorizationIsok; mutable ComputationInfo m_info; @@ -554,6 +577,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) m_info = Success; } +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename OrderingType, typename Rhs> @@ -565,7 +589,7 @@ struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; template<typename _MatrixType, typename OrderingType, typename Rhs> @@ -581,6 +605,7 @@ struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs> } }; } // end namespace internal +#endif // EIGEN_TEST_EVALUATORS template <typename SparseQRType, typename Derived> struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> > diff --git a/Eigen/src/SuperLUSupport/SuperLUSupport.h b/Eigen/src/SuperLUSupport/SuperLUSupport.h index bcb355760..fcecd4fcf 100644 --- a/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2008-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 @@ -288,8 +288,12 @@ MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat) * \brief The base class for the direct and incomplete LU factorization of SuperLU */ template<typename _MatrixType, typename Derived> -class SuperLUBase : internal::noncopyable +class SuperLUBase : public SparseSolverBase<Derived> { + protected: + typedef SparseSolverBase<Derived> Base; + using Base::derived; + using Base::m_isInitialized; public: typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; @@ -309,9 +313,6 @@ class SuperLUBase : internal::noncopyable clearFactors(); } - Derived& derived() { return *static_cast<Derived*>(this); } - const Derived& derived() const { return *static_cast<const Derived*>(this); } - inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } @@ -336,6 +337,7 @@ class SuperLUBase : internal::noncopyable derived().factorize(matrix); } +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -361,7 +363,8 @@ class SuperLUBase : internal::noncopyable && "SuperLU::solve(): invalid number of rows of the right hand side matrix b"); return internal::sparse_solve_retval<SuperLUBase, Rhs>(*this, b.derived()); } - +#endif // EIGEN_TEST_EVALUATORS + /** Performs a symbolic decomposition on the sparcity of \a matrix. * * This function is particularly useful when solving for several problems having the same structure. @@ -453,7 +456,6 @@ class SuperLUBase : internal::noncopyable mutable char m_sluEqued; mutable ComputationInfo m_info; - bool m_isInitialized; int m_factorizationIsOk; int m_analysisIsOk; mutable bool m_extractedDataAreDirty; @@ -491,6 +493,7 @@ class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > typedef TriangularView<LUMatrixType, Upper> UMatrixType; public: + using Base::_solve_impl; SuperLU() : Base() { init(); } @@ -528,7 +531,7 @@ class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ template<typename Rhs,typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; #endif // EIGEN_PARSED_BY_DOXYGEN inline const LMatrixType& matrixL() const @@ -637,7 +640,7 @@ void SuperLU<MatrixType>::factorize(const MatrixType& a) template<typename MatrixType> template<typename Rhs,typename Dest> -void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const +void SuperLU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const { eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); @@ -828,6 +831,7 @@ class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > typedef typename Base::Index Index; public: + using Base::_solve_impl; SuperILU() : Base() { init(); } @@ -863,7 +867,7 @@ class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ template<typename Rhs,typename Dest> - void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; + void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const; #endif // EIGEN_PARSED_BY_DOXYGEN protected: @@ -948,7 +952,7 @@ void SuperILU<MatrixType>::factorize(const MatrixType& a) template<typename MatrixType> template<typename Rhs,typename Dest> -void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const +void SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const { eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()"); @@ -991,6 +995,7 @@ void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) } #endif +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Derived, typename Rhs> @@ -1002,7 +1007,7 @@ struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec().derived()._solve(rhs(),dst); + dec().derived()._solve_impl(rhs(),dst); } }; @@ -1020,7 +1025,7 @@ struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs> }; } // end namespace internal - +#endif } // end namespace Eigen #endif // EIGEN_SUPERLUSUPPORT_H diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index 3a48cecf7..845c8076a 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -121,9 +121,13 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N * \sa \ref TutorialSparseDirectSolvers */ template<typename _MatrixType> -class UmfPackLU : internal::noncopyable +class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> > { + protected: + typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base; + using Base::m_isInitialized; public: + using Base::_solve_impl; typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -197,7 +201,8 @@ class UmfPackLU : internal::noncopyable analyzePattern(matrix); factorize(matrix); } - + +#ifndef EIGEN_TEST_EVALUATORS /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * \sa compute() @@ -223,6 +228,7 @@ class UmfPackLU : internal::noncopyable && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b"); return internal::sparse_solve_retval<UmfPackLU, Rhs>(*this, b.derived()); } +#endif /** Performs a symbolic decomposition on the sparcity of \a matrix. * @@ -274,7 +280,7 @@ class UmfPackLU : internal::noncopyable #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal */ template<typename BDerived,typename XDerived> - bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; + bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const; #endif Scalar determinant() const; @@ -328,7 +334,6 @@ class UmfPackLU : internal::noncopyable void* m_symbolic; mutable ComputationInfo m_info; - bool m_isInitialized; int m_factorizationIsOk; int m_analysisIsOk; mutable bool m_extractedDataAreDirty; @@ -376,7 +381,7 @@ typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() cons template<typename MatrixType> template<typename BDerived,typename XDerived> -bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const +bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const { const int rhsCols = b.cols(); eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet"); @@ -396,7 +401,7 @@ bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDe return true; } - +#ifndef EIGEN_TEST_EVALUATORS namespace internal { template<typename _MatrixType, typename Rhs> @@ -408,7 +413,7 @@ struct solve_retval<UmfPackLU<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - dec()._solve(rhs(),dst); + dec()._solve_impl(rhs(),dst); } }; @@ -426,7 +431,7 @@ struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs> }; } // end namespace internal - +#endif } // end namespace Eigen #endif // EIGEN_UMFPACKSUPPORT_H |