Home | History | Annotate | Download | only in Core 
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_SOLVETRIANGULAR_H
     11 #define EIGEN_SOLVETRIANGULAR_H
     12 
     13 namespace Eigen {
     14 
     15 namespace internal {
     16 
     17 // Forward declarations:
     18 // The following two routines are implemented in the products/TriangularSolver*.h files
     19 template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder>
     20 struct triangular_solve_vector;
     21 
     22 template <typename Scalar, typename Index, int Side, int Mode, bool Conjugate, int TriStorageOrder, int OtherStorageOrder>
     23 struct triangular_solve_matrix;
     24 
     25 // small helper struct extracting some traits on the underlying solver operation
     26 template<typename Lhs, typename Rhs, int Side>
     27 class trsolve_traits
     28 {
     29   private:
     30     enum {
     31       RhsIsVectorAtCompileTime = (Side==OnTheLeft ? Rhs::ColsAtCompileTime : Rhs::RowsAtCompileTime)==1
     32     };
     33   public:
     34     enum {
     35       Unrolling   = (RhsIsVectorAtCompileTime && Rhs::SizeAtCompileTime != Dynamic && Rhs::SizeAtCompileTime <= 8)
     36                   ? CompleteUnrolling : NoUnrolling,
     37       RhsVectors  = RhsIsVectorAtCompileTime ? 1 : Dynamic
     38     };
     39 };
     40 
     41 template<typename Lhs, typename Rhs,
     42   int Side, // can be OnTheLeft/OnTheRight
     43   int Mode, // can be Upper/Lower | UnitDiag
     44   int Unrolling = trsolve_traits<Lhs,Rhs,Side>::Unrolling,
     45   int RhsVectors = trsolve_traits<Lhs,Rhs,Side>::RhsVectors
     46   >
     47 struct triangular_solver_selector;
     48 
     49 template<typename Lhs, typename Rhs, int Side, int Mode>
     50 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,1>
     51 {
     52   typedef typename Lhs::Scalar LhsScalar;
     53   typedef typename Rhs::Scalar RhsScalar;
     54   typedef blas_traits<Lhs> LhsProductTraits;
     55   typedef typename LhsProductTraits::ExtractType ActualLhsType;
     56   typedef Map<Matrix<RhsScalar,Dynamic,1>, Aligned> MappedRhs;
     57   static void run(const Lhs& lhs, Rhs& rhs)
     58   {
     59     ActualLhsType actualLhs = LhsProductTraits::extract(lhs);
     60 
     61     // FIXME find a way to allow an inner stride if packet_traits<Scalar>::size==1
     62 
     63     bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1;
     64 
     65     ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhs,rhs.size(),
     66                                                   (useRhsDirectly ? rhs.data() : 0));
     67 
     68     if(!useRhsDirectly)
     69       MappedRhs(actualRhs,rhs.size()) = rhs;
     70 
     71     triangular_solve_vector<LhsScalar, RhsScalar, typename Lhs::Index, Side, Mode, LhsProductTraits::NeedToConjugate,
     72                             (int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor>
     73       ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs);
     74 
     75     if(!useRhsDirectly)
     76       rhs = MappedRhs(actualRhs, rhs.size());
     77   }
     78 };
     79 
     80 // the rhs is a matrix
     81 template<typename Lhs, typename Rhs, int Side, int Mode>
     82 struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,Dynamic>
     83 {
     84   typedef typename Rhs::Scalar Scalar;
     85   typedef typename Rhs::Index Index;
     86   typedef blas_traits<Lhs> LhsProductTraits;
     87   typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType;
     88 
     89   static void run(const Lhs& lhs, Rhs& rhs)
     90   {
     91     typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs);
     92 
     93     const Index size = lhs.rows();
     94     const Index othersize = Side==OnTheLeft? rhs.cols() : rhs.rows();
     95 
     96     typedef internal::gemm_blocking_space<(Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor,Scalar,Scalar,
     97               Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType;
     98 
     99     BlockingType blocking(rhs.rows(), rhs.cols(), size);
    100 
    101     triangular_solve_matrix<Scalar,Index,Side,Mode,LhsProductTraits::NeedToConjugate,(int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor,
    102                                (Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor>
    103       ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.outerStride(), blocking);
    104   }
    105 };
    106 
    107 /***************************************************************************
    108 * meta-unrolling implementation
    109 ***************************************************************************/
    110 
    111 template<typename Lhs, typename Rhs, int Mode, int Index, int Size,
    112          bool Stop = Index==Size>
    113 struct triangular_solver_unroller;
    114 
    115 template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
    116 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> {
    117   enum {
    118     IsLower = ((Mode&Lower)==Lower),
    119     I = IsLower ? Index : Size - Index - 1,
    120     S = IsLower ? 0     : I+1
    121   };
    122   static void run(const Lhs& lhs, Rhs& rhs)
    123   {
    124     if (Index>0)
    125       rhs.coeffRef(I) -= lhs.row(I).template segment<Index>(S).transpose()
    126                          .cwiseProduct(rhs.template segment<Index>(S)).sum();
    127 
    128     if(!(Mode & UnitDiag))
    129       rhs.coeffRef(I) /= lhs.coeff(I,I);
    130 
    131     triangular_solver_unroller<Lhs,Rhs,Mode,Index+1,Size>::run(lhs,rhs);
    132   }
    133 };
    134 
    135 template<typename Lhs, typename Rhs, int Mode, int Index, int Size>
    136 struct triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,true> {
    137   static void run(const Lhs&, Rhs&) {}
    138 };
    139 
    140 template<typename Lhs, typename Rhs, int Mode>
    141 struct triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,CompleteUnrolling,1> {
    142   static void run(const Lhs& lhs, Rhs& rhs)
    143   { triangular_solver_unroller<Lhs,Rhs,Mode,0,Rhs::SizeAtCompileTime>::run(lhs,rhs); }
    144 };
    145 
    146 template<typename Lhs, typename Rhs, int Mode>
    147 struct triangular_solver_selector<Lhs,Rhs,OnTheRight,Mode,CompleteUnrolling,1> {
    148   static void run(const Lhs& lhs, Rhs& rhs)
    149   {
    150     Transpose<const Lhs> trLhs(lhs);
    151     Transpose<Rhs> trRhs(rhs);
    152 
    153     triangular_solver_unroller<Transpose<const Lhs>,Transpose<Rhs>,
    154                               ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag),
    155                               0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs);
    156   }
    157 };
    158 
    159 } // end namespace internal
    160 
    161 /***************************************************************************
    162 * TriangularView methods
    163 ***************************************************************************/
    164 
    165 /** "in-place" version of TriangularView::solve() where the result is written in \a other
    166   *
    167   * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
    168   * This function will const_cast it, so constness isn't honored here.
    169   *
    170   * See TriangularView:solve() for the details.
    171   */
    172 template<typename MatrixType, unsigned int Mode>
    173 template<int Side, typename OtherDerived>
    174 void TriangularView<MatrixType,Mode>::solveInPlace(const MatrixBase<OtherDerived>& _other) const
    175 {
    176   OtherDerived& other = _other.const_cast_derived();
    177   eigen_assert( cols() == rows() && ((Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols())) );
    178   eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower)));
    179 
    180   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit  && OtherDerived::IsVectorAtCompileTime };
    181   typedef typename internal::conditional<copy,
    182     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy;
    183   OtherCopy otherCopy(other);
    184 
    185   internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type,
    186     Side, Mode>::run(nestedExpression(), otherCopy);
    187 
    188   if (copy)
    189     other = otherCopy;
    190 }
    191 
    192 /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular.
    193   *
    194   * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other if
    195   * \a Side==OnTheLeft (the default), or the right-inverse-multiply  \a other * inverse(\c *this) if
    196   * \a Side==OnTheRight.
    197   *
    198   * The matrix \c *this must be triangular and invertible (i.e., all the coefficients of the
    199   * diagonal must be non zero). It works as a forward (resp. backward) substitution if \c *this
    200   * is an upper (resp. lower) triangular matrix.
    201   *
    202   * Example: \include MatrixBase_marked.cpp
    203   * Output: \verbinclude MatrixBase_marked.out
    204   *
    205   * This function returns an expression of the inverse-multiply and can works in-place if it is assigned
    206   * to the same matrix or vector \a other.
    207   *
    208   * For users coming from BLAS, this function (and more specifically solveInPlace()) offer
    209   * all the operations supported by the \c *TRSV and \c *TRSM BLAS routines.
    210   *
    211   * \sa TriangularView::solveInPlace()
    212   */
    213 template<typename Derived, unsigned int Mode>
    214 template<int Side, typename Other>
    215 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other>
    216 TriangularView<Derived,Mode>::solve(const MatrixBase<Other>& other) const
    217 {
    218   return internal::triangular_solve_retval<Side,TriangularView,Other>(*this, other.derived());
    219 }
    220 
    221 namespace internal {
    222 
    223 
    224 template<int Side, typename TriangularType, typename Rhs>
    225 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> >
    226 {
    227   typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType;
    228 };
    229 
    230 template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval
    231  : public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> >
    232 {
    233   typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned;
    234   typedef ReturnByValue<triangular_solve_retval> Base;
    235   typedef typename Base::Index Index;
    236 
    237   triangular_solve_retval(const TriangularType& tri, const Rhs& rhs)
    238     : m_triangularMatrix(tri), m_rhs(rhs)
    239   {}
    240 
    241   inline Index rows() const { return m_rhs.rows(); }
    242   inline Index cols() const { return m_rhs.cols(); }
    243 
    244   template<typename Dest> inline void evalTo(Dest& dst) const
    245   {
    246     if(!(is_same<RhsNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_rhs)))
    247       dst = m_rhs;
    248     m_triangularMatrix.template solveInPlace<Side>(dst);
    249   }
    250 
    251   protected:
    252     const TriangularType& m_triangularMatrix;
    253     typename Rhs::Nested m_rhs;
    254 };
    255 
    256 } // namespace internal
    257 
    258 } // end namespace Eigen
    259 
    260 #endif // EIGEN_SOLVETRIANGULAR_H
    261