Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 1 | // This file is part of Eigen, a lightweight C++ template library |
| 2 | // for linear algebra. |
| 3 | // |
| 4 | // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> |
| 5 | // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> |
| 6 | // |
| 7 | // This Source Code Form is subject to the terms of the Mozilla |
| 8 | // Public License v. 2.0. If a copy of the MPL was not distributed |
| 9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| 10 | |
| 11 | #ifndef EIGEN_GENERAL_PRODUCT_H |
| 12 | #define EIGEN_GENERAL_PRODUCT_H |
| 13 | |
| 14 | namespace Eigen { |
| 15 | |
| 16 | /** \class GeneralProduct |
| 17 | * \ingroup Core_Module |
| 18 | * |
| 19 | * \brief Expression of the product of two general matrices or vectors |
| 20 | * |
| 21 | * \param LhsNested the type used to store the left-hand side |
| 22 | * \param RhsNested the type used to store the right-hand side |
| 23 | * \param ProductMode the type of the product |
| 24 | * |
| 25 | * This class represents an expression of the product of two general matrices. |
| 26 | * We call a general matrix, a dense matrix with full storage. For instance, |
| 27 | * This excludes triangular, selfadjoint, and sparse matrices. |
| 28 | * It is the return type of the operator* between general matrices. Its template |
| 29 | * arguments are determined automatically by ProductReturnType. Therefore, |
| 30 | * GeneralProduct should never be used direclty. To determine the result type of a |
| 31 | * function which involves a matrix product, use ProductReturnType::Type. |
| 32 | * |
| 33 | * \sa ProductReturnType, MatrixBase::operator*(const MatrixBase<OtherDerived>&) |
| 34 | */ |
| 35 | template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value> |
| 36 | class GeneralProduct; |
| 37 | |
| 38 | enum { |
| 39 | Large = 2, |
| 40 | Small = 3 |
| 41 | }; |
| 42 | |
| 43 | namespace internal { |
| 44 | |
| 45 | template<int Rows, int Cols, int Depth> struct product_type_selector; |
| 46 | |
| 47 | template<int Size, int MaxSize> struct product_size_category |
| 48 | { |
| 49 | enum { is_large = MaxSize == Dynamic || |
| 50 | Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD, |
| 51 | value = is_large ? Large |
| 52 | : Size == 1 ? 1 |
| 53 | : Small |
| 54 | }; |
| 55 | }; |
| 56 | |
| 57 | template<typename Lhs, typename Rhs> struct product_type |
| 58 | { |
| 59 | typedef typename remove_all<Lhs>::type _Lhs; |
| 60 | typedef typename remove_all<Rhs>::type _Rhs; |
| 61 | enum { |
| 62 | MaxRows = _Lhs::MaxRowsAtCompileTime, |
| 63 | Rows = _Lhs::RowsAtCompileTime, |
| 64 | MaxCols = _Rhs::MaxColsAtCompileTime, |
| 65 | Cols = _Rhs::ColsAtCompileTime, |
| 66 | MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::MaxColsAtCompileTime, |
| 67 | _Rhs::MaxRowsAtCompileTime), |
| 68 | Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime, |
| 69 | _Rhs::RowsAtCompileTime), |
| 70 | LargeThreshold = EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD |
| 71 | }; |
| 72 | |
| 73 | // the splitting into different lines of code here, introducing the _select enums and the typedef below, |
| 74 | // is to work around an internal compiler error with gcc 4.1 and 4.2. |
| 75 | private: |
| 76 | enum { |
| 77 | rows_select = product_size_category<Rows,MaxRows>::value, |
| 78 | cols_select = product_size_category<Cols,MaxCols>::value, |
| 79 | depth_select = product_size_category<Depth,MaxDepth>::value |
| 80 | }; |
| 81 | typedef product_type_selector<rows_select, cols_select, depth_select> selector; |
| 82 | |
| 83 | public: |
| 84 | enum { |
| 85 | value = selector::ret |
| 86 | }; |
| 87 | #ifdef EIGEN_DEBUG_PRODUCT |
| 88 | static void debug() |
| 89 | { |
| 90 | EIGEN_DEBUG_VAR(Rows); |
| 91 | EIGEN_DEBUG_VAR(Cols); |
| 92 | EIGEN_DEBUG_VAR(Depth); |
| 93 | EIGEN_DEBUG_VAR(rows_select); |
| 94 | EIGEN_DEBUG_VAR(cols_select); |
| 95 | EIGEN_DEBUG_VAR(depth_select); |
| 96 | EIGEN_DEBUG_VAR(value); |
| 97 | } |
| 98 | #endif |
| 99 | }; |
| 100 | |
| 101 | |
| 102 | /* The following allows to select the kind of product at compile time |
| 103 | * based on the three dimensions of the product. |
| 104 | * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */ |
| 105 | // FIXME I'm not sure the current mapping is the ideal one. |
| 106 | template<int M, int N> struct product_type_selector<M,N,1> { enum { ret = OuterProduct }; }; |
| 107 | template<int Depth> struct product_type_selector<1, 1, Depth> { enum { ret = InnerProduct }; }; |
| 108 | template<> struct product_type_selector<1, 1, 1> { enum { ret = InnerProduct }; }; |
| 109 | template<> struct product_type_selector<Small,1, Small> { enum { ret = CoeffBasedProductMode }; }; |
| 110 | template<> struct product_type_selector<1, Small,Small> { enum { ret = CoeffBasedProductMode }; }; |
| 111 | template<> struct product_type_selector<Small,Small,Small> { enum { ret = CoeffBasedProductMode }; }; |
| 112 | template<> struct product_type_selector<Small, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; }; |
| 113 | template<> struct product_type_selector<Small, Large, 1> { enum { ret = LazyCoeffBasedProductMode }; }; |
| 114 | template<> struct product_type_selector<Large, Small, 1> { enum { ret = LazyCoeffBasedProductMode }; }; |
| 115 | template<> struct product_type_selector<1, Large,Small> { enum { ret = CoeffBasedProductMode }; }; |
| 116 | template<> struct product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; }; |
| 117 | template<> struct product_type_selector<1, Small,Large> { enum { ret = CoeffBasedProductMode }; }; |
| 118 | template<> struct product_type_selector<Large,1, Small> { enum { ret = CoeffBasedProductMode }; }; |
| 119 | template<> struct product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; }; |
| 120 | template<> struct product_type_selector<Small,1, Large> { enum { ret = CoeffBasedProductMode }; }; |
| 121 | template<> struct product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; }; |
| 122 | template<> struct product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; }; |
| 123 | template<> struct product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; }; |
| 124 | template<> struct product_type_selector<Large,Large,Large> { enum { ret = GemmProduct }; }; |
| 125 | template<> struct product_type_selector<Large,Small,Small> { enum { ret = GemmProduct }; }; |
| 126 | template<> struct product_type_selector<Small,Large,Small> { enum { ret = GemmProduct }; }; |
| 127 | template<> struct product_type_selector<Large,Large,Small> { enum { ret = GemmProduct }; }; |
| 128 | |
| 129 | } // end namespace internal |
| 130 | |
| 131 | /** \class ProductReturnType |
| 132 | * \ingroup Core_Module |
| 133 | * |
| 134 | * \brief Helper class to get the correct and optimized returned type of operator* |
| 135 | * |
| 136 | * \param Lhs the type of the left-hand side |
| 137 | * \param Rhs the type of the right-hand side |
| 138 | * \param ProductMode the type of the product (determined automatically by internal::product_mode) |
| 139 | * |
| 140 | * This class defines the typename Type representing the optimized product expression |
| 141 | * between two matrix expressions. In practice, using ProductReturnType<Lhs,Rhs>::Type |
| 142 | * is the recommended way to define the result type of a function returning an expression |
| 143 | * which involve a matrix product. The class Product should never be |
| 144 | * used directly. |
| 145 | * |
| 146 | * \sa class Product, MatrixBase::operator*(const MatrixBase<OtherDerived>&) |
| 147 | */ |
| 148 | template<typename Lhs, typename Rhs, int ProductType> |
| 149 | struct ProductReturnType |
| 150 | { |
| 151 | // TODO use the nested type to reduce instanciations ???? |
| 152 | // typedef typename internal::nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested; |
| 153 | // typedef typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested; |
| 154 | |
| 155 | typedef GeneralProduct<Lhs/*Nested*/, Rhs/*Nested*/, ProductType> Type; |
| 156 | }; |
| 157 | |
| 158 | template<typename Lhs, typename Rhs> |
| 159 | struct ProductReturnType<Lhs,Rhs,CoeffBasedProductMode> |
| 160 | { |
| 161 | typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested; |
| 162 | typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested; |
| 163 | typedef CoeffBasedProduct<LhsNested, RhsNested, EvalBeforeAssigningBit | EvalBeforeNestingBit> Type; |
| 164 | }; |
| 165 | |
| 166 | template<typename Lhs, typename Rhs> |
| 167 | struct ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode> |
| 168 | { |
| 169 | typedef typename internal::nested<Lhs, Rhs::ColsAtCompileTime, typename internal::plain_matrix_type<Lhs>::type >::type LhsNested; |
| 170 | typedef typename internal::nested<Rhs, Lhs::RowsAtCompileTime, typename internal::plain_matrix_type<Rhs>::type >::type RhsNested; |
| 171 | typedef CoeffBasedProduct<LhsNested, RhsNested, NestByRefBit> Type; |
| 172 | }; |
| 173 | |
| 174 | // this is a workaround for sun CC |
| 175 | template<typename Lhs, typename Rhs> |
| 176 | struct LazyProductReturnType : public ProductReturnType<Lhs,Rhs,LazyCoeffBasedProductMode> |
| 177 | {}; |
| 178 | |
| 179 | /*********************************************************************** |
| 180 | * Implementation of Inner Vector Vector Product |
| 181 | ***********************************************************************/ |
| 182 | |
| 183 | // FIXME : maybe the "inner product" could return a Scalar |
| 184 | // instead of a 1x1 matrix ?? |
| 185 | // Pro: more natural for the user |
| 186 | // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix |
| 187 | // product ends up to a row-vector times col-vector product... To tackle this use |
| 188 | // case, we could have a specialization for Block<MatrixType,1,1> with: operator=(Scalar x); |
| 189 | |
| 190 | namespace internal { |
| 191 | |
| 192 | template<typename Lhs, typename Rhs> |
| 193 | struct traits<GeneralProduct<Lhs,Rhs,InnerProduct> > |
| 194 | : traits<Matrix<typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> > |
| 195 | {}; |
| 196 | |
| 197 | } |
| 198 | |
| 199 | template<typename Lhs, typename Rhs> |
| 200 | class GeneralProduct<Lhs, Rhs, InnerProduct> |
| 201 | : internal::no_assignment_operator, |
| 202 | public Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> |
| 203 | { |
| 204 | typedef Matrix<typename internal::scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType,1,1> Base; |
| 205 | public: |
| 206 | GeneralProduct(const Lhs& lhs, const Rhs& rhs) |
| 207 | { |
| 208 | EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value), |
| 209 | YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) |
| 210 | |
| 211 | Base::coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum(); |
| 212 | } |
| 213 | |
| 214 | /** Convertion to scalar */ |
| 215 | operator const typename Base::Scalar() const { |
| 216 | return Base::coeff(0,0); |
| 217 | } |
| 218 | }; |
| 219 | |
| 220 | /*********************************************************************** |
| 221 | * Implementation of Outer Vector Vector Product |
| 222 | ***********************************************************************/ |
| 223 | |
| 224 | namespace internal { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 225 | |
| 226 | // Column major |
| 227 | template<typename ProductType, typename Dest, typename Func> |
| 228 | EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const false_type&) |
| 229 | { |
| 230 | typedef typename Dest::Index Index; |
| 231 | // FIXME make sure lhs is sequentially stored |
| 232 | // FIXME not very good if rhs is real and lhs complex while alpha is real too |
| 233 | const Index cols = dest.cols(); |
| 234 | for (Index j=0; j<cols; ++j) |
| 235 | func(dest.col(j), prod.rhs().coeff(j) * prod.lhs()); |
| 236 | } |
| 237 | |
| 238 | // Row major |
| 239 | template<typename ProductType, typename Dest, typename Func> |
| 240 | EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& dest, const Func& func, const true_type&) { |
| 241 | typedef typename Dest::Index Index; |
| 242 | // FIXME make sure rhs is sequentially stored |
| 243 | // FIXME not very good if lhs is real and rhs complex while alpha is real too |
| 244 | const Index rows = dest.rows(); |
| 245 | for (Index i=0; i<rows; ++i) |
| 246 | func(dest.row(i), prod.lhs().coeff(i) * prod.rhs()); |
| 247 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 248 | |
| 249 | template<typename Lhs, typename Rhs> |
| 250 | struct traits<GeneralProduct<Lhs,Rhs,OuterProduct> > |
| 251 | : traits<ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> > |
| 252 | {}; |
| 253 | |
| 254 | } |
| 255 | |
| 256 | template<typename Lhs, typename Rhs> |
| 257 | class GeneralProduct<Lhs, Rhs, OuterProduct> |
| 258 | : public ProductBase<GeneralProduct<Lhs,Rhs,OuterProduct>, Lhs, Rhs> |
| 259 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 260 | template<typename T> struct IsRowMajor : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {}; |
| 261 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 262 | public: |
| 263 | EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) |
| 264 | |
| 265 | GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) |
| 266 | { |
| 267 | EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value), |
| 268 | YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) |
| 269 | } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 270 | |
| 271 | struct set { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() = src; } }; |
| 272 | struct add { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } }; |
| 273 | struct sub { template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } }; |
| 274 | struct adds { |
| 275 | Scalar m_scale; |
| 276 | adds(const Scalar& s) : m_scale(s) {} |
| 277 | template<typename Dst, typename Src> void operator()(const Dst& dst, const Src& src) const { |
| 278 | dst.const_cast_derived() += m_scale * src; |
| 279 | } |
| 280 | }; |
| 281 | |
| 282 | template<typename Dest> |
| 283 | inline void evalTo(Dest& dest) const { |
| 284 | internal::outer_product_selector_run(*this, dest, set(), IsRowMajor<Dest>()); |
| 285 | } |
| 286 | |
| 287 | template<typename Dest> |
| 288 | inline void addTo(Dest& dest) const { |
| 289 | internal::outer_product_selector_run(*this, dest, add(), IsRowMajor<Dest>()); |
| 290 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 291 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 292 | template<typename Dest> |
| 293 | inline void subTo(Dest& dest) const { |
| 294 | internal::outer_product_selector_run(*this, dest, sub(), IsRowMajor<Dest>()); |
| 295 | } |
| 296 | |
| 297 | template<typename Dest> void scaleAndAddTo(Dest& dest, const Scalar& alpha) const |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 298 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 299 | internal::outer_product_selector_run(*this, dest, adds(alpha), IsRowMajor<Dest>()); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 300 | } |
| 301 | }; |
| 302 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 303 | /*********************************************************************** |
| 304 | * Implementation of General Matrix Vector Product |
| 305 | ***********************************************************************/ |
| 306 | |
| 307 | /* According to the shape/flags of the matrix we have to distinghish 3 different cases: |
| 308 | * 1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine |
| 309 | * 2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine |
| 310 | * 3 - all other cases are handled using a simple loop along the outer-storage direction. |
| 311 | * Therefore we need a lower level meta selector. |
| 312 | * Furthermore, if the matrix is the rhs, then the product has to be transposed. |
| 313 | */ |
| 314 | namespace internal { |
| 315 | |
| 316 | template<typename Lhs, typename Rhs> |
| 317 | struct traits<GeneralProduct<Lhs,Rhs,GemvProduct> > |
| 318 | : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> > |
| 319 | {}; |
| 320 | |
| 321 | template<int Side, int StorageOrder, bool BlasCompatible> |
| 322 | struct gemv_selector; |
| 323 | |
| 324 | } // end namespace internal |
| 325 | |
| 326 | template<typename Lhs, typename Rhs> |
| 327 | class GeneralProduct<Lhs, Rhs, GemvProduct> |
| 328 | : public ProductBase<GeneralProduct<Lhs,Rhs,GemvProduct>, Lhs, Rhs> |
| 329 | { |
| 330 | public: |
| 331 | EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) |
| 332 | |
| 333 | typedef typename Lhs::Scalar LhsScalar; |
| 334 | typedef typename Rhs::Scalar RhsScalar; |
| 335 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 336 | GeneralProduct(const Lhs& a_lhs, const Rhs& a_rhs) : Base(a_lhs,a_rhs) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 337 | { |
| 338 | // EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::Scalar, typename Rhs::Scalar>::value), |
| 339 | // YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) |
| 340 | } |
| 341 | |
| 342 | enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight }; |
| 343 | typedef typename internal::conditional<int(Side)==OnTheRight,_LhsNested,_RhsNested>::type MatrixType; |
| 344 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 345 | template<typename Dest> void scaleAndAddTo(Dest& dst, const Scalar& alpha) const |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 346 | { |
| 347 | eigen_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols()); |
| 348 | internal::gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor, |
| 349 | bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)>::run(*this, dst, alpha); |
| 350 | } |
| 351 | }; |
| 352 | |
| 353 | namespace internal { |
| 354 | |
| 355 | // The vector is on the left => transposition |
| 356 | template<int StorageOrder, bool BlasCompatible> |
| 357 | struct gemv_selector<OnTheLeft,StorageOrder,BlasCompatible> |
| 358 | { |
| 359 | template<typename ProductType, typename Dest> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 360 | static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 361 | { |
| 362 | Transpose<Dest> destT(dest); |
| 363 | enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor }; |
| 364 | gemv_selector<OnTheRight,OtherStorageOrder,BlasCompatible> |
| 365 | ::run(GeneralProduct<Transpose<const typename ProductType::_RhsNested>,Transpose<const typename ProductType::_LhsNested>, GemvProduct> |
| 366 | (prod.rhs().transpose(), prod.lhs().transpose()), destT, alpha); |
| 367 | } |
| 368 | }; |
| 369 | |
| 370 | template<typename Scalar,int Size,int MaxSize,bool Cond> struct gemv_static_vector_if; |
| 371 | |
| 372 | template<typename Scalar,int Size,int MaxSize> |
| 373 | struct gemv_static_vector_if<Scalar,Size,MaxSize,false> |
| 374 | { |
| 375 | EIGEN_STRONG_INLINE Scalar* data() { eigen_internal_assert(false && "should never be called"); return 0; } |
| 376 | }; |
| 377 | |
| 378 | template<typename Scalar,int Size> |
| 379 | struct gemv_static_vector_if<Scalar,Size,Dynamic,true> |
| 380 | { |
| 381 | EIGEN_STRONG_INLINE Scalar* data() { return 0; } |
| 382 | }; |
| 383 | |
| 384 | template<typename Scalar,int Size,int MaxSize> |
| 385 | struct gemv_static_vector_if<Scalar,Size,MaxSize,true> |
| 386 | { |
| 387 | #if EIGEN_ALIGN_STATICALLY |
| 388 | internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data; |
| 389 | EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; } |
| 390 | #else |
| 391 | // Some architectures cannot align on the stack, |
| 392 | // => let's manually enforce alignment by allocating more data and return the address of the first aligned element. |
| 393 | enum { |
| 394 | ForceAlignment = internal::packet_traits<Scalar>::Vectorizable, |
| 395 | PacketSize = internal::packet_traits<Scalar>::size |
| 396 | }; |
| 397 | internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data; |
| 398 | EIGEN_STRONG_INLINE Scalar* data() { |
| 399 | return ForceAlignment |
| 400 | ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16) |
| 401 | : m_data.array; |
| 402 | } |
| 403 | #endif |
| 404 | }; |
| 405 | |
| 406 | template<> struct gemv_selector<OnTheRight,ColMajor,true> |
| 407 | { |
| 408 | template<typename ProductType, typename Dest> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 409 | static inline void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 410 | { |
| 411 | typedef typename ProductType::Index Index; |
| 412 | typedef typename ProductType::LhsScalar LhsScalar; |
| 413 | typedef typename ProductType::RhsScalar RhsScalar; |
| 414 | typedef typename ProductType::Scalar ResScalar; |
| 415 | typedef typename ProductType::RealScalar RealScalar; |
| 416 | typedef typename ProductType::ActualLhsType ActualLhsType; |
| 417 | typedef typename ProductType::ActualRhsType ActualRhsType; |
| 418 | typedef typename ProductType::LhsBlasTraits LhsBlasTraits; |
| 419 | typedef typename ProductType::RhsBlasTraits RhsBlasTraits; |
| 420 | typedef Map<Matrix<ResScalar,Dynamic,1>, Aligned> MappedDest; |
| 421 | |
| 422 | ActualLhsType actualLhs = LhsBlasTraits::extract(prod.lhs()); |
| 423 | ActualRhsType actualRhs = RhsBlasTraits::extract(prod.rhs()); |
| 424 | |
| 425 | ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) |
| 426 | * RhsBlasTraits::extractScalarFactor(prod.rhs()); |
| 427 | |
| 428 | enum { |
| 429 | // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1 |
| 430 | // on, the other hand it is good for the cache to pack the vector anyways... |
| 431 | EvalToDestAtCompileTime = Dest::InnerStrideAtCompileTime==1, |
| 432 | ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex), |
| 433 | MightCannotUseDest = (Dest::InnerStrideAtCompileTime!=1) || ComplexByReal |
| 434 | }; |
| 435 | |
| 436 | gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest; |
| 437 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 438 | bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0)); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 439 | bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible; |
| 440 | |
| 441 | RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha); |
| 442 | |
| 443 | ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(), |
| 444 | evalToDest ? dest.data() : static_dest.data()); |
| 445 | |
| 446 | if(!evalToDest) |
| 447 | { |
| 448 | #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN |
| 449 | int size = dest.size(); |
| 450 | EIGEN_DENSE_STORAGE_CTOR_PLUGIN |
| 451 | #endif |
| 452 | if(!alphaIsCompatible) |
| 453 | { |
| 454 | MappedDest(actualDestPtr, dest.size()).setZero(); |
| 455 | compatibleAlpha = RhsScalar(1); |
| 456 | } |
| 457 | else |
| 458 | MappedDest(actualDestPtr, dest.size()) = dest; |
| 459 | } |
| 460 | |
| 461 | general_matrix_vector_product |
| 462 | <Index,LhsScalar,ColMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run( |
| 463 | actualLhs.rows(), actualLhs.cols(), |
| 464 | actualLhs.data(), actualLhs.outerStride(), |
| 465 | actualRhs.data(), actualRhs.innerStride(), |
| 466 | actualDestPtr, 1, |
| 467 | compatibleAlpha); |
| 468 | |
| 469 | if (!evalToDest) |
| 470 | { |
| 471 | if(!alphaIsCompatible) |
| 472 | dest += actualAlpha * MappedDest(actualDestPtr, dest.size()); |
| 473 | else |
| 474 | dest = MappedDest(actualDestPtr, dest.size()); |
| 475 | } |
| 476 | } |
| 477 | }; |
| 478 | |
| 479 | template<> struct gemv_selector<OnTheRight,RowMajor,true> |
| 480 | { |
| 481 | template<typename ProductType, typename Dest> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 482 | static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 483 | { |
| 484 | typedef typename ProductType::LhsScalar LhsScalar; |
| 485 | typedef typename ProductType::RhsScalar RhsScalar; |
| 486 | typedef typename ProductType::Scalar ResScalar; |
| 487 | typedef typename ProductType::Index Index; |
| 488 | typedef typename ProductType::ActualLhsType ActualLhsType; |
| 489 | typedef typename ProductType::ActualRhsType ActualRhsType; |
| 490 | typedef typename ProductType::_ActualRhsType _ActualRhsType; |
| 491 | typedef typename ProductType::LhsBlasTraits LhsBlasTraits; |
| 492 | typedef typename ProductType::RhsBlasTraits RhsBlasTraits; |
| 493 | |
| 494 | typename add_const<ActualLhsType>::type actualLhs = LhsBlasTraits::extract(prod.lhs()); |
| 495 | typename add_const<ActualRhsType>::type actualRhs = RhsBlasTraits::extract(prod.rhs()); |
| 496 | |
| 497 | ResScalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) |
| 498 | * RhsBlasTraits::extractScalarFactor(prod.rhs()); |
| 499 | |
| 500 | enum { |
| 501 | // FIXME find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1 |
| 502 | // on, the other hand it is good for the cache to pack the vector anyways... |
| 503 | DirectlyUseRhs = _ActualRhsType::InnerStrideAtCompileTime==1 |
| 504 | }; |
| 505 | |
| 506 | gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs; |
| 507 | |
| 508 | ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(), |
| 509 | DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data()); |
| 510 | |
| 511 | if(!DirectlyUseRhs) |
| 512 | { |
| 513 | #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN |
| 514 | int size = actualRhs.size(); |
| 515 | EIGEN_DENSE_STORAGE_CTOR_PLUGIN |
| 516 | #endif |
| 517 | Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs; |
| 518 | } |
| 519 | |
| 520 | general_matrix_vector_product |
| 521 | <Index,LhsScalar,RowMajor,LhsBlasTraits::NeedToConjugate,RhsScalar,RhsBlasTraits::NeedToConjugate>::run( |
| 522 | actualLhs.rows(), actualLhs.cols(), |
| 523 | actualLhs.data(), actualLhs.outerStride(), |
| 524 | actualRhsPtr, 1, |
| 525 | dest.data(), dest.innerStride(), |
| 526 | actualAlpha); |
| 527 | } |
| 528 | }; |
| 529 | |
| 530 | template<> struct gemv_selector<OnTheRight,ColMajor,false> |
| 531 | { |
| 532 | template<typename ProductType, typename Dest> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 533 | static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 534 | { |
| 535 | typedef typename Dest::Index Index; |
| 536 | // TODO makes sure dest is sequentially stored in memory, otherwise use a temp |
| 537 | const Index size = prod.rhs().rows(); |
| 538 | for(Index k=0; k<size; ++k) |
| 539 | dest += (alpha*prod.rhs().coeff(k)) * prod.lhs().col(k); |
| 540 | } |
| 541 | }; |
| 542 | |
| 543 | template<> struct gemv_selector<OnTheRight,RowMajor,false> |
| 544 | { |
| 545 | template<typename ProductType, typename Dest> |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 546 | static void run(const ProductType& prod, Dest& dest, const typename ProductType::Scalar& alpha) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 547 | { |
| 548 | typedef typename Dest::Index Index; |
| 549 | // TODO makes sure rhs is sequentially stored in memory, otherwise use a temp |
| 550 | const Index rows = prod.rows(); |
| 551 | for(Index i=0; i<rows; ++i) |
| 552 | dest.coeffRef(i) += alpha * (prod.lhs().row(i).cwiseProduct(prod.rhs().transpose())).sum(); |
| 553 | } |
| 554 | }; |
| 555 | |
| 556 | } // end namespace internal |
| 557 | |
| 558 | /*************************************************************************** |
| 559 | * Implementation of matrix base methods |
| 560 | ***************************************************************************/ |
| 561 | |
| 562 | /** \returns the matrix product of \c *this and \a other. |
| 563 | * |
| 564 | * \note If instead of the matrix product you want the coefficient-wise product, see Cwise::operator*(). |
| 565 | * |
| 566 | * \sa lazyProduct(), operator*=(const MatrixBase&), Cwise::operator*() |
| 567 | */ |
| 568 | template<typename Derived> |
| 569 | template<typename OtherDerived> |
| 570 | inline const typename ProductReturnType<Derived, OtherDerived>::Type |
| 571 | MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const |
| 572 | { |
| 573 | // A note regarding the function declaration: In MSVC, this function will sometimes |
| 574 | // not be inlined since DenseStorage is an unwindable object for dynamic |
| 575 | // matrices and product types are holding a member to store the result. |
| 576 | // Thus it does not help tagging this function with EIGEN_STRONG_INLINE. |
| 577 | enum { |
| 578 | ProductIsValid = Derived::ColsAtCompileTime==Dynamic |
| 579 | || OtherDerived::RowsAtCompileTime==Dynamic |
| 580 | || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), |
| 581 | AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, |
| 582 | SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) |
| 583 | }; |
| 584 | // note to the lost user: |
| 585 | // * for a dot product use: v1.dot(v2) |
| 586 | // * for a coeff-wise product use: v1.cwiseProduct(v2) |
| 587 | EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), |
| 588 | INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) |
| 589 | EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), |
| 590 | INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) |
| 591 | EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) |
| 592 | #ifdef EIGEN_DEBUG_PRODUCT |
| 593 | internal::product_type<Derived,OtherDerived>::debug(); |
| 594 | #endif |
| 595 | return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); |
| 596 | } |
| 597 | |
| 598 | /** \returns an expression of the matrix product of \c *this and \a other without implicit evaluation. |
| 599 | * |
| 600 | * The returned product will behave like any other expressions: the coefficients of the product will be |
| 601 | * computed once at a time as requested. This might be useful in some extremely rare cases when only |
| 602 | * a small and no coherent fraction of the result's coefficients have to be computed. |
| 603 | * |
| 604 | * \warning This version of the matrix product can be much much slower. So use it only if you know |
| 605 | * what you are doing and that you measured a true speed improvement. |
| 606 | * |
| 607 | * \sa operator*(const MatrixBase&) |
| 608 | */ |
| 609 | template<typename Derived> |
| 610 | template<typename OtherDerived> |
| 611 | const typename LazyProductReturnType<Derived,OtherDerived>::Type |
| 612 | MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived> &other) const |
| 613 | { |
| 614 | enum { |
| 615 | ProductIsValid = Derived::ColsAtCompileTime==Dynamic |
| 616 | || OtherDerived::RowsAtCompileTime==Dynamic |
| 617 | || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime), |
| 618 | AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime, |
| 619 | SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived) |
| 620 | }; |
| 621 | // note to the lost user: |
| 622 | // * for a dot product use: v1.dot(v2) |
| 623 | // * for a coeff-wise product use: v1.cwiseProduct(v2) |
| 624 | EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes), |
| 625 | INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS) |
| 626 | EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors), |
| 627 | INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION) |
| 628 | EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT) |
| 629 | |
| 630 | return typename LazyProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); |
| 631 | } |
| 632 | |
| 633 | } // end namespace Eigen |
| 634 | |
| 635 | #endif // EIGEN_PRODUCT_H |