src/Core/Solve.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_SOLVE_H
0011 #define EIGEN_SOLVE_H
0012
0013 namespace Eigen {
0014
0015 template<typename Decomposition, typename RhsType, typename StorageKind> class SolveImpl;
0016
0017 /** \class Solve
0018   * \ingroup Core_Module
0019   *
0020   * \brief Pseudo expression representing a solving operation
0021   *
0022   * \tparam Decomposition the type of the matrix or decomposition 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 // this solve_traits class permits to determine the evaluation type with respect to storage kind (Dense vs Sparse)
0032 template<typename Decomposition, typename RhsType,typename StorageKind> struct solve_traits;
0033
0034 template<typename Decomposition, typename RhsType>
0035 struct solve_traits<Decomposition,RhsType,Dense>
0036 {
0037   typedef typename make_proper_matrix_type<typename RhsType::Scalar,
0038                  Decomposition::ColsAtCompileTime,
0039                  RhsType::ColsAtCompileTime,
0040                  RhsType::PlainObject::Options,
0041                  Decomposition::MaxColsAtCompileTime,
0042                  RhsType::MaxColsAtCompileTime>::type PlainObject;
0043 };
0044
0045 template<typename Decomposition, typename RhsType>
0046 struct traits<Solve<Decomposition, RhsType> >
0047   : traits<typename solve_traits<Decomposition,RhsType,typename internal::traits<RhsType>::StorageKind>::PlainObject>
0048 {
0049   typedef typename solve_traits<Decomposition,RhsType,typename internal::traits<RhsType>::StorageKind>::PlainObject PlainObject;
0050   typedef typename promote_index_type<typename Decomposition::StorageIndex, typename RhsType::StorageIndex>::type StorageIndex;
0051   typedef traits<PlainObject> BaseTraits;
0052   enum {
0053     Flags = BaseTraits::Flags & RowMajorBit,
0054     CoeffReadCost = HugeCost
0055   };
0056 };
0057
0058 }
0059
0060
0061 template<typename Decomposition, typename RhsType>
0062 class Solve : public SolveImpl<Decomposition,RhsType,typename internal::traits<RhsType>::StorageKind>
0063 {
0064 public:
0065   typedef typename internal::traits<Solve>::PlainObject PlainObject;
0066   typedef typename internal::traits<Solve>::StorageIndex StorageIndex;
0067
0068   Solve(const Decomposition &dec, const RhsType &rhs)
0069     : m_dec(dec), m_rhs(rhs)
0070   {}
0071
0072   EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_dec.cols(); }
0073   EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_rhs.cols(); }
0074
0075   EIGEN_DEVICE_FUNC const Decomposition& dec() const { return m_dec; }
0076   EIGEN_DEVICE_FUNC const RhsType&       rhs() const { return m_rhs; }
0077
0078 protected:
0079   const Decomposition &m_dec;
0080   const RhsType       &m_rhs;
0081 };
0082
0083
0084 // Specialization of the Solve expression for dense results
0085 template<typename Decomposition, typename RhsType>
0086 class SolveImpl<Decomposition,RhsType,Dense>
0087   : public MatrixBase<Solve<Decomposition,RhsType> >
0088 {
0089   typedef Solve<Decomposition,RhsType> Derived;
0090
0091 public:
0092
0093   typedef MatrixBase<Solve<Decomposition,RhsType> > Base;
0094   EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
0095
0096 private:
0097
0098   Scalar coeff(Index row, Index col) const;
0099   Scalar coeff(Index i) const;
0100 };
0101
0102 // Generic API dispatcher
0103 template<typename Decomposition, typename RhsType, typename StorageKind>
0104 class SolveImpl : public internal::generic_xpr_base<Solve<Decomposition,RhsType>, MatrixXpr, StorageKind>::type
0105 {
0106   public:
0107     typedef typename internal::generic_xpr_base<Solve<Decomposition,RhsType>, MatrixXpr, StorageKind>::type Base;
0108 };
0109
0110 namespace internal {
0111
0112 // Evaluator of Solve -> eval into a temporary
0113 template<typename Decomposition, typename RhsType>
0114 struct evaluator<Solve<Decomposition,RhsType> >
0115   : public evaluator<typename Solve<Decomposition,RhsType>::PlainObject>
0116 {
0117   typedef Solve<Decomposition,RhsType> SolveType;
0118   typedef typename SolveType::PlainObject PlainObject;
0119   typedef evaluator<PlainObject> Base;
0120
0121   enum { Flags = Base::Flags | EvalBeforeNestingBit };
0122
0123   EIGEN_DEVICE_FUNC explicit evaluator(const SolveType& solve)
0124     : m_result(solve.rows(), solve.cols())
0125   {
0126     ::new (static_cast<Base*>(this)) Base(m_result);
0127     solve.dec()._solve_impl(solve.rhs(), m_result);
0128   }
0129
0130 protected:
0131   PlainObject m_result;
0132 };
0133
0134 // Specialization for "dst = dec.solve(rhs)"
0135 // NOTE we need to specialize it for Dense2Dense to avoid ambiguous specialization error and a Sparse2Sparse specialization must exist somewhere
0136 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
0137 struct Assignment<DstXprType, Solve<DecType,RhsType>, internal::assign_op<Scalar,Scalar>, Dense2Dense>
0138 {
0139   typedef Solve<DecType,RhsType> SrcXprType;
0140   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
0141   {
0142     Index dstRows = src.rows();
0143     Index dstCols = src.cols();
0144     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0145       dst.resize(dstRows, dstCols);
0146
0147     src.dec()._solve_impl(src.rhs(), dst);
0148   }
0149 };
0150
0151 // Specialization for "dst = dec.transpose().solve(rhs)"
0152 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
0153 struct Assignment<DstXprType, Solve<Transpose<const DecType>,RhsType>, internal::assign_op<Scalar,Scalar>, Dense2Dense>
0154 {
0155   typedef Solve<Transpose<const DecType>,RhsType> SrcXprType;
0156   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
0157   {
0158     Index dstRows = src.rows();
0159     Index dstCols = src.cols();
0160     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0161       dst.resize(dstRows, dstCols);
0162
0163     src.dec().nestedExpression().template _solve_impl_transposed<false>(src.rhs(), dst);
0164   }
0165 };
0166
0167 // Specialization for "dst = dec.adjoint().solve(rhs)"
0168 template<typename DstXprType, typename DecType, typename RhsType, typename Scalar>
0169 struct Assignment<DstXprType, Solve<CwiseUnaryOp<internal::scalar_conjugate_op<typename DecType::Scalar>, const Transpose<const DecType> >,RhsType>,
0170                   internal::assign_op<Scalar,Scalar>, Dense2Dense>
0171 {
0172   typedef Solve<CwiseUnaryOp<internal::scalar_conjugate_op<typename DecType::Scalar>, const Transpose<const DecType> >,RhsType> SrcXprType;
0173   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
0174   {
0175     Index dstRows = src.rows();
0176     Index dstCols = src.cols();
0177     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0178       dst.resize(dstRows, dstCols);
0179
0180     src.dec().nestedExpression().nestedExpression().template _solve_impl_transposed<true>(src.rhs(), dst);
0181   }
0182 };
0183
0184 } // end namespace internal
0185
0186 } // end namespace Eigen
0187
0188 #endif // EIGEN_SOLVE_H