src/IterativeLinearSolvers/SolveWithGuess.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_SOLVEWITHGUESS_H
0011 #define EIGEN_SOLVEWITHGUESS_H
0012
0013 namespace RivetEigen {
0014
0015 template<typename Decomposition, typename RhsType, typename GuessType> class SolveWithGuess;
0016
0017 /** \class SolveWithGuess
0018   * \ingroup IterativeLinearSolvers_Module
0019   *
0020   * \brief Pseudo expression representing a solving operation
0021   *
0022   * \tparam Decomposition the type of the matrix or decomposion object
0023   * \tparam Rhstype the type of the right-hand side
0024   *
0025   * This class represents an expression of A.solve(B)
0026   * and most of the time this is the only way it is used.
0027   *
0028   */
0029 namespace internal {
0030
0031
0032 template<typename Decomposition, typename RhsType, typename GuessType>
0033 struct traits<SolveWithGuess<Decomposition, RhsType, GuessType> >
0034   : traits<Solve<Decomposition,RhsType> >
0035 {};
0036
0037 }
0038
0039
0040 template<typename Decomposition, typename RhsType, typename GuessType>
0041 class SolveWithGuess : public internal::generic_xpr_base<SolveWithGuess<Decomposition,RhsType,GuessType>, MatrixXpr, typename internal::traits<RhsType>::StorageKind>::type
0042 {
0043 public:
0044   typedef typename internal::traits<SolveWithGuess>::Scalar Scalar;
0045   typedef typename internal::traits<SolveWithGuess>::PlainObject PlainObject;
0046   typedef typename internal::generic_xpr_base<SolveWithGuess<Decomposition,RhsType,GuessType>, MatrixXpr, typename internal::traits<RhsType>::StorageKind>::type Base;
0047   typedef typename internal::ref_selector<SolveWithGuess>::type Nested;
0048
0049   SolveWithGuess(const Decomposition &dec, const RhsType &rhs, const GuessType &guess)
0050     : m_dec(dec), m_rhs(rhs), m_guess(guess)
0051   {}
0052
0053   EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0054   Index rows() const EIGEN_NOEXCEPT { return m_dec.cols(); }
0055   EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR
0056   Index cols() const EIGEN_NOEXCEPT { return m_rhs.cols(); }
0057
0058   EIGEN_DEVICE_FUNC const Decomposition& dec()   const { return m_dec; }
0059   EIGEN_DEVICE_FUNC const RhsType&       rhs()   const { return m_rhs; }
0060   EIGEN_DEVICE_FUNC const GuessType&     guess() const { return m_guess; }
0061
0062 protected:
0063   const Decomposition &m_dec;
0064   const RhsType       &m_rhs;
0065   const GuessType     &m_guess;
0066
0067 private:
0068   Scalar coeff(Index row, Index col) const;
0069   Scalar coeff(Index i) const;
0070 };
0071
0072 namespace internal {
0073
0074 // Evaluator of SolveWithGuess -> eval into a temporary
0075 template<typename Decomposition, typename RhsType, typename GuessType>
0076 struct evaluator<SolveWithGuess<Decomposition,RhsType, GuessType> >
0077   : public evaluator<typename SolveWithGuess<Decomposition,RhsType,GuessType>::PlainObject>
0078 {
0079   typedef SolveWithGuess<Decomposition,RhsType,GuessType> SolveType;
0080   typedef typename SolveType::PlainObject PlainObject;
0081   typedef evaluator<PlainObject> Base;
0082
0083   evaluator(const SolveType& solve)
0084     : m_result(solve.rows(), solve.cols())
0085   {
0086     ::new (static_cast<Base*>(this)) Base(m_result);
0087     m_result = solve.guess();
0088     solve.dec()._solve_with_guess_impl(solve.rhs(), m_result);
0089   }
0090
0091 protected:
0092   PlainObject m_result;
0093 };
0094
0095 // Specialization for "dst = dec.solveWithGuess(rhs)"
0096 // NOTE we need to specialize it for Dense2Dense to avoid ambiguous specialization error and a Sparse2Sparse specialization must exist somewhere
0097 template<typename DstXprType, typename DecType, typename RhsType, typename GuessType, typename Scalar>
0098 struct Assignment<DstXprType, SolveWithGuess<DecType,RhsType,GuessType>, internal::assign_op<Scalar,Scalar>, Dense2Dense>
0099 {
0100   typedef SolveWithGuess<DecType,RhsType,GuessType> SrcXprType;
0101   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
0102   {
0103     Index dstRows = src.rows();
0104     Index dstCols = src.cols();
0105     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0106       dst.resize(dstRows, dstCols);
0107
0108     dst = src.guess();
0109     src.dec()._solve_with_guess_impl(src.rhs(), dst/*, src.guess()*/);
0110   }
0111 };
0112
0113 } // end namespace internal
0114
0115 } // end namespace RivetEigen
0116
0117 #endif // EIGEN_SOLVEWITHGUESS_H