$treeview $search $mathjax
Eigen
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) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> 00005 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_GENERAL_PRODUCT_H 00012 #define EIGEN_GENERAL_PRODUCT_H 00013 00014 namespace Eigen { 00015 00035 template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value> 00036 class GeneralProduct; 00037 00038 enum { 00039 Large = 2, 00040 Small = 3 00041 }; 00042 00043 namespace internal { 00044 00045 template<int Rows, int Cols, int Depth> struct product_type_selector; 00046 00047 template<int Size, int MaxSize> struct product_size_category 00048 { 00049 enum { is_large = MaxSize == Dynamic || 00050 Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD, 00051 value = is_large ? Large 00052 : Size == 1 ? 1 00053 : Small 00054 }; 00055 }; 00056 00057 template<typename Lhs, typename Rhs> struct product_type 00058 { 00059 typedef typename remove_all<Lhs>::type _Lhs; 00060 typedef typename remove_all<Rhs>::type _Rhs; 00061 enum { 00062 MaxRows = _Lhs::MaxRowsAtCompileTime, 00063 Rows = _Lhs::RowsAtCompileTime, 00064 MaxCols = _Rhs::MaxColsAtCompileTime, 00065 Cols = _Rhs::ColsAtCompileTime, 00066 MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::MaxColsAtCompileTime, 00067 _Rhs::MaxRowsAtCompileTime), 00068 Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime, 00069 _Rhs::RowsAtCompileTime), 00070 LargeThreshold = EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD 00071 }; 00072 00073 // the splitting into different lines of code here, introducing the _select enums and the typedef below, 00074 // is to work around an internal compiler error with gcc 4.1 and 4.2. 00075 private: 00076 enum { 00077 rows_select = product_size_category<Rows,MaxRows>::value, 00078 cols_select = product_size_category<Cols,MaxCols>::value, 00079 depth_select = product_size_category<Depth,MaxDepth>::value 00080 }; 00081 typedef product_type_selector<rows_select, cols_select, depth_select> selector; 00082 00083 public: 00084 enum { 00085 value = selector::ret 00086 }; 00087 #ifdef EIGEN_DEBUG_PRODUCT 00088 static void debug() 00089 { 00090 EIGEN_DEBUG_VAR(Rows); 00091 EIGEN_DEBUG_VAR(Cols); 00092 EIGEN_DEBUG_VAR(Depth); 00093 EIGEN_DEBUG_VAR(rows_select); 00094 EIGEN_DEBUG_VAR(cols_select); 00095 EIGEN_DEBUG_VAR(depth_select); 00096 EIGEN_DEBUG_VAR(value); 00097 } 00098 #endif 00099 }; 00100 00101 00102 /* The following allows to select the kind of product at compile time 00103 * based on the three dimensions of the product. 00104 * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */ 00105 // FIXME I'm not sure the current mapping is the ideal one. 00106 template<int M, int N> struct product_type_selector<M,N,1> { enum { ret = OuterProduct }; }; 00107 template<int Depth> struct product_type_selector<1, 1, Depth> { enum { ret = InnerProduct }; }; 00108 template<> struct product_type_selector<1, 1, 1> { enum { ret = InnerProduct }; }; 00109 template<> struct product_type_selector<Small,1, Small> { enum { ret = CoeffBasedProductMode }; }; 00110 template<> struct product_type_selector<1, Small,Small> { enum { ret = CoeffBasedProductMode }; }; 00111 template<> struct product_type_selector<Small,Small,Small> { enum { ret = CoeffBasedProductMode }; }; 00112 template<> struct product_type_selector<Small, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; }; 00113 template<> struct product_type_selector<Small, Large, 1> { enum { ret = LazyCoeffBasedProductMode }; }; 00114 template<> struct product_type_selector<Large, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; }; 00115 template<> struct product_type_selector<1, Large,Small> { enum { ret = CoeffBasedProductMode }; }; 00116 template<> struct product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; }; 00117 template<> struct product_type_selector<1, Small,Large> { enum { ret = CoeffBasedProductMode }; }; 00118 template<> struct product_type_selector<Large,1, Small> { enum { ret = CoeffBasedProductMode }; }; 00119 template<> struct product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; }; 00120 template<> struct product_type_selector<Small,1, Large> { enum { ret = CoeffBasedProductMode }; }; 00121 template<> struct product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; }; 00122 template<> struct product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; }; 00123 template<> struct product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; }; 00124 template<> struct product_type_selector<Large,Large,Large> { enum { ret = GemmProduct }; }; 00125 template<> struct product_type_selector<Large,Small,Small> { enum { ret = GemmProduct }; }; 00126 template<> struct product_type_selector<Small,Large,Small> { enum { ret = GemmProduct }; }; 00127 template<> struct product_type_selector<Large,Large,Small> { enum { ret = GemmProduct }; }; 00128 00129 } // end namespace internal 00130 00148 template<typename Lhs, typename Rhs, int ProductType> 00149 struct ProductReturnType 00150 { 00151 // TODO use the nested type to reduce instanciations ???? 00152 // typedef typename internal::nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested; 00153 // typedef typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested; 00154 00155 typedef GeneralProduct<Lhs/*Nested*/, Rhs/*Nested*/, ProductType> Type; 00156 }; 00157 00158 template<typename Lhs, typename Rhs> 00159 struct ProductReturnType<Lhs,Rhs,CoeffBasedProductMode> 00160 { 00161 typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested; 00162 typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested; 00163 typedef CoeffBasedProduct<LhsNested, RhsNested, EvalBeforeAssigningBit | EvalBeforeNestingBit> Type; 00164 }; 00165 00166 template<typename Lhs, typename Rhs> 00167 struct ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode> 00168 { 00169 typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested; 00170 typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested; 00171 typedef CoeffBasedProduct<LhsNested, RhsNested, NestByRefBit> Type; 00172 }; 00173 00174 // this is a workaround for sun CC 00175 template<typename Lhs, typename Rhs> 00176 struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode> 00177 {}; 00178 00179 /*********************************************************************** 00180 * Implementation of Inner Vector Vector Product 00181 ***********************************************************************/ 00182 00183 // FIXME : maybe the "inner product" could return a Scalar 00184 // instead of a 1x1 matrix ?? 00185 // Pro: more natural for the user 00186 // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix 00187 // product ends up to a row-vector times col-vector product... To tackle this use 00188 // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x); 00189 00190 namespace internal { 00191 00192 template<typename Lhs, typename Rhs> 00193 struct traits<GeneralProduct<Lhs,Rhs,InnerProduct> > 00194 : traits<Matrix<typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> > 00195 {}; 00196 00197 } 00198 00199 template<typename Lhs, typename Rhs> 00200 class GeneralProduct<Lhs, Rhs, InnerProduct> 00201 : internal::no_assignment_operator, 00202 public Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> 00203 { 00204 typedef Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> Base; 00205 public: 00206 GeneralProduct(const Lhs& lhs, const Rhs& rhs) 00207 { 00208 EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value), 00209 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 00210 00211 Base::coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum(); 00212 } 00213 00215 operator const typename Base::Scalar() const { 00216 return Base::coeff(0,0); 00217 } 00218 }; 00219 00220 /*********************************************************************** 00221 * Implementation of Outer Vector Vector Product 00222 ***********************************************************************/ 00223 00224 namespace internal { 00225 00226 // Column major 00227 template<typename ProductType, typename Dest, typename Func> 00228 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const false_type&) 00229 { 00230 typedef typename Dest::Index Index; 00231 // FIXME make sure lhs is sequentially stored 00232 // FIXME not very good if rhs is real and lhs complex while alpha is real too 00233 const Index cols = dest.cols(); 00234 for (Index j=0; j<cols; ++j) 00235 func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs()); 00236 } 00237 00238 // Row major 00239 template<typename ProductType, typename Dest, typename Func> 00240 EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const true_type&) { 00241 typedef typename Dest::Index Index; 00242 // FIXME make sure rhs is sequentially stored 00243 // FIXME not very good if lhs is real and rhs complex while alpha is real too 00244 const Index rows = dest.rows(); 00245 for (Index i=0; i<rows; ++i) 00246 func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs()); 00247 } 00248 00249 template<typename Lhs, typename Rhs> 00250 struct traits<GeneralProduct<Lhs,Rhs,OuterProduct> > 00251 : traits<ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> > 00252 {}; 00253 00254 } 00255 00256 template<typename Lhs, typename Rhs> 00257 class GeneralProduct<Lhs, Rhs, OuterProduct> 00258 : public ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> 00259 { 00260 template<typename T> struct IsRowMajor : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {}; 00261 00262 public: 00263 EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) 00264 00265 GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) 00266 { 00267 EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value), 00268 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 00269 } 00270 00271 struct set { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } }; 00272 struct add { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } }; 00273 struct sub { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } }; 00274 struct adds { 00275 Scalar m_scale; 00276 adds(const Scalar& s) : m_scale(s) {} 00277 template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { 00278 dst.const_cast_derived() += m_scale * src; 00279 } 00280 }; 00281 00282 template<typename Dest> 00283 inline void evalTo(Dest& dest) const { 00284 internal::outer_product_selector_run(*this, dest, set(), IsRowMajor<Dest>()); 00285 } 00286 00287 template<typename Dest> 00288 inline void addTo(Dest& dest) const { 00289 internal::outer_product_selector_run(*this, dest, add(), IsRowMajor<Dest>()); 00290 } 00291 00292 template<typename Dest> 00293 inline void subTo(Dest& dest) const { 00294 internal::outer_product_selector_run(*this, dest, sub(), IsRowMajor<Dest>()); 00295 } 00296 00297 template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const 00298 { 00299 internal::outer_product_selector_run(*this, dest, adds(alpha), IsRowMajor<Dest>()); 00300 } 00301 }; 00302 00303 /*********************************************************************** 00304 * Implementation of General Matrix Vector Product 00305 ***********************************************************************/ 00306 00307 /* According to the shape/flags of the matrix we have to distinghish 3 different cases: 00308 * 1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine 00309 * 2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine 00310 * 3 - all other cases are handled using a simple loop along the outer-storage direction. 00311 * Therefore we need a lower level meta selector. 00312 * Furthermore, if the matrix is the rhs, then the product has to be transposed. 00313 */ 00314 namespace internal { 00315 00316 template<typename Lhs, typename Rhs> 00317 struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> > 00318 : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> > 00319 {}; 00320 00321 template<int Side, int StorageOrder, bool BlasCompatible> 00322 struct gemv_selector; 00323 00324 } // end namespace internal 00325 00326 template<typename Lhs, typename Rhs> 00327 class GeneralProduct<Lhs, Rhs, GemvProduct> 00328 : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> 00329 { 00330 public: 00331 EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) 00332 00333 typedef typename Lhs::Scalar LhsScalar; 00334 typedef typename Rhs::Scalar RhsScalar; 00335 00336 GeneralProduct(const Lhs& a_lhs, const Rhs& a_rhs) : Base(a_lhs,a_rhs) 00337 { 00338 // EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::Scalar, typename Rhs::Scalar>::value), 00339 // YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 00340 } 00341 00342 enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight }; 00343 typedef typename internal::conditional<int(Side)==OnTheRight,_LhsNested,_RhsNested>::type MatrixType; 00344 00345 template<typename Dest> void scaleAndAddTo(Dest& dst, const Scalar& alpha) const 00346 { 00347 eigen_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols()); 00348 internal::gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor, 00349 bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha); 00350 } 00351 }; 00352 00353 namespace internal { 00354 00355 // The vector is on the left => transposition 00356 template<int StorageOrder, bool BlasCompatible> 00357 struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible> 00358 { 00359 template<typename ProductType, typename Dest> 00360 static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) 00361 { 00362 Transpose<Dest> destT(dest); 00363 enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor }; 00364 gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible> 00365 ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct> 00366 (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha); 00367 } 00368 }; 00369 00370 template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if; 00371 00372 template<typename Scalar,int Size,int MaxSize> 00373 struct gemv_static_vector_if<Scalar,Size,MaxSize,false> 00374 { 00375 EIGEN_STRONG_INLINE Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; } 00376 }; 00377 00378 template<typename Scalar,int Size> 00379 struct gemv_static_vector_if<Scalar,Size,Dynamic,true> 00380 { 00381 EIGEN_STRONG_INLINE Scalar* data() { return 0; } 00382 }; 00383 00384 template<typename Scalar,int Size,int MaxSize> 00385 struct gemv_static_vector_if<Scalar,Size,MaxSize,true> 00386 { 00387 #if EIGEN_ALIGN_STATICALLY 00388 internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data; 00389 EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; } 00390 #else 00391 // Some architectures cannot align on the stack, 00392 // => let's manually enforce alignment by allocating more data and return the address of the first aligned element. 00393 enum { 00394 ForceAlignment = internal::packet_traits<Scalar>::Vectorizable, 00395 PacketSize = internal::packet_traits<Scalar>::size 00396 }; 00397 internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data; 00398 EIGEN_STRONG_INLINE Scalar* data() { 00399 return ForceAlignment 00400 ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16) 00401 : m_data.array; 00402 } 00403 #endif 00404 }; 00405 00406 template<> struct gemv_selector<OnTheRight,ColMajor,true> 00407 { 00408 template<typename ProductType, typename Dest> 00409 static inline void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) 00410 { 00411 typedef typename ProductType::Index Index; 00412 typedef typename ProductType::LhsScalar LhsScalar; 00413 typedef typename ProductType::RhsScalar RhsScalar; 00414 typedef typename ProductType::Scalar ResScalar; 00415 typedef typename ProductType::RealScalar RealScalar; 00416 typedef typename ProductType::ActualLhsType ActualLhsType; 00417 typedef typename ProductType::ActualRhsType ActualRhsType; 00418 typedef typename ProductType::LhsBlasTraits LhsBlasTraits; 00419 typedef typename ProductType::RhsBlasTraits RhsBlasTraits; 00420 typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest; 00421 00422 ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs()); 00423 ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs()); 00424 00425 ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) 00426 * RhsBlasTraits::extractScalarFactor(prod.rhs()); 00427 00428 enum { 00429 // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1 00430 // on, the other hand it is good for the cache to pack the vector anyways... 00431 EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1, 00432 ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex), 00433 MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal 00434 }; 00435 00436 gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest; 00437 00438 bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0)); 00439 bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible; 00440 00441 RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha); 00442 00443 ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(), 00444 evalToDest ? dest.data() : static_dest.data()); 00445 00446 if(!evalToDest) 00447 { 00448 #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN 00449 int size = dest.size(); 00450 EIGEN_DENSE_STORAGE_CTOR_PLUGIN 00451 #endif 00452 if(!alphaIsCompatible) 00453 { 00454 MappedDest(actualDestPtr, dest.size()).setZero(); 00455 compatibleAlpha = RhsScalar(1); 00456 } 00457 else 00458 MappedDest(actualDestPtr, dest.size()) = dest; 00459 } 00460 00461 general_matrix_vector_product 00462 <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run( 00463 actualLhs.rows(), actualLhs.cols(), 00464 actualLhs.data(), actualLhs.outerStride(), 00465 actualRhs.data(), actualRhs.innerStride(), 00466 actualDestPtr, 1, 00467 compatibleAlpha); 00468 00469 if (!evalToDest) 00470 { 00471 if(!alphaIsCompatible) 00472 dest += actualAlpha * MappedDest(actualDestPtr, dest.size()); 00473 else 00474 dest = MappedDest(actualDestPtr, dest.size()); 00475 } 00476 } 00477 }; 00478 00479 template<> struct gemv_selector<OnTheRight,RowMajor,true> 00480 { 00481 template<typename ProductType, typename Dest> 00482 static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) 00483 { 00484 typedef typename ProductType::LhsScalar LhsScalar; 00485 typedef typename ProductType::RhsScalar RhsScalar; 00486 typedef typename ProductType::Scalar ResScalar; 00487 typedef typename ProductType::Index Index; 00488 typedef typename ProductType::ActualLhsType ActualLhsType; 00489 typedef typename ProductType::ActualRhsType ActualRhsType; 00490 typedef typename ProductType::_ActualRhsType _ActualRhsType; 00491 typedef typename ProductType::LhsBlasTraits LhsBlasTraits; 00492 typedef typename ProductType::RhsBlasTraits RhsBlasTraits; 00493 00494 typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs()); 00495 typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs()); 00496 00497 ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) 00498 * RhsBlasTraits::extractScalarFactor(prod.rhs()); 00499 00500 enum { 00501 // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1 00502 // on, the other hand it is good for the cache to pack the vector anyways... 00503 DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1 00504 }; 00505 00506 gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs; 00507 00508 ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(), 00509 DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data()); 00510 00511 if(!DirectlyUseRhs) 00512 { 00513 #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN 00514 int size = actualRhs.size(); 00515 EIGEN_DENSE_STORAGE_CTOR_PLUGIN 00516 #endif 00517 Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs; 00518 } 00519 00520 general_matrix_vector_product 00521 <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run( 00522 actualLhs.rows(), actualLhs.cols(), 00523 actualLhs.data(), actualLhs.outerStride(), 00524 actualRhsPtr, 1, 00525 dest.data(), dest.innerStride(), 00526 actualAlpha); 00527 } 00528 }; 00529 00530 template<> struct gemv_selector<OnTheRight,ColMajor,false> 00531 { 00532 template<typename ProductType, typename Dest> 00533 static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) 00534 { 00535 typedef typename Dest::Index Index; 00536 // TODO makes sure dest is sequentially stored in memory, otherwise use a temp 00537 const Index size = prod.rhs().rows(); 00538 for(Index k=0; k<size; ++k) 00539 dest += (alpha*prod.rhs().coeff(k)) * prod.lhs().col(k); 00540 } 00541 }; 00542 00543 template<> struct gemv_selector<OnTheRight,RowMajor,false> 00544 { 00545 template<typename ProductType, typename Dest> 00546 static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) 00547 { 00548 typedef typename Dest::Index Index; 00549 // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp 00550 const Index rows = prod.rows(); 00551 for(Index i=0; i<rows; ++i) 00552 dest.coeffRef(i) += alpha * (prod.lhs().row(i).cwiseProduct(prod.rhs().transpose())).sum(); 00553 } 00554 }; 00555 00556 } // end namespace internal 00557 00558 /*************************************************************************** 00559 * Implementation of matrix base methods 00560 ***************************************************************************/ 00561 00568 template<typename Derived> 00569 template<typename OtherDerived> 00570 inline const typename ProductReturnType<Derived, OtherDerived>::Type 00571 MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const 00572 { 00573 // A note regarding the function declaration: In MSVC, this function will sometimes 00574 // not be inlined since DenseStorage is an unwindable object for dynamic 00575 // matrices and product types are holding a member to store the result. 00576 // Thus it does not help tagging this function with EIGEN_STRONG_INLINE. 00577 enum { 00578 ProductIsValid = Derived::ColsAtCompileTime==Dynamic 00579 || OtherDerived::RowsAtCompileTime==Dynamic 00580 || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), 00581 AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, 00582 SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) 00583 }; 00584 // note to the lost user: 00585 // * for a dot product use: v1.dot(v2) 00586 // * for a coeff-wise product use: v1.cwiseProduct(v2) 00587 EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), 00588 INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) 00589 EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), 00590 INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) 00591 EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) 00592 #ifdef EIGEN_DEBUG_PRODUCT 00593 internal::product_type<Derived,OtherDerived>::debug(); 00594 #endif 00595 return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); 00596 } 00597 00609 template<typename Derived> 00610 template<typename OtherDerived> 00611 const typename LazyProductReturnType<Derived,OtherDerived>::Type 00612 MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &other) const 00613 { 00614 enum { 00615 ProductIsValid = Derived::ColsAtCompileTime==Dynamic 00616 || OtherDerived::RowsAtCompileTime==Dynamic 00617 || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), 00618 AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, 00619 SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) 00620 }; 00621 // note to the lost user: 00622 // * for a dot product use: v1.dot(v2) 00623 // * for a coeff-wise product use: v1.cwiseProduct(v2) 00624 EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), 00625 INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) 00626 EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), 00627 INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) 00628 EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) 00629 00630 return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); 00631 } 00632 00633 } // end namespace Eigen 00634 00635 #endif // EIGEN_PRODUCT_H