Squashed 'third_party/eigen/' content from commit 61d72f6
Change-Id: Iccc90fa0b55ab44037f018046d2fcffd90d9d025
git-subtree-dir: third_party/eigen
git-subtree-split: 61d72f6383cfa842868c53e30e087b0258177257
diff --git a/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
new file mode 100644
index 0000000..2036922
--- /dev/null
+++ b/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h
@@ -0,0 +1,254 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 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_ITERATIVE_SOLVER_BASE_H
+#define EIGEN_ITERATIVE_SOLVER_BASE_H
+
+namespace Eigen {
+
+/** \ingroup IterativeLinearSolvers_Module
+ * \brief Base class for linear iterative solvers
+ *
+ * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
+ */
+template< typename Derived>
+class IterativeSolverBase : internal::noncopyable
+{
+public:
+ typedef typename internal::traits<Derived>::MatrixType MatrixType;
+ typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::RealScalar RealScalar;
+
+public:
+
+ Derived& derived() { return *static_cast<Derived*>(this); }
+ const Derived& derived() const { return *static_cast<const Derived*>(this); }
+
+ /** Default constructor. */
+ IterativeSolverBase()
+ : mp_matrix(0)
+ {
+ init();
+ }
+
+ /** 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
+ * matrix A, or modify a copy of A.
+ */
+ IterativeSolverBase(const MatrixType& A)
+ {
+ init();
+ compute(A);
+ }
+
+ ~IterativeSolverBase() {}
+
+ /** Initializes the iterative solver for the sparcity pattern of the matrix \a A for further solving \c Ax=b problems.
+ *
+ * Currently, this function mostly call analyzePattern on the preconditioner. In the future
+ * we might, for instance, implement column reodering for faster matrix vector products.
+ */
+ Derived& analyzePattern(const MatrixType& A)
+ {
+ m_preconditioner.analyzePattern(A);
+ m_isInitialized = true;
+ m_analysisIsOk = true;
+ m_info = Success;
+ 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 call factorize on the preconditioner.
+ *
+ * \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
+ * matrix A, or modify a copy of A.
+ */
+ Derived& factorize(const MatrixType& A)
+ {
+ eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+ mp_matrix = &A;
+ m_preconditioner.factorize(A);
+ m_factorizationIsOk = true;
+ m_info = Success;
+ return derived();
+ }
+
+ /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems.
+ *
+ * Currently, this function mostly initialized/compute the preconditioner. In the future
+ * we might, for instance, implement column reodering for faster matrix vector products.
+ *
+ * \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
+ * matrix A, or modify a copy of A.
+ */
+ Derived& compute(const MatrixType& A)
+ {
+ mp_matrix = &A;
+ m_preconditioner.compute(A);
+ m_isInitialized = true;
+ m_analysisIsOk = true;
+ m_factorizationIsOk = true;
+ m_info = Success;
+ return derived();
+ }
+
+ /** \internal */
+ Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
+ /** \internal */
+ Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
+
+ /** \returns the tolerance threshold used by the stopping criteria */
+ RealScalar tolerance() const { return m_tolerance; }
+
+ /** Sets the tolerance threshold used by the stopping criteria */
+ Derived& setTolerance(const RealScalar& tolerance)
+ {
+ m_tolerance = tolerance;
+ return derived();
+ }
+
+ /** \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 */
+ int maxIterations() const
+ {
+ return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
+ }
+
+ /** Sets the max number of iterations */
+ Derived& setMaxIterations(int maxIters)
+ {
+ m_maxIterations = maxIters;
+ return derived();
+ }
+
+ /** \returns the number of iterations performed during the last solve */
+ int iterations() const
+ {
+ eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
+ return m_iterations;
+ }
+
+ /** \returns the tolerance error reached during the last solve */
+ RealScalar error() const
+ {
+ eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
+ 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.
+ *
+ * \sa compute()
+ */
+ template<typename Rhs>
+ inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
+ solve(const SparseMatrixBase<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::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
+ }
+
+ /** \returns Success if the iterations converged, and NoConvergence otherwise. */
+ ComputationInfo info() const
+ {
+ eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
+ return m_info;
+ }
+
+ /** \internal */
+ template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
+ void _solve_sparse(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);
+ }
+ }
+
+protected:
+ void init()
+ {
+ m_isInitialized = false;
+ m_analysisIsOk = false;
+ m_factorizationIsOk = false;
+ m_maxIterations = -1;
+ m_tolerance = NumTraits<Scalar>::epsilon();
+ }
+ const MatrixType* mp_matrix;
+ Preconditioner m_preconditioner;
+
+ int m_maxIterations;
+ RealScalar m_tolerance;
+
+ 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);
+ }
+};
+
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_ITERATIVE_SOLVER_BASE_H