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@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 13namespace Eigen { 14 15namespace internal { 16 17// Forward declarations: 18// The following two routines are implemented in the products/TriangularSolver*.h files 19template<typename LhsScalar, typename RhsScalar, typename Index, int Side, int Mode, bool Conjugate, int StorageOrder> 20struct triangular_solve_vector; 21 22template <typename Scalar, typename Index, int Side, int Mode, bool Conjugate, int TriStorageOrder, int OtherStorageOrder> 23struct triangular_solve_matrix; 24 25// small helper struct extracting some traits on the underlying solver operation 26template<typename Lhs, typename Rhs, int Side> 27class 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 41template<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 > 47struct triangular_solver_selector; 48 49template<typename Lhs, typename Rhs, int Side, int Mode> 50struct 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, 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 81template<typename Lhs, typename Rhs, int Side, int Mode> 82struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,Dynamic> 83{ 84 typedef typename Rhs::Scalar Scalar; 85 typedef blas_traits<Lhs> LhsProductTraits; 86 typedef typename LhsProductTraits::DirectLinearAccessType ActualLhsType; 87 88 static void run(const Lhs& lhs, Rhs& rhs) 89 { 90 typename internal::add_const_on_value_type<ActualLhsType>::type actualLhs = LhsProductTraits::extract(lhs); 91 92 const Index size = lhs.rows(); 93 const Index othersize = Side==OnTheLeft? rhs.cols() : rhs.rows(); 94 95 typedef internal::gemm_blocking_space<(Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor,Scalar,Scalar, 96 Rhs::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime, Lhs::MaxRowsAtCompileTime,4> BlockingType; 97 98 BlockingType blocking(rhs.rows(), rhs.cols(), size, 1, false); 99 100 triangular_solve_matrix<Scalar,Index,Side,Mode,LhsProductTraits::NeedToConjugate,(int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor, 101 (Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor> 102 ::run(size, othersize, &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &rhs.coeffRef(0,0), rhs.outerStride(), blocking); 103 } 104}; 105 106/*************************************************************************** 107* meta-unrolling implementation 108***************************************************************************/ 109 110template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size, 111 bool Stop = LoopIndex==Size> 112struct triangular_solver_unroller; 113 114template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size> 115struct triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex,Size,false> { 116 enum { 117 IsLower = ((Mode&Lower)==Lower), 118 DiagIndex = IsLower ? LoopIndex : Size - LoopIndex - 1, 119 StartIndex = IsLower ? 0 : DiagIndex+1 120 }; 121 static void run(const Lhs& lhs, Rhs& rhs) 122 { 123 if (LoopIndex>0) 124 rhs.coeffRef(DiagIndex) -= lhs.row(DiagIndex).template segment<LoopIndex>(StartIndex).transpose() 125 .cwiseProduct(rhs.template segment<LoopIndex>(StartIndex)).sum(); 126 127 if(!(Mode & UnitDiag)) 128 rhs.coeffRef(DiagIndex) /= lhs.coeff(DiagIndex,DiagIndex); 129 130 triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex+1,Size>::run(lhs,rhs); 131 } 132}; 133 134template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size> 135struct triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex,Size,true> { 136 static void run(const Lhs&, Rhs&) {} 137}; 138 139template<typename Lhs, typename Rhs, int Mode> 140struct triangular_solver_selector<Lhs,Rhs,OnTheLeft,Mode,CompleteUnrolling,1> { 141 static void run(const Lhs& lhs, Rhs& rhs) 142 { triangular_solver_unroller<Lhs,Rhs,Mode,0,Rhs::SizeAtCompileTime>::run(lhs,rhs); } 143}; 144 145template<typename Lhs, typename Rhs, int Mode> 146struct triangular_solver_selector<Lhs,Rhs,OnTheRight,Mode,CompleteUnrolling,1> { 147 static void run(const Lhs& lhs, Rhs& rhs) 148 { 149 Transpose<const Lhs> trLhs(lhs); 150 Transpose<Rhs> trRhs(rhs); 151 152 triangular_solver_unroller<Transpose<const Lhs>,Transpose<Rhs>, 153 ((Mode&Upper)==Upper ? Lower : Upper) | (Mode&UnitDiag), 154 0,Rhs::SizeAtCompileTime>::run(trLhs,trRhs); 155 } 156}; 157 158} // end namespace internal 159 160/*************************************************************************** 161* TriangularView methods 162***************************************************************************/ 163 164#ifndef EIGEN_PARSED_BY_DOXYGEN 165template<typename MatrixType, unsigned int Mode> 166template<int Side, typename OtherDerived> 167void TriangularViewImpl<MatrixType,Mode,Dense>::solveInPlace(const MatrixBase<OtherDerived>& _other) const 168{ 169 OtherDerived& other = _other.const_cast_derived(); 170 eigen_assert( derived().cols() == derived().rows() && ((Side==OnTheLeft && derived().cols() == other.rows()) || (Side==OnTheRight && derived().cols() == other.cols())) ); 171 eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper|Lower))); 172 173 enum { copy = (internal::traits<OtherDerived>::Flags & RowMajorBit) && OtherDerived::IsVectorAtCompileTime && OtherDerived::SizeAtCompileTime!=1}; 174 typedef typename internal::conditional<copy, 175 typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::type OtherCopy; 176 OtherCopy otherCopy(other); 177 178 internal::triangular_solver_selector<MatrixType, typename internal::remove_reference<OtherCopy>::type, 179 Side, Mode>::run(derived().nestedExpression(), otherCopy); 180 181 if (copy) 182 other = otherCopy; 183} 184 185template<typename Derived, unsigned int Mode> 186template<int Side, typename Other> 187const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other> 188TriangularViewImpl<Derived,Mode,Dense>::solve(const MatrixBase<Other>& other) const 189{ 190 return internal::triangular_solve_retval<Side,TriangularViewType,Other>(derived(), other.derived()); 191} 192#endif 193 194namespace internal { 195 196 197template<int Side, typename TriangularType, typename Rhs> 198struct traits<triangular_solve_retval<Side, TriangularType, Rhs> > 199{ 200 typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType; 201}; 202 203template<int Side, typename TriangularType, typename Rhs> struct triangular_solve_retval 204 : public ReturnByValue<triangular_solve_retval<Side, TriangularType, Rhs> > 205{ 206 typedef typename remove_all<typename Rhs::Nested>::type RhsNestedCleaned; 207 typedef ReturnByValue<triangular_solve_retval> Base; 208 209 triangular_solve_retval(const TriangularType& tri, const Rhs& rhs) 210 : m_triangularMatrix(tri), m_rhs(rhs) 211 {} 212 213 inline Index rows() const { return m_rhs.rows(); } 214 inline Index cols() const { return m_rhs.cols(); } 215 216 template<typename Dest> inline void evalTo(Dest& dst) const 217 { 218 if(!is_same_dense(dst,m_rhs)) 219 dst = m_rhs; 220 m_triangularMatrix.template solveInPlace<Side>(dst); 221 } 222 223 protected: 224 const TriangularType& m_triangularMatrix; 225 typename Rhs::Nested m_rhs; 226}; 227 228} // namespace internal 229 230} // end namespace Eigen 231 232#endif // EIGEN_SOLVETRIANGULAR_H 233