$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2009 Guillaume Saupin <guillaume.saupin@cea.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SKYLINEPRODUCT_H 00011 #define EIGEN_SKYLINEPRODUCT_H 00012 00013 namespace Eigen { 00014 00015 template<typename Lhs, typename Rhs, int ProductMode> 00016 struct SkylineProductReturnType { 00017 typedef const typename internal::nested<Lhs, Rhs::RowsAtCompileTime>::type LhsNested; 00018 typedef const typename internal::nested<Rhs, Lhs::RowsAtCompileTime>::type RhsNested; 00019 00020 typedef SkylineProduct<LhsNested, RhsNested, ProductMode> Type; 00021 }; 00022 00023 template<typename LhsNested, typename RhsNested, int ProductMode> 00024 struct internal::traits<SkylineProduct<LhsNested, RhsNested, ProductMode> > { 00025 // clean the nested types: 00026 typedef typename internal::remove_all<LhsNested>::type _LhsNested; 00027 typedef typename internal::remove_all<RhsNested>::type _RhsNested; 00028 typedef typename _LhsNested::Scalar Scalar; 00029 00030 enum { 00031 LhsCoeffReadCost = _LhsNested::CoeffReadCost, 00032 RhsCoeffReadCost = _RhsNested::CoeffReadCost, 00033 LhsFlags = _LhsNested::Flags, 00034 RhsFlags = _RhsNested::Flags, 00035 00036 RowsAtCompileTime = _LhsNested::RowsAtCompileTime, 00037 ColsAtCompileTime = _RhsNested::ColsAtCompileTime, 00038 InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime), 00039 00040 MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime, 00041 MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime, 00042 00043 EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit), 00044 ResultIsSkyline = ProductMode == SkylineTimeSkylineProduct, 00045 00046 RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSkyline ? 0 : SkylineBit)), 00047 00048 Flags = (int(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits) 00049 | EvalBeforeAssigningBit 00050 | EvalBeforeNestingBit, 00051 00052 CoeffReadCost = Dynamic 00053 }; 00054 00055 typedef typename internal::conditional<ResultIsSkyline, 00056 SkylineMatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> >, 00057 MatrixBase<SkylineProduct<LhsNested, RhsNested, ProductMode> > >::type Base; 00058 }; 00059 00060 namespace internal { 00061 template<typename LhsNested, typename RhsNested, int ProductMode> 00062 class SkylineProduct : no_assignment_operator, 00063 public traits<SkylineProduct<LhsNested, RhsNested, ProductMode> >::Base { 00064 public: 00065 00066 EIGEN_GENERIC_PUBLIC_INTERFACE(SkylineProduct) 00067 00068 private: 00069 00070 typedef typename traits<SkylineProduct>::_LhsNested _LhsNested; 00071 typedef typename traits<SkylineProduct>::_RhsNested _RhsNested; 00072 00073 public: 00074 00075 template<typename Lhs, typename Rhs> 00076 EIGEN_STRONG_INLINE SkylineProduct(const Lhs& lhs, const Rhs& rhs) 00077 : m_lhs(lhs), m_rhs(rhs) { 00078 eigen_assert(lhs.cols() == rhs.rows()); 00079 00080 enum { 00081 ProductIsValid = _LhsNested::ColsAtCompileTime == Dynamic 00082 || _RhsNested::RowsAtCompileTime == Dynamic 00083 || int(_LhsNested::ColsAtCompileTime) == int(_RhsNested::RowsAtCompileTime), 00084 AreVectors = _LhsNested::IsVectorAtCompileTime && _RhsNested::IsVectorAtCompileTime, 00085 SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(_LhsNested, _RhsNested) 00086 }; 00087 // note to the lost user: 00088 // * for a dot product use: v1.dot(v2) 00089 // * for a coeff-wise product use: v1.cwise()*v2 00090 EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), 00091 INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) 00092 EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), 00093 INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) 00094 EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) 00095 } 00096 00097 EIGEN_STRONG_INLINE Index rows() const { 00098 return m_lhs.rows(); 00099 } 00100 00101 EIGEN_STRONG_INLINE Index cols() const { 00102 return m_rhs.cols(); 00103 } 00104 00105 EIGEN_STRONG_INLINE const _LhsNested& lhs() const { 00106 return m_lhs; 00107 } 00108 00109 EIGEN_STRONG_INLINE const _RhsNested& rhs() const { 00110 return m_rhs; 00111 } 00112 00113 protected: 00114 LhsNested m_lhs; 00115 RhsNested m_rhs; 00116 }; 00117 00118 // dense = skyline * dense 00119 // Note that here we force no inlining and separate the setZero() because GCC messes up otherwise 00120 00121 template<typename Lhs, typename Rhs, typename Dest> 00122 EIGEN_DONT_INLINE void skyline_row_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) { 00123 typedef typename remove_all<Lhs>::type _Lhs; 00124 typedef typename remove_all<Rhs>::type _Rhs; 00125 typedef typename traits<Lhs>::Scalar Scalar; 00126 00127 enum { 00128 LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit, 00129 LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit, 00130 ProcessFirstHalf = LhsIsSelfAdjoint 00131 && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0) 00132 || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor) 00133 || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)), 00134 ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf) 00135 }; 00136 00137 //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix. 00138 for (Index col = 0; col < rhs.cols(); col++) { 00139 for (Index row = 0; row < lhs.rows(); row++) { 00140 dst(row, col) = lhs.coeffDiag(row) * rhs(row, col); 00141 } 00142 } 00143 //Use matrix lower triangular part 00144 for (Index row = 0; row < lhs.rows(); row++) { 00145 typename _Lhs::InnerLowerIterator lIt(lhs, row); 00146 const Index stop = lIt.col() + lIt.size(); 00147 for (Index col = 0; col < rhs.cols(); col++) { 00148 00149 Index k = lIt.col(); 00150 Scalar tmp = 0; 00151 while (k < stop) { 00152 tmp += 00153 lIt.value() * 00154 rhs(k++, col); 00155 ++lIt; 00156 } 00157 dst(row, col) += tmp; 00158 lIt += -lIt.size(); 00159 } 00160 00161 } 00162 00163 //Use matrix upper triangular part 00164 for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) { 00165 typename _Lhs::InnerUpperIterator uIt(lhs, lhscol); 00166 const Index stop = uIt.size() + uIt.row(); 00167 for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) { 00168 00169 00170 const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol); 00171 Index k = uIt.row(); 00172 while (k < stop) { 00173 dst(k++, rhscol) += 00174 uIt.value() * 00175 rhsCoeff; 00176 ++uIt; 00177 } 00178 uIt += -uIt.size(); 00179 } 00180 } 00181 00182 } 00183 00184 template<typename Lhs, typename Rhs, typename Dest> 00185 EIGEN_DONT_INLINE void skyline_col_major_time_dense_product(const Lhs& lhs, const Rhs& rhs, Dest& dst) { 00186 typedef typename remove_all<Lhs>::type _Lhs; 00187 typedef typename remove_all<Rhs>::type _Rhs; 00188 typedef typename traits<Lhs>::Scalar Scalar; 00189 00190 enum { 00191 LhsIsRowMajor = (_Lhs::Flags & RowMajorBit) == RowMajorBit, 00192 LhsIsSelfAdjoint = (_Lhs::Flags & SelfAdjointBit) == SelfAdjointBit, 00193 ProcessFirstHalf = LhsIsSelfAdjoint 00194 && (((_Lhs::Flags & (UpperTriangularBit | LowerTriangularBit)) == 0) 00195 || ((_Lhs::Flags & UpperTriangularBit) && !LhsIsRowMajor) 00196 || ((_Lhs::Flags & LowerTriangularBit) && LhsIsRowMajor)), 00197 ProcessSecondHalf = LhsIsSelfAdjoint && (!ProcessFirstHalf) 00198 }; 00199 00200 //Use matrix diagonal part <- Improvement : use inner iterator on dense matrix. 00201 for (Index col = 0; col < rhs.cols(); col++) { 00202 for (Index row = 0; row < lhs.rows(); row++) { 00203 dst(row, col) = lhs.coeffDiag(row) * rhs(row, col); 00204 } 00205 } 00206 00207 //Use matrix upper triangular part 00208 for (Index row = 0; row < lhs.rows(); row++) { 00209 typename _Lhs::InnerUpperIterator uIt(lhs, row); 00210 const Index stop = uIt.col() + uIt.size(); 00211 for (Index col = 0; col < rhs.cols(); col++) { 00212 00213 Index k = uIt.col(); 00214 Scalar tmp = 0; 00215 while (k < stop) { 00216 tmp += 00217 uIt.value() * 00218 rhs(k++, col); 00219 ++uIt; 00220 } 00221 00222 00223 dst(row, col) += tmp; 00224 uIt += -uIt.size(); 00225 } 00226 } 00227 00228 //Use matrix lower triangular part 00229 for (Index lhscol = 0; lhscol < lhs.cols(); lhscol++) { 00230 typename _Lhs::InnerLowerIterator lIt(lhs, lhscol); 00231 const Index stop = lIt.size() + lIt.row(); 00232 for (Index rhscol = 0; rhscol < rhs.cols(); rhscol++) { 00233 00234 const Scalar rhsCoeff = rhs.coeff(lhscol, rhscol); 00235 Index k = lIt.row(); 00236 while (k < stop) { 00237 dst(k++, rhscol) += 00238 lIt.value() * 00239 rhsCoeff; 00240 ++lIt; 00241 } 00242 lIt += -lIt.size(); 00243 } 00244 } 00245 00246 } 00247 00248 template<typename Lhs, typename Rhs, typename ResultType, 00249 int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit> 00250 struct skyline_product_selector; 00251 00252 template<typename Lhs, typename Rhs, typename ResultType> 00253 struct skyline_product_selector<Lhs, Rhs, ResultType, RowMajor> { 00254 typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar; 00255 00256 static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) { 00257 skyline_row_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res); 00258 } 00259 }; 00260 00261 template<typename Lhs, typename Rhs, typename ResultType> 00262 struct skyline_product_selector<Lhs, Rhs, ResultType, ColMajor> { 00263 typedef typename traits<typename remove_all<Lhs>::type>::Scalar Scalar; 00264 00265 static void run(const Lhs& lhs, const Rhs& rhs, ResultType & res) { 00266 skyline_col_major_time_dense_product<Lhs, Rhs, ResultType > (lhs, rhs, res); 00267 } 00268 }; 00269 00270 } // end namespace internal 00271 00272 // template<typename Derived> 00273 // template<typename Lhs, typename Rhs > 00274 // Derived & MatrixBase<Derived>::lazyAssign(const SkylineProduct<Lhs, Rhs, SkylineTimeDenseProduct>& product) { 00275 // typedef typename internal::remove_all<Lhs>::type _Lhs; 00276 // internal::skyline_product_selector<typename internal::remove_all<Lhs>::type, 00277 // typename internal::remove_all<Rhs>::type, 00278 // Derived>::run(product.lhs(), product.rhs(), derived()); 00279 // 00280 // return derived(); 00281 // } 00282 00283 // skyline * dense 00284 00285 template<typename Derived> 00286 template<typename OtherDerived > 00287 EIGEN_STRONG_INLINE const typename SkylineProductReturnType<Derived, OtherDerived>::Type 00288 SkylineMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const { 00289 00290 return typename SkylineProductReturnType<Derived, OtherDerived>::Type(derived(), other.derived()); 00291 } 00292 00293 } // end namespace Eigen 00294 00295 #endif // EIGEN_SKYLINEPRODUCT_H