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 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, 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 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 110 template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size, 111 bool Stop = LoopIndex==Size> 112 struct triangular_solver_unroller; 113 114 template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size> 115 struct 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 134 template<typename Lhs, typename Rhs, int Mode, int LoopIndex, int Size> 135 struct triangular_solver_unroller<Lhs,Rhs,Mode,LoopIndex,Size,true> { 136 static void run(const Lhs&, Rhs&) {} 137 }; 138 139 template<typename Lhs, typename Rhs, int Mode> 140 struct 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 145 template<typename Lhs, typename Rhs, int Mode> 146 struct 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 165 template<typename MatrixType, unsigned int Mode> 166 template<int Side, typename OtherDerived> 167 void 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 185 template<typename Derived, unsigned int Mode> 186 template<int Side, typename Other> 187 const internal::triangular_solve_retval<Side,TriangularView<Derived,Mode>,Other> 188 TriangularViewImpl<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 194 namespace internal { 195 196 197 template<int Side, typename TriangularType, typename Rhs> 198 struct traits<triangular_solve_retval<Side, TriangularType, Rhs> > 199 { 200 typedef typename internal::plain_matrix_type_column_major<Rhs>::type ReturnType; 201 }; 202 203 template<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